From 07a78bf335b0b5039a1fecd7cd74a0fe1d1605ea Mon Sep 17 00:00:00 2001
From: Ali Romman <226877337+aliromman@users.noreply.github.com>
Date: Sun, 11 Oct 2026 12:13:36 -0500
Subject: [PATCH] Prove the nine-call lower bound for three abelian bit
 registers

---
 .../bit-three-register-lower-bound/README.md  | 632 ++++++++++++++++++
 .../catalogue.py                              | 185 +++++
 .../controls.py                               |  89 +++
 .../expected.json                             |  23 +
 .../bit-three-register-lower-bound/f101.cpp   | 449 +++++++++++++
 .../rational.cpp                              | 403 +++++++++++
 .../bit-three-register-lower-bound/replay.py  |  98 +++
 7 files changed, 1879 insertions(+)
 create mode 100644 research/bit-three-register-lower-bound/README.md
 create mode 100644 research/bit-three-register-lower-bound/catalogue.py
 create mode 100644 research/bit-three-register-lower-bound/controls.py
 create mode 100644 research/bit-three-register-lower-bound/expected.json
 create mode 100644 research/bit-three-register-lower-bound/f101.cpp
 create mode 100644 research/bit-three-register-lower-bound/rational.cpp
 create mode 100644 research/bit-three-register-lower-bound/replay.py

diff --git a/research/bit-three-register-lower-bound/README.md b/research/bit-three-register-lower-bound/README.md
new file mode 100644
index 0000000000..25e41c370c
--- /dev/null
+++ b/research/bit-three-register-lower-bound/README.md
@@ -0,0 +1,632 @@
+# Nine calls are necessary in the abelian three-register bit model
+
+Eight rational rank-one calls cannot produce three dimension-three bit
+interchanges in the model below. Nine calls are necessary and sufficient.
+The proof combines algebraic reductions with a complete finite calculation.
+It checks 2,276 scalar multiplicity families and all 21,317,520 distinct
+orders. An adaptive exact search excludes every rational frame relation
+space. There is no bound on frame numerators or denominators.
+
+This is a computer-assisted theorem for three registers. It does not cover
+extra registers or arbitrary free address maps between paid calls. It gives
+no improvement in $\kappa$ or the integer-multiplication exponent. It is not
+a proof-assistant verification.
+
+## 1. Model, attribution, and statement
+
+The additive group $M_3(\mathbb Q)$ has formal symbols $[A]$ with
+$[A][B]=[A+B]$. Work in its group algebra over $\mathbb F_2$. A paid call is
+
+$$
+D_N([A])=I+([A]-1)N,\qquad
+\operatorname{rank}_{\mathbb Q}A=1,\qquad
+N=xy^T,\quad x,y\in\mathbb F_2^3,\quad y^Tx=1.
+$$
+
+There are 28 possible binary rank-one idempotents $N$. We ask whether
+
+$$
+\prod_{j=1}^{8}D_{N_j}([A_j])
+ =\operatorname{diag}([C_1],[C_2],[C_3]),
+\qquad C_w\in M_3(\mathbb Q)
+\text{ invertible and lower triangular}. \tag{1}
+$$
+
+Factor order matters. Rational frame entries are not reduced modulo two.
+Equation (1) requires equality of every formal coefficient. An identity at
+one finite address size does not suffice.
+
+This is the scalar-idempotent abelian register model in
+[eumemic's PR #288, section 1.2, pinned at
+`475a088b5862f75f0b10596cfd6842d19531a5ec`](https://github.com/eumemic/integer-mult-bounds/blob/475a088b5862f75f0b10596cfd6842d19531a5ec/research/negative-results-paradigm-search/README.md).
+That source also gives the dimension-two flag bound used below. Free
+scalar XOR maps are represented by the choices of $N_j$. A final output
+permutation can be removed by relabeling roles. In the normalized identity,
+augmentation $[A]\mapsto1$ forces the remaining scalar output map to be
+the identity. Arbitrary free address maps between calls are outside this
+model.
+
+**Theorem.** For three registers over $\mathbb F_2$, nonzero rational
+rank-one paid frames, binary rank-one scalar idempotents, and three
+invertible lower-triangular rational frame targets, the minimum number of
+paid factors is nine.
+
+The pinned PR #288, section 1.2(f), leaves the three-register cases with
+seven or eight calls open after symbolic bounds and bounded searches.
+This note excludes both cases for the stated model, with no bound on the
+rational frames. A related earlier result is
+[huxint's PR #203](https://github.com/CrocSwap/integer-mult-bounds/pull/203),
+at [commit `d2873e20750d0949d2ad75cbc4cd82b5b1a74e84`](https://github.com/huxint/integer-mult-bounds/commit/d2873e20750d0949d2ad75cbc4cd82b5b1a74e84):
+paid frames with a common triangular flag cannot save rank. Our paid
+frames need not have such a flag. We use the flag and projection methods
+of these predecessors in the reductions below and do not claim those
+methods as new. This comparison concerns the pinned sources; no broader
+claim of literature priority is made. The proof, implementations, and
+checks were prepared with substantial OpenAI Codex assistance.
+
+## 2. Projection and rank tools
+
+Restrict each group-algebra argument to the subgroup generated by its
+finitely many exponents. This group is free abelian. Its group algebra is
+a Laurent polynomial domain and has a fraction field.
+
+For $k$ nontrivial factors, write $N_j=x_jy_j^T$. Over that fraction field,
+the ordered product has the form
+
+$$
+\prod_{j=1}^k(I+a_jx_jy_j^T)=I+XTY^T,
+\qquad X=(x_1\ \cdots\ x_k),\quad Y=(y_1\ \cdots\ y_k),
+$$
+
+where $T$ is upper triangular with nonzero diagonal $a_1,\ldots,a_k$.
+If the product is $I$, the rank inequality gives
+
+$$
+\operatorname{rank}X+\operatorname{rank}Y\le k. \tag{2}
+$$
+
+If the difference from $I$ instead has rank $W$, both active scalar spans
+have dimension at least $W$. Rational linear projections of the frames
+give group-algebra homomorphisms, so these bounds apply after projection.
+
+## 3. Frame span at most four
+
+**Low-span lemma.** In address dimension $m$ with $W$ registers, if the
+paid rank-one frames have rational span dimension at most $m$, invertible
+targets require at least $mW$ calls.
+
+Choose a basis $B_i=u_iv_i^T$ of that span from the paid frames. Every
+target belongs to the span because it occurs in the expanded product's
+support. A span dimension below $m$ cannot contain a rank-$m$ matrix. At
+dimension $m$, write one invertible target as
+$C=\sum_i c_iB_i=X\operatorname{diag}(c_i)Y^T$. Both $X,Y$ are invertible
+and each $c_i$ is nonzero. In these left and right bases, the frame span
+is diagonal. Each paid rank-one frame is on one basis line. Each target
+has a nonzero coefficient on every line. Projection onto one coefficient
+gives target-minus-identity rank $W$, so that line needs at least $W$
+calls. The $m$ sets of calls are disjoint.
+
+**Span-four lemma.** In dimension three, a paid frame span of dimension at
+most four requires at least $3W$ calls for $W$ invertible lower-triangular
+targets.
+
+Let $S$ be the frame span and choose its basis from paid rank-one matrices.
+The determinant of a linear combination is a nonzero homogeneous
+polynomial of degree three. It has a nonzero mixed coefficient at three
+distinct basis matrices: a repeated rank-one summand gives zero in the
+determinant expansion. The left vectors of those three matrices are
+independent, as are their right vectors. Common invertible rational row
+and column maps make them $E_{00},E_{11},E_{22}$. The low-span lemma handles
+dimension three. At dimension four, the last basis matrix is $B=uv^T$,
+and every matrix in $S$ has the form $D+tB$, with $D$ diagonal.
+
+Suppose all coordinates of $u,v$ are nonzero. For distinct $i,k,l$, the
+minor of $D+tB$ on rows $i,k$ and columns $i,l$ is $D_{ii}t u_kv_l$.
+Thus a rank-one matrix with $t\ne0$ has $D=0$. All paid frames lie on the
+four independent lines $E_{00},E_{11},E_{22},B$. Project onto one line's
+coefficient. Target-minus-identity rank equals the number of targets with
+nonzero coefficient there, and each retained call contributes at most one.
+Every invertible target needs at least three nonzero coefficients, by rank
+subadditivity. The sum of the four line bounds is therefore at least $3W$.
+
+Otherwise $u$ or $v$ has a zero coordinate. In the directed graph of $B$,
+each vertex of a cycle lies in
+$\operatorname{supp}(u)\cap\operatorname{supp}(v)$. Each strongly connected
+component consequently has size at most two. A simultaneous row and column
+permutation makes $S$ block upper triangular with diagonal blocks of size
+at most two. A rank-one matrix in this space has at most one nonzero
+diagonal block. Two such blocks would force a nonzero block below the
+diagonal. Assign each paid call to its nonzero diagonal block, if any.
+
+Every target has invertible diagonal blocks. The original ratios
+$C_1^{-1}C_w$ are simultaneously lower triangular. After the common row
+and column maps, the ratios are conjugate to these matrices and are block
+upper triangular. Each diagonal block is an invariant subquotient of a
+family with a complete common rational flag. Intersect that flag with the
+invariant subspaces, pass to quotients, and remove repeated spaces. Each
+block retains a complete common flag. Common invertible row and column
+maps within that block then make all its targets lower triangular.
+
+Here is the dimension-two flag bound. Let $E_2$ count calls assigned to
+one two-dimensional block, and let $U_2$ count its frames with nonzero
+upper entry. Write $r_X,r_Y$ for their scalar ranks. Upper-entry projection
+has zero target and gives $r_X+r_Y\le U_2$ by (2). For diagonal coordinate
+$i$, let $l_i$ count lower frames with nonzero entry there. Diagonal
+projection gives $W\le r_X+l_i$ and $W\le r_Y+l_i$. A lower rank-one
+frame has at most one nonzero diagonal entry. Adding the four inequalities
+gives
+
+$$
+4W\le2(r_X+r_Y+l_0+l_1)\le2E_2.
+$$
+
+A one-dimensional block needs $W$ calls by target rank. Block dimensions
+sum to three and assignments are disjoint, which proves the lemma. This
+step uses the lower-triangular target hypothesis.
+
+## 4. Frame span at least seven and scalar irreducibility
+
+**At-most-one-relation lemma.** Eight paid frames of span dimension seven
+or eight cannot give three invertible targets. Lower-triangular targets
+are not needed for this lemma.
+
+Introduce independent invertible variables $t_j$ and set
+$P(t)=\prod_j(I+(t_j-1)N_j)$. Each entry is multilinear. For independent
+frames, all subset sums differ. Every off-diagonal coefficient of $P$
+must therefore vanish. Differentiation at $t=(1,\ldots,1)$ makes each
+$N_j$ diagonal. Each paid frame is assigned to one target, and rank
+subadditivity gives the contradiction $9\le8$.
+
+At span seven, the integer relation lattice has a primitive generator
+$r\in\mathbb Z^8$. If some $|r_j|\ge2$, subset sums still cannot collide.
+Otherwise $r_j\in\{-1,0,1\}$. Set $S=\{j:r_j\ne0\}$. The substitution
+$t_j\mapsto[A_j]$ has Laurent ideal kernel generated by
+
+$$
+g=\prod_{r_j=1}t_j-\prod_{r_j=-1}t_j.
+$$
+
+A primitive vector extends to an integer basis, which proves this kernel
+description. Each off-diagonal entry of $P$ is a multiple of $g$.
+Differentiation at one shows that all projectors on $S$ have one common
+off-diagonal pattern; all others are diagonal.
+
+For $N=xy^T$ with $y^Tx=1$, the supports of $x,y$ meet in one or three
+coordinates. In the latter case $N=J$, the all-ones matrix. In the former,
+its directed off-diagonal graph has no cycle, since each cycle vertex
+would lie in both supports. Unless the common pattern is that of $J$,
+one coordinate order makes every projector triangular. Each projector
+has exactly one nonzero diagonal entry, so the disjoint rank bound applies.
+
+For the pattern of $J$, all projectors on $S$ equal $J$. Set outside
+variables to one; the sole relation permits this. The remaining product
+is $I+([\sum_{j\in S}A_j]-1)J$. Diagonality forces this sum to vanish.
+The primitive relation thus has all plus signs on $S$, up to overall sign.
+Keep one outside factor $E_{kk}$, and let $B$ be the sum of the support
+frames before it. With $z=[B]$, the conjugated projector is
+
+$$
+(I+(z-1)J)E_{kk}(I+(z^{-1}-1)J).
+$$
+
+For distinct $a,b\ne k$, its $(a,b)$ entry is $z+z^{-1}$ and must vanish.
+The rational address group is torsion-free, so $B=0$. The only zero-sum
+subsets of $S$ are empty and full. Every outside factor is before or after
+all support factors. The support forms an identity block, and the remaining
+coordinate factors give the same rank contradiction.
+
+**Two-register lemma.** Two invertible rational $3\times3$ targets require
+at least six calls over $\mathbb F_2$.
+
+For at most five calls, span at most three is excluded by the low-span
+lemma. Higher span has at most one relation. The preceding Laurent
+argument forces diagonal projectors, or one common off-diagonal pattern
+on the relation support. A binary $2\times2$ rank-one idempotent cannot
+have both off-diagonal entries equal to one: rank one would force both
+diagonal entries to one, but that all-ones matrix squares to zero. Thus all
+projectors are triangular in one order. Each frame belongs to one target,
+whose ranks sum to six.
+
+**Irreducibility lemma.** The projectors in an eight-call solution have no
+common proper nonzero binary invariant subspace.
+
+Suppose they preserve $V$. The diagonal target preserves $V$ too. Group
+coordinate axes by equal $C_w$. Each distinct target monomial has the
+coordinate projection onto its group as coefficient. Independence of
+monomials shows that these projections preserve $V$. Choose bases within
+the groups adapted to their intersections with $V$. The target stays
+diagonal, and the projectors become block triangular with block sizes one
+and two. For a rank-one idempotent, exactly one diagonal block is a
+rank-one idempotent and the other is zero; this follows from $N=xy^T$ and
+$y^Tx=1$. The one-register block needs three calls by rank subadditivity.
+The two-register block needs six. These assignments are disjoint.
+
+## 5. Frame span six
+
+### 5.1. Common-$J$ calculation
+
+Consider projectors in $\{J,E_{00},E_{11},E_{22}\}$ with at most two
+independent rational frame relations. Suppose the output is diagonal or
+has only one possible off-diagonal entry, at $01$. We show that the $J$
+frames make zero net contribution to the target diagonals.
+
+Put $y_j=t_j-1$. For the generated free abelian group $H$, the associated
+graded algebra at the augmentation ideal is the polynomial algebra on
+$H/2H$. The integer relation lattice is saturated, so the kernel from
+$\mathbb F_2^E$ has the rational relation dimension. The degree-one
+off-diagonal term is $J_{\rm off}L$, where $L=\sum_{N_j=J}y_j$. If $J$
+occurs, the output condition forces the nonzero form $L$ to vanish on
+$H/2H$.
+
+Set $K_i=\operatorname{offdiag}(JE_{ii}+E_{ii}J)$. The matrices $K_0,K_1$
+are independent modulo $\operatorname{span}(J_{\rm off},E_{01})$:
+$K_0,K_1,J_{\rm off}$ are symmetric and independent, whereas $E_{01}$ is
+not symmetric. Define
+
+$$
+q_i=\sum_{j:N_j=E_{ii}}y_j\sum_{k<j:N_k=J}y_k.
+$$
+
+After imposing $L=0$, the degree-two conditions are
+$q_0+q_1=q_0+q_2=0$. In this matrix quotient, a group relation gives no
+second-order correction because its first-order term is already zero.
+At most one additional linear relation remains. Thus the nonzero
+quadratics must have one common linear divisor.
+
+Eliminate the last $J$ variable by $L=0$. If there are $h$ such factors,
+their $h-1$ proper prefix sums $p_a$ are independent. Each displayed
+quadratic is bilinear between the remaining $J$ variables and the
+coordinate-factor variables. A nonzero bilinear polynomial with a linear
+divisor has rank one, with one factor in each variable set. Compare
+bidegrees $(2,0)$ and $(0,2)$ to see this. Its nonzero columns here are
+prefix vectors $p_a$, so all contributing coordinate factors occupy one
+internal gap between $J$ factors.
+
+If both quadratics are nonzero, an internal $E_{00}$ makes their gaps
+equal. Without one, they use disjoint nonempty sets of $E_{11}$ and
+$E_{22}$ variables; different gaps would leave no common linear divisor.
+If one quadratic is zero, the other gives the same single-gap conclusion.
+If both vanish, there is no internal coordinate factor. There are therefore
+at most two $J$ runs separated by a diagonal block.
+
+Combine the runs and remove exterior diagonal factors. The product is
+$U(b)DU(c)$, where $U(z)=I+(z-1)J$,
+$D=\operatorname{diag}(d_0,d_1,d_2)$, and $b,c,d_i$ are monomials. Its
+off-diagonal entry $ij$ is
+
+$$
+(b-1)d_j+(c-1)d_i+(b-1)(c-1)(d_0+d_1+d_2).
+$$
+
+If $b=1$ or $c=1$, the required zero entries force both to be one.
+Otherwise entries $10,20$ give $d_1=d_2$, and entries $02,12$ give
+$d_0=d_1$. Any required zero entry then gives $bc=1$. The output equals
+the product of the coordinate factors. The possible $01$ entry also
+vanishes. Target exponents are disjoint sums of coordinate frames.
+
+### 5.2. Two relations and the finite classification
+
+At span six the saturated integer relation lattice has rank two. Call a
+relation a *cube relation* if all its entries belong to $\{-1,0,1\}$.
+Cube relations span the full rational relation space. Otherwise quotient
+$\mathbb Z^8$ by their saturated span, of dimension at most one. This lift
+preserves every equality of subset sums, hence every coefficient of the
+product. Lift each target by a subset with its nonzero coefficient. The
+algebraic at-most-one-relation argument assigns disjoint frame sets to the
+lifted targets, possibly after canceling a $J$ block. Map back to the
+actual rank-one frames to obtain $9\le8$.
+
+Choose independent cube relations $r,s$. If their reductions modulo two
+are independent, some minor is odd. A $2\times2$ cube minor has absolute
+value at most two, so that minor is $\pm1$ and the basis is saturated.
+Otherwise the reductions agree. Replace $r,s$ by $(r+s)/2,(r-s)/2$.
+These are independent cube vectors with disjoint nonempty supports; again
+they have a unit minor. Thus a saturated cube-valued basis exists.
+
+First-order off-diagonal coefficients have the form
+
+$$
+\operatorname{offdiag}(N_j)=(r_j\bmod2)B+(s_j\bmod2)C.
+$$
+
+Pattern dimension zero gives coordinate projectors. Dimension one gives
+an acyclic common pattern or $J$; triangularity or section 5.1 excludes
+it. At dimension two, $B,C$ are independent. A zero pattern then means
+$r_j=s_j=0$ over the integers. Let $S$ index nonzero patterns and $T$ the
+remaining coordinate projectors. Both complete relations are supported
+on $S$, so the actual group splits as $H=H_S\oplus\mathbb Z^T$.
+
+The finite classification enumerates all 512 binary matrices, selects the
+28 rank-one idempotents and their 19 nonzero patterns, and forms all 129
+two-dimensional pattern spaces. For each nonzero-pattern pool it checks
+all 21 complete binary flags. There are 123 pools with a common flag and
+six exceptions. The exceptions are exactly the coordinate permutations of
+
+$$
+P=\begin{pmatrix}1&1&0\\0&0&0\\0&0&0\end{pmatrix},\qquad
+Q=\begin{pmatrix}0&1&0\\0&1&0\\0&0&0\end{pmatrix},\qquad
+R=\begin{pmatrix}0&0&0\\1&0&1\\1&0&1\end{pmatrix}.
+$$
+
+Every proper subset of an exception has a common flag. Two independent
+occurring patterns generate every dimension-two case, so this enumeration
+is exhaustive. A separate implementation with tuple matrices and actual
+vector subspaces reproduces the classification.
+
+For a common-flag pool, project onto $H_S$ and use its scalar flag basis.
+Each rank-one idempotent has exactly one nonzero diagonal entry. Product
+diagonals therefore assign disjoint sets of $S$ frames. Conjugation
+preserves the characteristic polynomial. Unique factorization over the
+fraction field identifies these diagonal monomials, with multiplicity,
+with the three original projected target monomials. Their rational matrix
+ranks sum to at most $|S|$.
+
+For an exceptional pool, put
+
+$$
+S_0=\begin{pmatrix}1&1&0\\0&1&0\\0&0&1\end{pmatrix}.
+$$
+
+Direct multiplication gives $S_0PS_0=E_{00}$, $S_0QS_0=E_{11}$, and
+$S_0RS_0=J$. The conjugated target has only a possible $01$ off-diagonal
+entry. Section 5.1 assigns disjoint $S$ frame sums to its diagonals. The
+same characteristic-polynomial argument gives the rank bound $|S|$.
+
+Project onto each free $T$ generator. Only its coordinate-projector factor
+remains, so that frame belongs to one original target. With
+$C_w=C_w^S+C_w^T$, rank subadditivity now gives
+
+$$
+9=\sum_w\operatorname{rank}C_w
+\le\sum_w\operatorname{rank}C_w^S+\sum_w\operatorname{rank}C_w^T
+\le |S|+|T|=8.
+$$
+
+This excludes span six. Sections 3–5 leave exactly span five.
+
+## 6. Scalar distance
+
+For each nonzero binary row $r$, at least three indices have $rx_j\ne0$.
+For each nonzero binary column $v$, at least three have $y_j^Tv\ne0$.
+
+Suppose at most two indices meet a row $r$, and let $V$ be the rational
+span of their frames. Each matrix in $V$ has rank at most two, so none of
+the invertible targets lies in $V$. Choose a rational linear functional
+$\ell$ that kills $V$ and is nonzero on every target. Such a functional
+exists: a finite union of proper subspaces cannot cover the annihilator
+of $V$ over $\mathbb Q$.
+
+Under this projection the selected factors become $I$. Each other factor
+fixes $r$ on the left because $rN_j=0$. Their product fixes $r$, but the
+projected diagonal target cannot: all its diagonal entries differ from
+one. This contradiction proves the row claim. Transpose and reverse the
+word to obtain the column claim.
+
+## 7. Complete scalar catalogue for span five
+
+Let $L\subset\mathbb Z^8$ be the saturated relation lattice, of rank
+three. Cube relations span $L$ rationally. Otherwise let $L_0$ be their
+rational span intersected with $\mathbb Z^8$; its rank is at most two.
+Lift the word to $\mathbb Z^8/L_0$. Every subset equality, and thus every
+coefficient, is preserved. The low-corank arguments of sections 4–5 assign
+the lifted target exponents to disjoint frame sets, possibly after
+canceling common-$J$ factors. Mapping back gives total target rank at
+most eight, a contradiction.
+
+Three independent cube vectors can be chosen whose reductions form a
+basis of $L/2L$. Start with any independent cube triple. If its integer
+span has even index in $L$, a nonempty subset has an even coordinate sum.
+It has two or three members: a nonzero cube vector is not itself even.
+Half the sum is an integral cube vector in $L$. Replace one member by that
+half-sum. Independence is retained and the index is halved. Repetition
+ends at odd index. Saturation embeds $L/2L$ into $\mathbb F_2^8$, so the
+three reductions are independent. This is a parity basis, not necessarily
+an integer basis of $L$.
+
+Writing the cube vectors as $r_1,r_2,r_3$, the first-order condition is
+
+$$
+\operatorname{offdiag}(N_j)=\sum_{a=1}^{3}(r_{aj}\bmod2)B_a.
+$$
+
+Thus the off-diagonal pattern span has dimension at most three. If its
+dimension is three, the $B_a$ are independent. At each zero-pattern index,
+all three cube entries are zero over $\mathbb Z$. Every relation therefore
+vanishes at these indices. With $S$ the noncoordinate indices and $T$ the
+coordinate indices, $H=H_S\oplus\mathbb Z^T$.
+
+In this dimension-three case, the noncoordinate projectors cannot have a
+common complete binary flag. Such a flag, followed by projection onto
+$H_S$, would assign disjoint $S$ frame sums to the diagonals. The
+characteristic-polynomial argument identifies them with the projected
+target monomials. Each free $T$ frame belongs to one original target, so
+the same rank budget would give $9\le8$. **This filter is used only at
+off-diagonal dimension three.**
+
+Enumerate the 28 projectors and 19 nonzero patterns. Closing subsets of
+these patterns under addition gives $1,19,129,363$ possible pattern spaces
+of dimensions $0,1,2,3$. Coordinate permutations and transpose reduce
+their projector pools to 67 classes. Coordinate permutations relabel
+target registers. Transpose reverses factor order and preserves a diagonal
+target; it is valid because every order is checked. Arbitrary changes of
+scalar basis are not used to identify full target problems.
+
+For each pool, enumerate all eight-factor multiplicity vectors whose
+support has that pool's actual off-diagonal span. Retain irreducible
+supports, row and column distance at least three, and the dimension-three
+common-flag condition. Identify multiplicities only under the pool's
+stabilizer. This leaves 2,276 records: 4 in pattern dimension one, 137 in
+dimension two, and 2,135 in dimension three. No additional subquotient or
+noncoordinate-count filter is used.
+
+The included `catalogue.py` selects projectors from all 512 binary matrices.
+It uses tuple matrices, explicit invariant subspaces, and sorted multisets.
+Its catalogue was checked against a separate implementation with
+bit-matrix rank and integer-composition routines. Every record agrees,
+not only the count; the replay checks the regenerated catalogue digest.
+Multiplicities $m_1,\ldots,m_k$ have
+$8!/\prod_i m_i!$ distinct orders. Their total is 21,317,520. This is a
+complete catalogue of scalar cases after proved reductions, not a bounded
+search over physical frames.
+
+## 8. Exhaust all rational frame relation spaces
+
+Fix one ordered scalar word. Expand
+
+$$
+P(t)=\prod_{j=1}^{8}\bigl((I+N_j)+t_jN_j\bigr)
+=\sum_{a\in\{0,1\}^8}M_a t^a.
+$$
+
+There are at most 256 subset exponents with ordered binary matrix
+coefficients. The unknown physical relation space is a rational subspace
+$R\subset\mathbb Q^8$ of dimension three. A search state is
+$R_0\subseteq R$, initially zero, generated by at most three cube vectors.
+Group subset indicators $a,b$ exactly when $a-b\in R_0$ and add their
+matrix coefficients over $\mathbb F_2$. This gives the image in the group
+algebra of $\mathbb Z^8/(R_0\cap\mathbb Z^8)$.
+
+If a current exponent group has a nonzero off-diagonal bit, any valid
+diagonal completion must merge it with another group containing that bit.
+Branch over every such partner and add the difference of their subset
+representatives to $R_0$. At least one branch follows every valid
+completion. Distinct groups have difference outside $R_0$, so each branch
+raises its dimension. The depth is at most three. Equal canonical relation
+spaces need be checked only once per word.
+
+Reject a branch if $e_j\in R_0$, since then $A_j=0$. Also reject it when
+the binary span of its primitive integral relation rows and the
+first-order off-diagonal coefficient rows has dimension greater than
+three. Each such row lies in $L/2L$: saturation makes every integral vector
+of $R$ a member of $L$. This test reduces integral relation vectors,
+not rational physical frame entries. At dimension three, a remaining
+off-diagonal coefficient rejects the branch.
+
+When all off-diagonal coefficients vanish, each diagonal entry is a single
+monomial. Every factor has monomial determinant, so the diagonal entries
+multiply to a unit. Laurent polynomial units over $\mathbb F_2$ are
+exactly group monomials. Choose a subset representative $c_w$ of each
+diagonal target. For every pair of coordinates $i,j$, test whether
+
+$$
+c_w\in R_0+\operatorname{span}_{\mathbb Q}(e_i,e_j). \tag{3}
+$$
+
+If (3) holds for any target, its physical matrix lies in the span of two
+paid rank-one frames and has rank at most two. Additional relations cannot
+remove this obstruction. Reject that terminal. A terminal that passed
+this check would only be a candidate, not a construction.
+
+The rational engine uses primitive integer echelon rows. Their entries
+are signed minors of at most three cube rows, divided by a common divisor,
+so their absolute values are at most four. The common denominator in a
+quotient key divides $\operatorname{lcm}(1,2,3,4)=12$. The terminal test
+adds two coordinate rows and then one subset row. After reduction, the
+first added coordinate row has entries at most four in absolute value;
+the second has entries at most 32. Reducing the subset row by the three
+original rows gives a coarse bound $8^3$. Reduction by the two added rows
+gives $8^3\cdot8\cdot64=262{,}144$. Normalization cannot increase these
+bounds. These bounds also cover search-row insertion and keep all products
+and differences within signed 32-bit range. No floating-point arithmetic
+is used.
+
+The rational run checks all 21,317,520 words. It visits 94,392,894 relation
+states and reaches 415,398 diagonal terminals. Every terminal has the
+rank-at-most-two obstruction. There are zero candidates. The branching
+argument therefore excludes every admissible rational relation space for
+every retained scalar word.
+
+## 9. A second full implementation with exact finite-field arithmetic
+
+The independent full implementation uses elimination over
+$\mathbb F_{101}$. In these tests this is an exact representation of the
+rational calculation, not a probabilistic evaluation.
+
+Each tested matrix has at most three original cube rows, possibly one
+further cube or binary target row, and at most two coordinate rows. Expand
+a minor along the coordinate rows. At most four rows with entries of
+absolute value at most one remain. Hadamard's bound makes the determinant
+at most 16 in absolute value. Thus it vanishes over $\mathbb Q$ exactly
+when it vanishes modulo 101. This proves exact ranks and membership for
+subset collisions, zero frames, and two-frame target spans.
+
+Parity pruning first recovers the rational reduced-echelon rows. Cramer's
+rule bounds each entry, in lowest terms, by $a/b$ with $|a|\le4$ and
+$1\le b\le4$. For three cube rows, the determinant bound four follows
+by maximizing the multilinear determinant at sign matrices: their
+determinants are divisible by four and bounded by $3\sqrt3$. Smaller
+minors obey the same bound. Permitted fractions have distinct residues
+modulo 101, since $|ad-bc|\le32<101$. Clear the recovered denominators
+and divide by the row gcd to obtain its primitive integral parity.
+Taking the parity of a finite-field residue itself would be invalid.
+
+The two full implementations each check 2,276 records and 21,317,520
+orders. Each reaches 415,398 diagonal terminals and zero survivors.
+Every per-record word and terminal count agrees. The finite-field engine
+reports no nonmonomial-diagonal anomaly. It has separate matrix expansion,
+quotient grouping, and arithmetic representation.
+
+## 10. Minimum nine and verification scope
+
+Span five was the last possible case, so eight calls are impossible.
+Fewer calls are excluded by splitting any nonzero rank-one call:
+
+$$
+D_N([A])=D_N([A/2])D_N([A/2]).
+$$
+
+Both new frames remain nonzero and rank one. Repeated splitting would
+turn any shorter positive-length circuit into an eight-call circuit.
+Zero calls cannot give invertible frame targets.
+
+Nine calls suffice for every allowed target triple. For each $w$, write
+
+$$
+C_w=\sum_{i=0}^{2}(C_we_i)e_i^T
+$$
+
+and use these three nonzero rank-one frames with scalar projector $E_{ww}$.
+They combine to the $w$-th target and fix the other registers. This proves
+the theorem.
+
+The proof depends on the symbolic reductions, complete catalogue, exact
+collision branching, and completed replay. Neither rational frames nor
+their coefficient sizes are sampled. The lower-triangular target condition
+is used in the span-four reduction. The theorem makes no claim for extra
+registers, other paid primitives, or arbitrary free address maps between
+calls, and it does not improve the multiplication exponent.
+
+## 11. Replay
+
+From this directory, run:
+
+```sh
+python3 replay.py
+```
+
+To retain the generated input and per-case results, run:
+
+```sh
+python3 replay.py --work-dir /tmp/bit-three-register-replay
+```
+
+The replay requires Python 3.10 or later and Clang or GCC with C++17 and
+at least 32-bit `int`. Python uses only the standard library. It
+regenerates the finite catalogue, compiles and runs both complete relation
+engines, checks the per-family results, and runs the included controls.
+The required result is 2,276 scalar families, 21,317,520 orders, 415,398
+diagonal terminals, and zero survivors in each full engine.
+
+The nine-call fixture checks $C_1=C_2=C_3=I$. A separate eight-call fixture
+has noncommuting scalar factors and rational frame span five, but rank-two
+targets. The Python controls check its exact formal product and these matrix
+ranks. Another Python control detects an invalid reduction of rational
+address entries modulo two. The finite-field search core separately accepts a constructed
+positive control and rejects its one-factor mutation, which forces a
+rank-two target. The positive core control deliberately uses a rank-two
+scalar projector. It checks the search core's acceptance path and is not
+an admissible construction in the theorem's model.
+
+Separate implementations and source review support this computer-assisted
+proof. They are not Lean verification or a formal verification of the full
+integer-multiplication theorem.
diff --git a/research/bit-three-register-lower-bound/catalogue.py b/research/bit-three-register-lower-bound/catalogue.py
new file mode 100644
index 0000000000..fc53513b9d
--- /dev/null
+++ b/research/bit-three-register-lower-bound/catalogue.py
@@ -0,0 +1,185 @@
+"""Enumerate the complete scalar catalogue with explicit binary subspaces."""
+
+from collections import Counter
+from itertools import combinations, combinations_with_replacement, permutations, product
+from math import factorial, prod
+from pathlib import Path
+import argparse
+import json
+
+VECTORS = tuple(product(range(2), repeat=3))
+POSITIONS = ((0, 1), (0, 2), (1, 0), (1, 2), (2, 0), (2, 1))
+
+
+def xor(a, b):
+    return tuple(x ^ y for x, y in zip(a, b))
+
+
+def closure(vectors, zero):
+    space = {zero}
+    for vector in vectors:
+        space.update([xor(vector, other) for other in space])
+    return frozenset(space)
+
+
+def multiply(a, b):
+    return tuple(sum(a[3 * i + k] * b[3 * k + j] for k in range(3)) % 2
+                 for i in range(3) for j in range(3))
+
+
+def action(matrix, vector):
+    return tuple(sum(matrix[3 * i + j] * vector[j] for j in range(3)) % 2
+                 for i in range(3))
+
+
+def encode(matrix):
+    return sum(value << i for i, value in enumerate(matrix))
+
+
+def off(matrix):
+    return tuple(matrix[3 * i + j] for i, j in POSITIONS)
+
+
+def transpose(matrix):
+    return tuple(matrix[3 * j + i] for i in range(3) for j in range(3))
+
+
+def conjugate(matrix, permutation):
+    return tuple(matrix[3 * permutation[i] + permutation[j]]
+                 for i in range(3) for j in range(3))
+
+
+def catalogue():
+    projectors = tuple(sorted(
+        (matrix for matrix in product(range(2), repeat=9)
+         if len({action(matrix, v) for v in VECTORS}) == 2
+         and multiply(matrix, matrix) == matrix), key=encode))
+    by_code = {encode(matrix): matrix for matrix in projectors}
+    patterns = set(map(off, projectors)) - {(0,) * 6}
+    lines = [frozenset(((0, 0, 0), v)) for v in VECTORS[1:]]
+    planes = {closure((v, w), (0, 0, 0)) for v, w in combinations(VECTORS[1:], 2)}
+    subspaces = set(lines) | planes
+    flags = [(line, plane) for line in lines for plane in planes if line <= plane]
+    invariant = {matrix: {space for space in subspaces
+                         if all(action(matrix, v) in space for v in space)}
+                 for matrix in projectors}
+
+    def preserved(pool):
+        return set.intersection(subspaces, *(invariant[matrix] for matrix in pool))
+
+    def has_flag(pool):
+        common = preserved(pool)
+        return any(line in common and plane in common for line, plane in flags)
+
+    def transforms(pool):
+        return [tuple(encode(conjugate(transpose(matrix) if turn else matrix, permutation))
+                      for matrix in pool)
+                for permutation in permutations(range(3)) for turn in (False, True)]
+
+    spaces = {frozenset(((0,) * 6,))}
+    for _ in range(3):
+        spaces |= {closure(tuple(space) + (pattern,), (0,) * 6)
+                   for space in spaces for pattern in patterns}
+    records = []
+    for space in spaces:
+        pool = tuple(matrix for matrix in projectors if off(matrix) in space)
+        noncoordinate = tuple(matrix for matrix in pool if any(off(matrix)))
+        records.append((len(space).bit_length() - 1, pool, noncoordinate))
+
+    # This complete classification is the finite step in the span-six proof.
+    exceptions = []
+    for dimension, _, pool in records:
+        if dimension == 2 and not has_flag(pool):
+            assert len(pool) == 3
+            assert all(has_flag(pool[:i] + pool[i + 1:]) for i in range(3))
+            exceptions.append(tuple(sorted(map(encode, pool))))
+    representative = tuple(by_code[code] for code in (3, 18, 360))
+    expected_exceptions = {tuple(sorted(encode(conjugate(matrix, permutation))
+                                        for matrix in representative))
+                           for permutation in permutations(range(3))}
+    assert set(exceptions) == expected_exceptions and len(exceptions) == 6
+    basis = (1, 1, 0, 0, 1, 0, 0, 0, 1)
+    assert [encode(multiply(multiply(basis, matrix), basis))
+            for matrix in representative] == [1, 16, 511]
+
+    canonical = {min(tuple(sorted(row)) for row in transforms(pool)): dimension
+                 for dimension, pool, _ in records}
+    counts = Counter()
+    passing = []
+    dimensions = Counter()
+    for key, dimension in sorted(canonical.items(), key=lambda item: (item[1], item[0])):
+        if dimension == 0:
+            continue
+        pool = tuple(by_code[code] for code in key)
+        stabilizers = [tuple(key.index(code) for code in row)
+                       for row in transforms(pool) if set(row) == set(key)]
+        support_cache = {}
+        # Precompute the seven row and seven column incidence tests.
+        incidence = []
+        for vector in VECTORS[1:]:
+            for turn in (False, True):
+                images = []
+                for matrix in pool:
+                    matrix = transpose(matrix) if turn else matrix
+                    image = next(action(matrix, v) for v in VECTORS if any(action(matrix, v)))
+                    images.append(sum(a * b for a, b in zip(vector, image)) % 2)
+                incidence.append(images)
+        for gates in combinations_with_replacement(range(len(pool)), 8):
+            support = tuple(sorted(set(gates)))
+            if support not in support_cache:
+                occurring = tuple(pool[i] for i in support)
+                full = len(closure(map(off, occurring), (0,) * 6)) == 2 ** dimension
+                irreducible = not preserved(occurring)
+                flag = dimension < 3 or not has_flag(tuple(m for m in occurring if any(off(m))))
+                support_cache[support] = full, irreducible, flag
+            full, irreducible, flag = support_cache[support]
+            if not full:
+                continue
+            counts['full_rank_multiplicities'] += 1
+            if not irreducible:
+                continue
+            counts['irreducible'] += 1
+            if not flag:
+                continue
+            counts['nonzero_flag_pass'] += 1
+            multiplicity = tuple(gates.count(i) for i in range(len(pool)))
+            if any(sum(a * b for a, b in zip(multiplicity, row)) < 3 for row in incidence):
+                continue
+            counts['distance_pass'] += 1
+            variants = []
+            for symmetry in stabilizers:
+                transformed = [0] * len(pool)
+                for i, j in enumerate(symmetry):
+                    transformed[j] = multiplicity[i]
+                variants.append(tuple(transformed))
+            if multiplicity != min(variants):
+                continue
+            counts['canonical_pass'] += 1
+            dimensions[dimension] += 1
+            passing.append(tuple(key[i] for i in gates))
+    passing.sort()
+    assert len(passing) == len(set(passing))
+    summary = {
+        'projectors': len(projectors), 'nonzero_patterns': len(patterns),
+        'pattern_spaces': dict(sorted(Counter(dimension for dimension, _, _ in records).items())),
+        'complete_flags': len(flags), 'plane_exceptions': len(exceptions),
+        'canonical_pools': len(canonical), 'filter_counts': dict(counts),
+        'families_by_dimension': dict(sorted(dimensions.items())), 'families': len(passing),
+        'orders': sum(factorial(8) // prod(factorial(c) for c in Counter(word).values())
+                      for word in passing),
+    }
+    return passing, summary
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('output', type=Path)
+    args = parser.parse_args()
+    words, summary = catalogue()
+    args.output.write_text(str(len(words)) + '\n' + ''.join(
+        ' '.join(map(str, (index, *word))) + '\n' for index, word in enumerate(words)))
+    print(json.dumps(summary, indent=2))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/research/bit-three-register-lower-bound/controls.py b/research/bit-three-register-lower-bound/controls.py
new file mode 100644
index 0000000000..6165600831
--- /dev/null
+++ b/research/bit-three-register-lower-bound/controls.py
@@ -0,0 +1,89 @@
+"""Check formal identities and the full-rank target boundary independently."""
+
+from fractions import Fraction
+from itertools import combinations
+
+from catalogue import multiply
+
+IDENTITY = (1, 0, 0, 0, 1, 0, 0, 0, 1)
+ZERO = (0,) * 9
+
+
+def unit(row, column):
+    return tuple(int((i, j) == (row, column)) for i in range(3) for j in range(3))
+
+
+def rank(rows):
+    rows = [list(map(Fraction, row)) for row in rows]
+    result = 0
+    for column in range(len(rows[0])):
+        pivot = next((i for i in range(result, len(rows)) if rows[i][column]), None)
+        if pivot is None:
+            continue
+        rows[result], rows[pivot] = rows[pivot], rows[result]
+        scale = rows[result][column]
+        rows[result] = [x / scale for x in rows[result]]
+        for i in range(result + 1, len(rows)):
+            scale = rows[i][column]
+            rows[i] = [x - scale * y for x, y in zip(rows[i], rows[result])]
+        result += 1
+        if result == len(rows):
+            break
+    return result
+
+
+def matrix_rank(matrix):
+    return rank([matrix[3 * row:3 * row + 3] for row in range(3)])
+
+
+def formal_product(operations):
+    terms = {ZERO: IDENTITY}
+    for frame, projector in operations:
+        updated = {}
+        complement = tuple(x ^ y for x, y in zip(IDENTITY, projector))
+        for address, coefficient in terms.items():
+            for shift, scalar in ((ZERO, complement), (frame, projector)):
+                exponent = tuple(x + y for x, y in zip(address, shift))
+                matrix = multiply(coefficient, scalar)
+                previous = updated.get(exponent, ZERO)
+                updated[exponent] = tuple(x ^ y for x, y in zip(previous, matrix))
+        terms = {address: matrix for address, matrix in updated.items() if any(matrix)}
+    return terms
+
+
+def scalar(code):
+    return tuple((code >> i) & 1 for i in range(9))
+
+
+def check_formal_controls():
+    baseline = [(unit(i, i), unit(w, w)) for w in range(3) for i in range(3)]
+    assert formal_product(baseline) == {IDENTITY: IDENTITY}
+    # This rank-one address change is hidden by an invalid reduction modulo two.
+    alias = [(tuple(3 * x for x in baseline[0][0]), baseline[0][1]), *baseline[1:]]
+    assert formal_product(alias) != {IDENTITY: IDENTITY}
+    assert formal_product([(tuple(x % 2 for x in frame), projector)
+                           for frame, projector in alias]) == {IDENTITY: IDENTITY}
+
+    e00, e11, e22, e10, e01 = (unit(i, j) for i, j in ((0, 0), (1, 1), (2, 2), (1, 0), (0, 1)))
+    plus = lambda a, b: tuple(x + y for x, y in zip(a, b))
+    minus = lambda a, b: tuple(x - y for x, y in zip(a, b))
+    frames = [e22, e00, e11, plus(e00, e10), minus(e11, e10),
+              plus(e00, e01), minus(e11, e01), tuple(-x for x in e22)]
+    projectors = list(map(scalar, (3, 1, 1, 16, 16, 256, 256, 3)))
+    assert all(matrix_rank(frame) == 1 for frame in frames)
+    assert all(matrix_rank(projector) == 1 and multiply(projector, projector) == projector
+               for projector in projectors)
+    assert any(multiply(a, b) != multiply(b, a) for a, b in combinations(projectors, 2))
+    assert rank(frames) == 5
+    target = plus(e00, e11)
+    assert matrix_rank(target) == 2
+    assert formal_product(list(zip(frames, projectors))) == {target: IDENTITY}
+    return {'nine_call_identity': True, 'modular_alias_rejected': True,
+            'noncommuting_eight_call_identity': True, 'frame_span': 5, 'target_rank': 2}
+
+
+# These are search-core controls, not admissible constructions. Code 17 is
+# a rank-two scalar projector. A real accepted core example detects a checker
+# that rejects every terminal; the mutation removes one target contribution.
+CORE_WORDS = ((3, 3, 17, 17, 17, 256, 256, 256),
+              (3, 3, 17, 17, 1, 256, 256, 256))
diff --git a/research/bit-three-register-lower-bound/expected.json b/research/bit-three-register-lower-bound/expected.json
new file mode 100644
index 0000000000..d2c4f09cb5
--- /dev/null
+++ b/research/bit-three-register-lower-bound/expected.json
@@ -0,0 +1,23 @@
+{
+  "projectors": 28,
+  "nonzero_patterns": 19,
+  "pattern_spaces": {
+    "0": 1,
+    "1": 19,
+    "2": 129,
+    "3": 363
+  },
+  "complete_flags": 21,
+  "plane_exceptions": 6,
+  "canonical_pools": 67,
+  "families_by_dimension": {
+    "1": 4,
+    "2": 137,
+    "3": 2135
+  },
+  "families": 2276,
+  "orders": 21317520,
+  "diagonal_terminals": 415398,
+  "survivors": 0,
+  "catalogue_sha256": "8be9e2582199e127b061942cf0fba9156b7560312f20a87dc184fe0010d7cead"
+}
diff --git a/research/bit-three-register-lower-bound/f101.cpp b/research/bit-three-register-lower-bound/f101.cpp
new file mode 100644
index 0000000000..e3fad6c59c
--- /dev/null
+++ b/research/bit-three-register-lower-bound/f101.cpp
@@ -0,0 +1,449 @@
+// Check subset collisions with exact row reduction over F101.
+#include <algorithm>
+#include <array>
+#include <chrono>
+#include <cstdint>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <numeric>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <unordered_set>
+#include <utility>
+#include <vector>
+
+using Row = std::array<unsigned char, 8>;
+
+int prod[101][101], inverse[101], num[101], den[101], mm[512][512];
+
+int minus_product(int a, int t, int b) {
+    int q = a - prod[t][b];
+    return q < 0 ? q + 101 : q;
+}
+
+struct Basis {
+    std::array<Row, 3> a{};
+    std::array<unsigned char, 3> p{};
+    int n = 0;
+
+    bool operator==(const Basis& b) const { return a == b.a; }
+};
+
+struct BasisHash {
+    size_t operator()(const Basis& b) const {
+        size_t h = 0xcbf29ce484222325ULL;
+        for (const auto& r : b.a)
+            for (int x : r)
+                h = (h ^ x) * 0x100000001b3ULL;
+        return h;
+    }
+};
+
+Row subset(int m) {
+    Row v{};
+    for (int i = 0; i < 8; ++i)
+        v[i] = (m >> i) & 1;
+    return v;
+}
+
+Row residual(Row v, const Basis& b) {
+    for (int i = 0; i < b.n; ++i) {
+        int t = v[b.p[i]];
+        if (t)
+            for (int j = b.p[i]; j < 8; ++j)
+                v[j] = minus_product(v[j], t, b.a[i][j]);
+    }
+    return v;
+}
+
+Basis insert(Basis b, Row v) {
+    v = residual(v, b);
+    int p = 0;
+    while (p < 8 && !v[p])
+        ++p;
+    if (p == 8)
+        return b;
+    if (b.n == 3)
+        throw std::logic_error("A state has too many relations");
+
+    int q = inverse[v[p]];
+    for (int j = p; j < 8; ++j)
+        v[j] = prod[v[j]][q];
+    for (int i = 0; i < b.n; ++i) {
+        int t = b.a[i][p];
+        if (t)
+            for (int j = p; j < 8; ++j)
+                b.a[i][j] = minus_product(b.a[i][j], t, v[j]);
+    }
+    int at = b.n;
+    while (at && b.p[at - 1] > p) {
+        b.a[at] = b.a[at - 1];
+        b.p[at] = b.p[at - 1];
+        --at;
+    }
+    b.a[at] = v;
+    b.p[at] = p;
+    ++b.n;
+    return b;
+}
+
+int binrank(const std::vector<int>& rows) {
+    int b[8]{}, n = 0;
+    for (int v : rows) {
+        for (int k = 7; k >= 0; --k) {
+            if (!((v >> k) & 1))
+                continue;
+            if (b[k]) {
+                v ^= b[k];
+            } else {
+                b[k] = v;
+                ++n;
+                break;
+            }
+        }
+    }
+    return n;
+}
+
+int primitive_parity(const Row& r) {
+    int d = 1;
+    for (int x : r) {
+        if (!den[x])
+            throw std::logic_error("Rational reconstruction exceeds the proved bound");
+        d = std::lcm(d, den[x]);
+    }
+    int v[8], g = 0;
+    for (int i = 0; i < 8; ++i) {
+        v[i] = num[r[i]] * (d / den[r[i]]);
+        g = std::gcd(g, std::abs(v[i]));
+    }
+    int bits = 0;
+    for (int i = 0; i < 8; ++i)
+        bits |= (std::abs(v[i] / g) & 1) << i;
+    return bits;
+}
+
+struct GeneralBasis {
+    std::array<Row, 6> a{};
+    std::array<int, 6> p{};
+    int n = 0;
+
+    explicit GeneralBasis(const Basis& b) {
+        n = b.n;
+        for (int i = 0; i < n; ++i) {
+            a[i] = b.a[i];
+            p[i] = b.p[i];
+        }
+    }
+
+    bool add(Row v) {
+        for (int i = 0; i < n; ++i) {
+            int t = v[p[i]];
+            if (t)
+                for (int j = p[i]; j < 8; ++j)
+                    v[j] = minus_product(v[j], t, a[i][j]);
+        }
+        int j = 0;
+        while (j < 8 && !v[j])
+            ++j;
+        if (j == 8)
+            return false;
+        int q = inverse[v[j]];
+        for (int k = j; k < 8; ++k)
+            v[k] = prod[v[k]][q];
+        int at = n;
+        while (at && p[at - 1] > j) {
+            a[at] = a[at - 1];
+            p[at] = p[at - 1];
+            --at;
+        }
+        a[at] = v;
+        p[at] = j;
+        ++n;
+        return true;
+    }
+};
+
+bool target_small(int mask, const Basis& b) {
+    if (__builtin_popcount(static_cast<unsigned>(mask)) <= 2)
+        return true;
+    for (int i = 0; i < 8; ++i) {
+        for (int j = i + 1; j < 8; ++j) {
+            GeneralBasis g(b);
+            g.add(subset(1 << i));
+            g.add(subset(1 << j));
+            if (!g.add(subset(mask)))
+                return true;
+        }
+    }
+    return false;
+}
+
+struct Group {
+    uint64_t key;
+    int mask, coefficient;
+};
+
+struct Search {
+    std::array<int, 8> word;
+    std::vector<std::pair<int, int>> terms;
+    std::vector<int> off;
+    std::unordered_set<Basis, BasisHash> visited;
+    uint64_t states[4]{}, parity = 0, zero = 0, terminals = 0;
+    uint64_t robust = 0, not_monomial = 0;
+
+    void start() {
+        int coefficients[256]{};
+        coefficients[0] = 273;
+        for (int j = 0; j < 8; ++j) {
+            for (int m = 0; m < (1 << j); ++m) {
+                int c = coefficients[m];
+                coefficients[m] = mm[c][273 ^ word[j]];
+                coefficients[m | (1 << j)] = mm[c][word[j]];
+            }
+        }
+        for (int m = 0; m < 256; ++m)
+            if (coefficients[m])
+                terms.emplace_back(m, coefficients[m]);
+        for (int bit : {1, 2, 3, 5, 6, 7}) {
+            int v = 0;
+            for (int j = 0; j < 8; ++j)
+                v |= ((word[j] >> bit) & 1) << j;
+            off.push_back(v);
+        }
+        go(Basis{});
+    }
+
+    void go(const Basis& b) {
+        if (!visited.insert(b).second)
+            return;
+        ++states[b.n];
+        auto ps = off;
+        for (int i = 0; i < b.n; ++i)
+            ps.push_back(primitive_parity(b.a[i]));
+        if (binrank(ps) > 3) {
+            ++parity;
+            return;
+        }
+        for (int i = 0; i < b.n; ++i) {
+            int entries = 0;
+            for (int v : b.a[i])
+                entries += v != 0;
+            if (entries == 1) {
+                ++zero;
+                return;
+            }
+        }
+
+        std::unordered_map<uint64_t, int> positions;
+        std::vector<Group> groups;
+        positions.reserve(terms.size());
+        groups.reserve(terms.size());
+        for (auto [mask, c] : terms) {
+            auto v = residual(subset(mask), b);
+            uint64_t key = 0;
+            for (int j = 0; j < 8; ++j)
+                key |= uint64_t(v[j]) << (7 * j);
+            auto [it, fresh] = positions.emplace(key, int(groups.size()));
+            if (fresh)
+                groups.push_back({key, mask, c});
+            else
+                groups[it->second].coefficient ^= c;
+        }
+
+        std::vector<Group> active;
+        int counts[9]{};
+        for (auto g : groups) {
+            if (g.coefficient) {
+                active.push_back(g);
+                for (int bit : {1, 2, 3, 5, 6, 7})
+                    counts[bit] += (g.coefficient >> bit) & 1;
+            }
+        }
+        int chosen = -1;
+        for (int bit : {1, 2, 3, 5, 6, 7})
+            if (counts[bit] && (chosen < 0 || counts[bit] < counts[chosen]))
+                chosen = bit;
+
+        if (chosen < 0) {
+            int target[3] = {-1, -1, -1};
+            for (auto g : active) {
+                for (int w = 0; w < 3; ++w) {
+                    if ((g.coefficient >> (4 * w)) & 1) {
+                        if (target[w] >= 0) {
+                            ++not_monomial;
+                            return;
+                        }
+                        target[w] = g.mask;
+                    }
+                }
+            }
+            for (int m : target) {
+                if (m < 0) {
+                    ++not_monomial;
+                    return;
+                }
+            }
+            ++terminals;
+            for (int m : target)
+                if (target_small(m, b))
+                    return;
+            ++robust;
+            std::cerr << "ROBUST";
+            for (int n : word)
+                std::cerr << ' ' << n;
+            std::cerr << '\n';
+            return;
+        }
+
+        if (b.n == 3)
+            return;
+        int first = 0;
+        while (!((active[first].coefficient >> chosen) & 1))
+            ++first;
+        for (int j = first + 1; j < int(active.size()); ++j) {
+            if ((active[j].coefficient >> chosen) & 1) {
+                Row v{};
+                for (int k = 0; k < 8; ++k) {
+                    int t = ((active[first].mask >> k) & 1) -
+                            ((active[j].mask >> k) & 1);
+                    v[k] = t < 0 ? 100 : t;
+                }
+                auto c = insert(b, v);
+                if (c.n != b.n + 1)
+                    throw std::logic_error("A partner is already in the same group");
+                go(c);
+            }
+        }
+    }
+};
+
+void setup() {
+    for (int a = 0; a < 101; ++a)
+        for (int b = 0; b < 101; ++b)
+            prod[a][b] = (a * b) % 101;
+    for (int a = 1; a < 101; ++a)
+        for (int b = 1; b < 101; ++b)
+            if (prod[a][b] == 1)
+                inverse[a] = b;
+    for (int d = 1; d <= 5; ++d) {
+        for (int n = -5; n <= 5; ++n) {
+            if (std::gcd(std::abs(n), d) == 1) {
+                int x = prod[(n + 101) % 101][inverse[d]];
+                if (den[x] && (num[x] != n || den[x] != d))
+                    throw std::logic_error("Rational reconstruction is not unique");
+                num[x] = n;
+                den[x] = d;
+            }
+        }
+    }
+    for (int a = 0; a < 512; ++a) {
+        for (int b = 0; b < 512; ++b) {
+            int c = 0;
+            for (int i = 0; i < 3; ++i) {
+                for (int j = 0; j < 3; ++j) {
+                    int t = 0;
+                    for (int k = 0; k < 3; ++k)
+                        t ^= ((a >> (3 * i + k)) & 1) & ((b >> (3 * k + j)) & 1);
+                    c |= t << (3 * i + j);
+                }
+            }
+            mm[a][b] = c;
+        }
+    }
+}
+
+struct InputCase {
+    int id;
+    std::array<int, 8> word;
+};
+
+int main(int argc, char** argv) {
+    try {
+        if (argc != 3)
+            throw std::runtime_error("Usage: f101 INPUT RESULTS");
+        setup();
+        std::ifstream input(argv[1]);
+        if (!input)
+            throw std::runtime_error("Cannot open input: " + std::string(argv[1]));
+        int count;
+        if (!(input >> count) || count <= 0)
+            throw std::runtime_error("The case count must be a positive integer");
+
+        std::vector<InputCase> cases;
+        std::unordered_set<int> ids;
+        for (int c = 0; c < count; ++c) {
+            InputCase entry;
+            if (!(input >> entry.id))
+                throw std::runtime_error("Incomplete input: missing case ID");
+            if (!ids.insert(entry.id).second)
+                throw std::runtime_error("Duplicate case ID: " + std::to_string(entry.id));
+            for (int& n : entry.word) {
+                if (!(input >> n))
+                    throw std::runtime_error("Incomplete input: missing projector");
+                // Controls can use rank-two idempotents. Their acceptance is
+                // not a physical construction. The catalogue checks rank one.
+                if (n <= 0 || n >= 512 || mm[n][n] != n)
+                    throw std::runtime_error("Each projector must be a nonzero 3-by-3 binary idempotent");
+            }
+            cases.push_back(entry);
+        }
+        input >> std::ws;
+        auto trailing = input.peek();
+        if (input.bad())
+            throw std::runtime_error("Cannot read input: " + std::string(argv[1]));
+        if (trailing != std::ifstream::traits_type::eof())
+            throw std::runtime_error("Input has extra data after the last case");
+
+        std::ofstream output(argv[2]);
+        if (!output)
+            throw std::runtime_error("Cannot open output: " + std::string(argv[2]));
+        uint64_t allwords = 0, allrob = 0;
+        auto begun = std::chrono::steady_clock::now();
+        for (int c = 0; c < count; ++c) {
+            auto ns = cases[c].word;
+            std::sort(ns.begin(), ns.end());
+            uint64_t totals[4]{}, terminal = 0, rob = 0, parity = 0;
+            uint64_t zero = 0, words = 0, invalid = 0;
+            do {
+                Search s;
+                s.word = ns;
+                s.start();
+                for (int k = 0; k < 4; ++k)
+                    totals[k] += s.states[k];
+                terminal += s.terminals;
+                rob += s.robust;
+                parity += s.parity;
+                zero += s.zero;
+                invalid += s.not_monomial;
+                ++words;
+            } while (std::next_permutation(ns.begin(), ns.end()));
+            output << "{\"case\":" << cases[c].id << ",\"words\":" << words
+                   << ",\"states\":[" << totals[0] << ',' << totals[1] << ','
+                   << totals[2] << ',' << totals[3] << "],\"terminals\":" << terminal
+                   << ",\"robust\":" << rob << ",\"parity_prunes\":" << parity
+                   << ",\"zero_frame_prunes\":" << zero
+                   << ",\"nonmonomial_diagonal\":" << invalid << "}\n";
+            output.flush();
+            if (!output)
+                throw std::runtime_error("Cannot write output: " + std::string(argv[2]));
+            allwords += words;
+            allrob += rob;
+            if ((c + 1) % 10 == 0 || c + 1 == count) {
+                double seconds = std::chrono::duration<double>(
+                    std::chrono::steady_clock::now() - begun).count();
+                std::cerr << c + 1 << " / " << count << " families, " << allwords
+                          << " words, " << allrob << " robust, " << seconds << " seconds\n";
+            }
+        }
+        output.close();
+        if (!output)
+            throw std::runtime_error("Cannot close output: " + std::string(argv[2]));
+        return 0;
+    } catch (const std::exception& error) {
+        std::cerr << error.what() << '\n';
+        return 1;
+    }
+}
diff --git a/research/bit-three-register-lower-bound/rational.cpp b/research/bit-three-register-lower-bound/rational.cpp
new file mode 100644
index 0000000000..09bdb3a90f
--- /dev/null
+++ b/research/bit-three-register-lower-bound/rational.cpp
@@ -0,0 +1,403 @@
+// Enumerate the required rational subset-collision spaces of rank at most three.
+#include <algorithm>
+#include <array>
+#include <cstdint>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <numeric>
+#include <set>
+#include <stdexcept>
+#include <string>
+#include <unordered_map>
+#include <unordered_set>
+#include <vector>
+
+using V = std::array<int, 8>;
+
+// Use a primitive integer row with a positive leading entry.
+void normalize(V& v) {
+    int divisor = 0;
+    for (int x : v) divisor = std::gcd(divisor, std::abs(x));
+    if (!divisor) return;
+    int leading = 0;
+    while (!v[leading]) ++leading;
+    if (v[leading] < 0) divisor = -divisor;
+    for (int& x : v) x /= divisor;
+}
+
+int pivot(const V& v) {
+    for (int i = 0; i < 8; ++i) {
+        if (v[i]) return i;
+    }
+    return 8;
+}
+
+struct Rows {
+    std::array<V, 3> rows{};
+    int size = 0;
+    bool operator==(const Rows& other) const { return rows == other.rows; }
+};
+
+struct HashRows {
+    size_t operator()(const Rows& basis) const {
+        size_t hash = 0;
+        for (const auto& row : basis.rows) {
+            for (int x : row) hash = (hash ^ size_t(x + 97)) * 1099511628211ull;
+        }
+        return hash;
+    }
+};
+
+V reduce(V v, const Rows& basis) {
+    for (int a = 0; a < basis.size; ++a) {
+        const V& row = basis.rows[a];
+        int p = pivot(row), coefficient = v[p];
+        if (coefficient) {
+            for (int i = 0; i < 8; ++i) {
+                v[i] = row[p] * v[i] - coefficient * row[i];
+            }
+            normalize(v);
+        }
+    }
+    return v;
+}
+
+Rows append(Rows basis, V v) {
+    v = reduce(v, basis);
+    if (pivot(v) == 8) return basis;
+    normalize(v);
+    int p = pivot(v);
+    for (int a = 0; a < basis.size; ++a) {
+        int coefficient = basis.rows[a][p];
+        if (coefficient) {
+            for (int i = 0; i < 8; ++i) {
+                basis.rows[a][i] = v[p] * basis.rows[a][i] - coefficient * v[i];
+            }
+            normalize(basis.rows[a]);
+        }
+    }
+    basis.rows[basis.size++] = v;
+    std::sort(basis.rows.begin(), basis.rows.begin() + basis.size,
+              [](const V& a, const V& b) { return pivot(a) < pivot(b); });
+    return basis;
+}
+
+V vec(int mask) {
+    V v{};
+    for (int i = 0; i < 8; ++i) v[i] = (mask >> i) & 1;
+    return v;
+}
+
+int binaryrank(std::vector<int> rows) {
+    int basis[8]{}, rank = 0;
+    for (int v : rows) {
+        for (int i = 7; i >= 0; --i) {
+            if (!(v >> i & 1)) continue;
+            if (basis[i]) {
+                v ^= basis[i];
+            } else {
+                basis[i] = v;
+                ++rank;
+                break;
+            }
+        }
+    }
+    return rank;
+}
+
+// Binary 3 by 3 matrices use row-major bit order. The identity is 273.
+int mul(int a, int b) {
+    int result = 0;
+    for (int i = 0; i < 3; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            int entry = 0;
+            for (int k = 0; k < 3; ++k) {
+                entry ^= ((a >> (3 * i + k)) & 1) * ((b >> (3 * k + j)) & 1);
+            }
+            result |= entry << (3 * i + j);
+        }
+    }
+    return result;
+}
+
+struct KeyHash {
+    size_t operator()(const V& v) const {
+        size_t hash = 0;
+        for (int x : v) hash = (hash ^ size_t(x + 4097)) * 1099511628211ull;
+        return hash;
+    }
+};
+
+struct Group {
+    V key;
+    int mask, matrix;
+};
+
+struct Solver {
+    std::array<int, 8> ns;
+    std::vector<std::pair<int, int>> terms;
+    std::vector<int> offrows;
+    std::unordered_set<Rows, HashRows> seen;
+    uint64_t states[4]{}, parityprune = 0, zeroprune = 0, terminal = 0, robust = 0;
+    std::ofstream* cert = nullptr;
+    int wordidx = 0;
+
+    void init() {
+        terms = {{0, 273}};
+        for (int j = 0; j < 8; ++j) {
+            std::vector<std::pair<int, int>> next;
+            for (auto [mask, matrix] : terms) {
+                int constant = mul(matrix, 273 ^ ns[j]);
+                int variable = mul(matrix, ns[j]);
+                if (constant) next.emplace_back(mask, constant);
+                if (variable) next.emplace_back(mask | (1 << j), variable);
+            }
+            terms = next;
+        }
+        for (int bit : {1, 2, 3, 5, 6, 7}) {
+            int row = 0;
+            for (int i = 0; i < 8; ++i) row |= ((ns[i] >> bit) & 1) << i;
+            offrows.push_back(row);
+        }
+    }
+
+    // Check whether a target is a linear combination of at most two physical frames.
+    bool target_rank_small(int mask, const Rows& relations) {
+        for (int i = 0; i < 8; ++i) {
+            for (int j = i + 1; j < 8; ++j) {
+                std::vector<V> basis(relations.rows.begin(),
+                                     relations.rows.begin() + relations.size);
+                auto add = [&](V v) {
+                    for (const auto& row : basis) {
+                        int p = pivot(row), coefficient = v[p];
+                        if (coefficient) {
+                            for (int k = 0; k < 8; ++k) {
+                                v[k] = row[p] * v[k] - coefficient * row[k];
+                            }
+                            normalize(v);
+                        }
+                    }
+                    if (pivot(v) < 8) {
+                        normalize(v);
+                        basis.push_back(v);
+                        std::sort(basis.begin(), basis.end(),
+                                  [](const V& a, const V& b) { return pivot(a) < pivot(b); });
+                    }
+                };
+                add(vec(1 << i));
+                add(vec(1 << j));
+                int before = basis.size();
+                add(vec(mask));
+                if (int(basis.size()) == before) return true;
+            }
+        }
+        return false;
+    }
+
+    void go(const Rows& relations) {
+        if (!seen.insert(relations).second) return;
+        ++states[relations.size];
+
+        auto parity = offrows;
+        for (int a = 0; a < relations.size; ++a) {
+            int row = 0;
+            for (int j = 0; j < 8; ++j) {
+                row |= (std::abs(relations.rows[a][j]) & 1) << j;
+            }
+            parity.push_back(row);
+        }
+        if (binaryrank(parity) > 3) {
+            ++parityprune;
+            return;
+        }
+        for (int i = 0; i < 8; ++i) {
+            if (pivot(reduce(vec(1 << i), relations)) == 8) {
+                ++zeroprune;
+                return;
+            }
+        }
+
+        // Form exact quotient keys with one common integer denominator.
+        int denominator = 1;
+        for (int a = 0; a < relations.size; ++a) {
+            denominator = std::lcm(denominator, relations.rows[a][pivot(relations.rows[a])]);
+        }
+        std::unordered_map<V, std::pair<int, int>, KeyHash> groups;
+        for (auto [mask, matrix] : terms) {
+            V key = vec(mask);
+            for (int& x : key) x *= denominator;
+            for (int a = 0; a < relations.size; ++a) {
+                int p = pivot(relations.rows[a]);
+                if (mask >> p & 1) {
+                    for (int j = 0; j < 8; ++j) {
+                        key[j] -= (denominator / relations.rows[a][p]) * relations.rows[a][j];
+                    }
+                }
+            }
+            auto [it, inserted] = groups.emplace(key, std::make_pair(mask, 0));
+            it->second.second ^= matrix;
+        }
+        std::vector<Group> active;
+        bool has_offdiagonal = false;
+        for (const auto& [key, term] : groups) {
+            if (term.second) {
+                active.push_back({key, term.first, term.second});
+                has_offdiagonal |= (term.second & ~273) != 0;
+            }
+        }
+
+        if (!has_offdiagonal) {
+            int targets[3] = {-1, -1, -1};
+            for (const auto& group : active) {
+                for (int w = 0; w < 3; ++w) {
+                    if (group.matrix >> (4 * w) & 1) {
+                        if (targets[w] != -1) return;
+                        targets[w] = group.mask;
+                    }
+                }
+            }
+            for (int target : targets) {
+                if (target < 0) return;
+            }
+            ++terminal;
+            bool survives = true;
+            for (int target : targets) survives &= !target_rank_small(target, relations);
+            if (survives) {
+                ++robust;
+                if (cert) {
+                    *cert << wordidx;
+                    for (int projector : ns) *cert << ' ' << projector;
+                    *cert << " : ";
+                    for (int a = 0; a < relations.size; ++a) {
+                        for (int x : relations.rows[a]) *cert << x << ',';
+                        *cert << " ; ";
+                    }
+                    *cert << " : ";
+                    for (int target : targets) *cert << target << ' ';
+                    *cert << '\n';
+                }
+            }
+            return;
+        }
+        if (relations.size == 3) return;
+
+        // A nonzero off-diagonal entry must cancel with a partner group.
+        int best = -1, best_bit = -1, best_count = 1000;
+        for (int a = 0; a < int(active.size()); ++a) {
+            for (int bit : {1, 2, 3, 5, 6, 7}) {
+                if (!(active[a].matrix >> bit & 1)) continue;
+                int count = 0;
+                for (int b = 0; b < int(active.size()); ++b) {
+                    if (a != b && (active[b].matrix >> bit & 1)) ++count;
+                }
+                if (count < best_count) {
+                    best = a;
+                    best_bit = bit;
+                    best_count = count;
+                }
+            }
+        }
+        for (int b = 0; b < int(active.size()); ++b) {
+            if (b != best && (active[b].matrix >> best_bit & 1)) {
+                V difference{};
+                for (int i = 0; i < 8; ++i) {
+                    difference[i] = ((active[best].mask >> i) & 1) - ((active[b].mask >> i) & 1);
+                }
+                Rows next = append(relations, difference);
+                if (next.size == relations.size) {
+                    throw std::logic_error("Partner already in this group");
+                }
+                go(next);
+            }
+        }
+    }
+};
+
+int main(int argc, char** argv) {
+    try {
+        if (argc != 5 || std::string(argv[1]) != "--batch") {
+            throw std::runtime_error("Usage: rational --batch INPUT RESULTS CANDIDATES");
+        }
+        std::ifstream input(argv[2]);
+        if (!input) throw std::runtime_error("Cannot open batch input");
+        int case_count;
+        if (!(input >> case_count) || case_count <= 0) {
+            throw std::runtime_error("Invalid batch case count");
+        }
+        std::ofstream results(argv[3]), candidates(argv[4]);
+        if (!results || !candidates) throw std::runtime_error("Cannot open result paths");
+
+        std::set<int> case_ids;
+        uint64_t total_words = 0, total_robust = 0;
+        for (int index = 0; index < case_count; ++index) {
+            int case_id;
+            std::array<int, 8> projectors;
+            if (!(input >> case_id) || !case_ids.insert(case_id).second) {
+                throw std::runtime_error("Invalid or duplicate case id");
+            }
+            for (int& projector : projectors) {
+                if (!(input >> projector)) throw std::runtime_error("Incomplete projector row");
+            }
+            for (int projector : projectors) {
+                if (projector <= 0 || projector >= 512 || mul(projector, projector) != projector) {
+                    throw std::runtime_error("Invalid projector");
+                }
+                int nonzero_row = 0;
+                for (int i = 0; i < 3; ++i) {
+                    int row = (projector >> (3 * i)) & 7;
+                    if (row) {
+                        if (nonzero_row && nonzero_row != row) {
+                            throw std::runtime_error("Projector rank is not one");
+                        }
+                        nonzero_row = row;
+                    }
+                }
+            }
+
+            std::sort(projectors.begin(), projectors.end());
+            uint64_t states[4]{}, terminals = 0, robust = 0, parity = 0, zero = 0;
+            int words = 0;
+            do {
+                Solver solver;
+                solver.ns = projectors;
+                solver.wordidx = words;
+                solver.cert = &candidates;
+                solver.init();
+                solver.go(Rows{});
+                for (int i = 0; i < 4; ++i) states[i] += solver.states[i];
+                terminals += solver.terminal;
+                robust += solver.robust;
+                parity += solver.parityprune;
+                zero += solver.zeroprune;
+                ++words;
+            } while (std::next_permutation(projectors.begin(), projectors.end()));
+
+            results << "{\"case\":" << case_id << ",\"words\":" << words
+                    << ",\"states\":[" << states[0] << ',' << states[1] << ','
+                    << states[2] << ',' << states[3] << "],\"terminals\":" << terminals
+                    << ",\"robust\":" << robust << ",\"parity_prunes\":" << parity
+                    << ",\"zero_frame_prunes\":" << zero << "}\n";
+            results.flush();
+            if (!results || !candidates) throw std::runtime_error("Cannot write result");
+            total_words += words;
+            total_robust += robust;
+            if ((index + 1) % 10 == 0) {
+                std::cerr << index + 1 << " / " << case_count << " cases; "
+                          << total_words << " words; robust " << total_robust << '\n';
+            }
+        }
+        std::string extra;
+        if (input >> extra) throw std::runtime_error("Trailing batch data");
+        candidates.close();
+        if (!candidates) throw std::runtime_error("Cannot close candidate file");
+        results.close();
+        if (!results) throw std::runtime_error("Cannot close result file");
+        std::cout << "{\"completed_cases\":" << case_count << ",\"words\":" << total_words
+                  << ",\"robust\":" << total_robust << "}\n";
+    } catch (const std::exception& error) {
+        std::cerr << error.what() << '\n';
+        return 1;
+    }
+    return 0;
+}
diff --git a/research/bit-three-register-lower-bound/replay.py b/research/bit-three-register-lower-bound/replay.py
new file mode 100644
index 0000000000..8e289bd7bc
--- /dev/null
+++ b/research/bit-three-register-lower-bound/replay.py
@@ -0,0 +1,98 @@
+"""Regenerate the scalar catalogue and run both complete exact searches."""
+
+from collections import Counter
+from contextlib import nullcontext
+from hashlib import sha256
+from math import factorial, prod
+from pathlib import Path
+from tempfile import TemporaryDirectory
+import argparse
+import json
+import os
+import subprocess
+import time
+
+from catalogue import catalogue
+from controls import CORE_WORDS, check_formal_controls
+
+ROOT = Path(__file__).resolve().parent
+
+
+def run(command, log):
+    with log.open('w') as output:
+        result = subprocess.run(list(map(str, command)), stdout=output, stderr=subprocess.STDOUT)
+    if result.returncode:
+        raise RuntimeError(f'Command failed: {command}\n{log.read_text()[-4000:]}')
+
+
+def read_results(path, words):
+    rows = [json.loads(line) for line in path.read_text().splitlines()]
+    by_case = {row['case']: row for row in rows}
+    assert len(rows) == len(by_case) == len(words)
+    assert set(by_case) == set(range(len(words)))
+    for case, word in enumerate(words):
+        expected_orders = factorial(8) // prod(factorial(count) for count in Counter(word).values())
+        assert by_case[case]['words'] == expected_orders
+        assert by_case[case].get('nonmonomial_diagonal', 0) == 0
+    return by_case
+
+
+def replay(directory):
+    directory.mkdir(parents=True, exist_ok=True)
+    expected = json.loads((ROOT / 'expected.json').read_text())
+    words, census = catalogue()
+    # JSON normalizes the dimension keys to strings, as in the small record.
+    census = json.loads(json.dumps(census))
+    assert all(census[key] == value for key, value in expected.items() if key in census)
+    input_path = directory / 'catalogue.txt'
+    input_path.write_text(str(len(words)) + '\n' + ''.join(
+        ' '.join(map(str, (case, *word))) + '\n' for case, word in enumerate(words)))
+    assert sha256(input_path.read_bytes()).hexdigest() == expected['catalogue_sha256']
+    formal_controls = check_formal_controls()
+    compiler = os.environ.get('CXX', 'c++')
+    for name in ('rational', 'f101'):
+        run([compiler, '-O3', '-std=c++17', ROOT / f'{name}.cpp', '-o', directory / name],
+            directory / f'{name}-compile.log')
+
+    control_input = directory / 'controls.txt'
+    control_input.write_text('2\n' + ''.join(
+        ' '.join(map(str, (case, *word))) + '\n' for case, word in enumerate(CORE_WORDS)))
+    control_output = directory / 'controls.jsonl'
+    run([directory / 'f101', control_input, control_output], directory / 'controls.log')
+    controls = read_results(control_output, CORE_WORDS)
+    assert controls[0]['robust'] > 0 and controls[1]['robust'] == 0
+
+    first = directory / 'rational-results.jsonl'
+    second = directory / 'f101-results.jsonl'
+    candidates = directory / 'candidates.txt'
+    run([directory / 'rational', '--batch', input_path, first, candidates], directory / 'rational.log')
+    run([directory / 'f101', input_path, second], directory / 'f101.log')
+    rational, finite = read_results(first, words), read_results(second, words)
+    assert candidates.read_text() == ''
+    for case in rational:
+        assert all(rational[case][key] == finite[case][key] for key in ('words', 'terminals', 'robust'))
+    counts = {'families': len(words), 'orders': sum(row['words'] for row in rational.values()),
+              'diagonal_terminals': sum(row['terminals'] for row in rational.values()),
+              'survivors': sum(row['robust'] for row in rational.values())}
+    assert all(value == expected[key] for key, value in counts.items())
+    return {**counts, 'all_per_case_counts_agree': True, 'formal_controls': formal_controls,
+            'core_positive_accepted': True, 'core_rank_two_mutation_rejected': True,
+            'catalogue_sha256': sha256(input_path.read_bytes()).hexdigest(),
+            'source_sha256': {name: sha256((ROOT / name).read_bytes()).hexdigest()
+                              for name in ('catalogue.py', 'rational.cpp', 'f101.cpp')}}
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('--work-dir', type=Path, help='Keep generated inputs, binaries, and per-case results here')
+    args = parser.parse_args()
+    context = nullcontext(args.work_dir.resolve()) if args.work_dir else TemporaryDirectory(prefix='bit-three-register-')
+    start = time.monotonic()
+    with context as directory:
+        result = replay(Path(directory))
+        result['seconds'] = time.monotonic() - start
+        print(json.dumps(result, indent=2))
+
+
+if __name__ == '__main__':
+    main()