diff --git a/paulimer/README.md b/paulimer/README.md
index 66215631..ecb7a154 100644
--- a/paulimer/README.md
+++ b/paulimer/README.md
@@ -21,6 +21,9 @@ the building blocks for stabilizer quantum mechanics and quantum error correctio
 - **Clifford Unitaries**: Efficient representation enabling fast operations
   - [`CliffordUnitary`]: O(n²) Pauli conjugation via binary symplectic matrix
   - Supports all standard Clifford gates (H, S, CNOT, etc.)
+  - Decomposition into Clifford transvections (`π/4` Pauli exponents), including a
+    strict-minimum-length variant, via [`clifford_to_transvections`] and
+    [`clifford_to_transvections_minimal`]
 
 Based on algorithms from [arXiv:2309.08676](https://arxiv.org/abs/2309.08676).
 
@@ -171,6 +174,8 @@ Key documentation:
 - [`SparsePauli`](src/pauli/sparse.rs) - Sparse Pauli representation for large systems
 - [`PauliGroup`](src/pauli_group.rs) - Subgroup operations and stabilizer groups
 - [`CliffordUnitary`](src/clifford.rs) - Clifford gates and Pauli conjugation
+- [Transvection decomposition](src/clifford/transvection.rs) - Decomposing Cliffords into `π/4`
+  Pauli exponents (`clifford_to_transvections`, `clifford_to_transvections_minimal`)
 - [Trait documentation](src/lib.rs) - `Pauli`, `Clifford`, and other core traits
 
 ## Contributing
diff --git a/paulimer/bindings/python/README.md b/paulimer/bindings/python/README.md
index 281d240c..f4cbc8c9 100644
--- a/paulimer/bindings/python/README.md
+++ b/paulimer/bindings/python/README.md
@@ -22,6 +22,11 @@ print(q * q)  # Identity
 h = paulimer.CliffordUnitary.from_name("Hadamard", [0], qubit_count=1)
 print(h.image_of(paulimer.DensePauli("X")))  # Z
 
+# Decompose a Clifford into pi/4 Pauli exponents (Clifford transvections)
+cnot = paulimer.CliffordUnitary.from_name("ControlledX", [0, 1], qubit_count=2)
+factors = cnot.to_transvections_minimal()
+print(factors)  # minimal-length list of transvection Paulis reproducing the symplectic action
+
 # Stabilizer simulation
 sim = paulimer.OutcomeCompleteSimulation(2)
 sim.apply_unitary(paulimer.UnitaryOpcode.Hadamard, [0])
@@ -32,7 +37,8 @@ sim.measure(paulimer.SparsePauli("Z0"))
 ## Features
 
 - **DensePauli / SparsePauli** - Pauli operators with phase tracking and multiplication
-- **CliffordUnitary** - Clifford gates with conjugation and composition
+- **CliffordUnitary** - Clifford gates with conjugation, composition, and decomposition into `π/4`
+  Pauli exponents (`to_transvections`, `to_transvections_minimal`)
 - **PauliGroup** - Group operations including membership testing and factorization
 - **Stabilizer Simulation** - Noiseless (OutcomeComplete, OutcomeFree, OutcomeSpecific) and noisy (Faulty) modes
 
diff --git a/paulimer/bindings/python/examples/clifford-transvection-decomposition.ipynb b/paulimer/bindings/python/examples/clifford-transvection-decomposition.ipynb
new file mode 100644
index 00000000..b2e95c8f
--- /dev/null
+++ b/paulimer/bindings/python/examples/clifford-transvection-decomposition.ipynb
@@ -0,0 +1,300 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a95d9b9a",
+   "metadata": {},
+   "source": [
+    "# Decomposing Cliffords into transvections (π/4 Pauli exponents)\n",
+    "\n",
+    "Every Clifford unitary can be written as an ordered product of **Clifford transvections** — the\n",
+    "`π/4` Pauli exponents $\\exp\\!\\big(i\\tfrac{\\pi}{4} P_v\\big)$. Conjugation by such an exponent acts on\n",
+    "Pauli operators as a **symplectic transvection**\n",
+    "\n",
+    "$$\n",
+    "x \\;\\mapsto\\; x + \\langle x, v\\rangle\\, v,\n",
+    "$$\n",
+    "\n",
+    "where $\\langle\\cdot,\\cdot\\rangle$ is the symplectic (commutation) form. `paulimer` exposes two\n",
+    "decompositions, following the transvection framework of\n",
+    "[arXiv:2102.11380](https://arxiv.org/abs/2102.11380) (Pllaha, Volanto & Tirkkonen,\n",
+    "*Decomposition of Clifford Gates*):\n",
+    "\n",
+    "- [`CliffordUnitary.to_transvections`](../paulimer.pyi) — a greedy reduction that always returns a\n",
+    "  **linear** number of factors ($O(n)$),\n",
+    "- [`CliffordUnitary.to_transvections_minimal`](../paulimer.pyi) — the **strict minimum** number of\n",
+    "  factors.\n",
+    "\n",
+    "Both reproduce the Clifford's **symplectic (conjugation) action** only; the Pauli-image signs and\n",
+    "the global phase are *not* preserved (the sign of a transvection does not change its symplectic\n",
+    "action)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "9e2b876e",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-07-05T18:09:15.236725Z",
+     "iopub.status.busy": "2026-07-05T18:09:15.236584Z",
+     "iopub.status.idle": "2026-07-05T18:09:15.240675Z",
+     "shell.execute_reply": "2026-07-05T18:09:15.239473Z"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import paulimer\n",
+    "from paulimer import CliffordUnitary, SparsePauli, DensePauli"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b42aee0d",
+   "metadata": {},
+   "source": [
+    "## A single transvection\n",
+    "\n",
+    "A `π/4` Pauli exponent *is* a Clifford transvection, so the simplest Cliffords decompose into a\n",
+    "single factor. The phase gate $S = \\exp(-i\\tfrac{\\pi}{4} Z)$ and the Hadamard are both single\n",
+    "transvections (recall the returned sign is irrelevant to the symplectic action):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "9bbcb823",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-07-05T18:09:15.242261Z",
+     "iopub.status.busy": "2026-07-05T18:09:15.242209Z",
+     "iopub.status.idle": "2026-07-05T18:09:15.244600Z",
+     "shell.execute_reply": "2026-07-05T18:09:15.243912Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "S        -> [Z]\n",
+      "Hadamard -> [-𝑖Y]\n"
+     ]
+    }
+   ],
+   "source": [
+    "s_gate = CliffordUnitary.from_name(\"SqrtZ\", [0], qubit_count=1)\n",
+    "hadamard = CliffordUnitary.from_name(\"Hadamard\", [0], qubit_count=1)\n",
+    "\n",
+    "print(\"S        ->\", s_gate.to_transvections_minimal())\n",
+    "print(\"Hadamard ->\", hadamard.to_transvections_minimal())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c78b1510",
+   "metadata": {},
+   "source": [
+    "## Rebuilding a Clifford and checking the symplectic action\n",
+    "\n",
+    "Applying the returned transvections in order with\n",
+    "[`left_mul_pauli_exp`](../paulimer.pyi) reconstructs the original **symplectic matrix**. We compare\n",
+    "`symplectic_matrix` (not the full signed tableau, since signs and global phase are not tracked by\n",
+    "this decomposition).\n",
+    "\n",
+    "The minimal factor count is either $r$ or $r+1$, where the **residue rank**\n",
+    "\n",
+    "$$\n",
+    "r \\;=\\; 2n - \\dim \\operatorname{Fix}(F)\n",
+    "$$\n",
+    "\n",
+    "is the codimension of the space of Pauli operators fixed under conjugation. In `paulimer`,\n",
+    "$\\dim\\operatorname{Fix}(F)$ is the size of the Clifford's centralizer."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "ed66171a",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-07-05T18:09:15.246310Z",
+     "iopub.status.busy": "2026-07-05T18:09:15.246264Z",
+     "iopub.status.idle": "2026-07-05T18:09:15.248826Z",
+     "shell.execute_reply": "2026-07-05T18:09:15.248499Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "residue rank r      = 1\n",
+      "number of factors   = 1\n",
+      "symplectic action ok: True\n"
+     ]
+    }
+   ],
+   "source": [
+    "def residue_rank(clifford):\n",
+    "    return 2 * clifford.qubit_count - len(clifford.centralizer())\n",
+    "\n",
+    "\n",
+    "def rebuild(factors, qubit_count):\n",
+    "    rebuilt = CliffordUnitary.identity(qubit_count)\n",
+    "    for pauli in factors:\n",
+    "        rebuilt.left_mul_pauli_exp(pauli)\n",
+    "    return rebuilt\n",
+    "\n",
+    "\n",
+    "factors = s_gate.to_transvections_minimal()\n",
+    "rebuilt = rebuild(factors, s_gate.qubit_count)\n",
+    "print(\"residue rank r      =\", residue_rank(s_gate))\n",
+    "print(\"number of factors   =\", len(factors))\n",
+    "print(\"symplectic action ok:\", rebuilt.symplectic_matrix == s_gate.symplectic_matrix)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eab0ecb9",
+   "metadata": {},
+   "source": [
+    "## Greedy versus minimal, and the $r+1$ case\n",
+    "\n",
+    "For many Cliffords the greedy and minimal decompositions agree, but not always. The CNOT gate has\n",
+    "residue rank $r = 2$ yet needs $r + 1 = 3$ transvections: its symplectic action is *hyperbolic*\n",
+    "($\\langle v, vF\\rangle = 0$ for all $v$), which forces one extra factor."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "id": "00096f40",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-07-05T18:09:15.250074Z",
+     "iopub.status.busy": "2026-07-05T18:09:15.249935Z",
+     "iopub.status.idle": "2026-07-05T18:09:15.252141Z",
+     "shell.execute_reply": "2026-07-05T18:09:15.251741Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "residue rank r = 2\n",
+      "greedy  : [Z, ZX, IX]  ( 3 factors )\n",
+      "minimal : [Z, ZX, IX]  ( 3 factors )\n",
+      "symplectic action ok: True\n"
+     ]
+    }
+   ],
+   "source": [
+    "cnot = CliffordUnitary.from_name(\"ControlledX\", [0, 1], qubit_count=2)\n",
+    "\n",
+    "greedy = cnot.to_transvections()\n",
+    "minimal = cnot.to_transvections_minimal()\n",
+    "print(\"residue rank r =\", residue_rank(cnot))\n",
+    "print(\"greedy  :\", greedy, \" (\", len(greedy), \"factors )\")\n",
+    "print(\"minimal :\", minimal, \" (\", len(minimal), \"factors )\")\n",
+    "print(\"symplectic action ok:\", rebuild(minimal, 2).symplectic_matrix == cnot.symplectic_matrix)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6d7dfd52",
+   "metadata": {},
+   "source": [
+    "## A subtle case: non-hyperbolic maps that still need $r+1$\n",
+    "\n",
+    "The 2021 paper claims that *every* non-hyperbolic Clifford decomposes into exactly $r$ transvections.\n",
+    "That is **not correct over $\\mathbb{F}_2$**: some non-hyperbolic maps still require $r + 1$. The\n",
+    "smallest example already occurs on two qubits — the symplectic action built below (a product of the\n",
+    "transvections $X_0, X_1, X_0X_1, Z_0$) has residue rank $r = 3$, is non-hyperbolic, yet needs $4$\n",
+    "transvections. `to_transvections_minimal` returns the correct minimum. See\n",
+    "[`docs/transvection-minimality-correction.md`](../../../docs/transvection-minimality-correction.md)\n",
+    "for the full analysis and a machine-checked proof."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "3a68cc03",
+   "metadata": {
+    "execution": {
+     "iopub.execute_input": "2026-07-05T18:09:15.253271Z",
+     "iopub.status.busy": "2026-07-05T18:09:15.253223Z",
+     "iopub.status.idle": "2026-07-05T18:09:15.256499Z",
+     "shell.execute_reply": "2026-07-05T18:09:15.255166Z"
+    }
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "residue rank r      = 3\n",
+      "minimal factors     = [IX, XX, -𝑖Y, X] ( 4 factors )\n",
+      "needs r + 1          : True\n",
+      "symplectic action ok: True\n"
+     ]
+    }
+   ],
+   "source": [
+    "example = CliffordUnitary.identity(2)\n",
+    "for pauli in [\"X0\", \"X1\", \"X0 X1\", \"Z0\"]:\n",
+    "    example.left_mul_pauli_exp(SparsePauli(pauli))\n",
+    "\n",
+    "minimal = example.to_transvections_minimal()\n",
+    "r = residue_rank(example)\n",
+    "print(\"residue rank r      =\", r)\n",
+    "print(\"minimal factors     =\", minimal, \"(\", len(minimal), \"factors )\")\n",
+    "print(\"needs r + 1          :\", len(minimal) == r + 1)\n",
+    "print(\"symplectic action ok:\", rebuild(minimal, 2).symplectic_matrix == example.symplectic_matrix)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "992a2da1",
+   "metadata": {},
+   "source": [
+    "## Summary\n",
+    "\n",
+    "- Clifford transvections are `π/4` Pauli exponents; `to_transvections` /\n",
+    "  `to_transvections_minimal` decompose any Clifford into them, reproducing its symplectic action\n",
+    "  with $O(n)$ factors.\n",
+    "- The minimal count is $r$ or $r + 1$, where $r = 2n - \\dim\\operatorname{Fix}(F)$.\n",
+    "- Only the symplectic action is reproduced — Pauli-image signs and the global phase are not.\n",
+    "\n",
+    "### References\n",
+    "\n",
+    "- T. Pllaha, K. Volanto, O. Tirkkonen, *Decomposition of Clifford Gates*, GLOBECOM 2021,\n",
+    "  [arXiv:2102.11380](https://arxiv.org/abs/2102.11380).\n",
+    "- [`docs/transvection-minimality-correction.md`](../../../docs/transvection-minimality-correction.md)\n",
+    "  — a correction to the paper's minimality claim, with a verified counterexample."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "paulimer",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/paulimer/bindings/python/paulimer.pyi b/paulimer/bindings/python/paulimer.pyi
index ee79347d..6346ca79 100644
--- a/paulimer/bindings/python/paulimer.pyi
+++ b/paulimer/bindings/python/paulimer.pyi
@@ -582,6 +582,32 @@ class CliffordUnitary:
         """Get the symplectic matrix representation."""
         ...
 
+    def to_transvections(self) -> list[SparsePauli]:
+        """Decompose into an ordered product of Clifford transvections (pi/4 Pauli exponents).
+
+        Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+        ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces this Clifford's symplectic
+        (conjugation) action, using a linear number of factors. Pauli-image signs and the global
+        phase are not reproduced.
+        """
+        ...
+
+    def to_transvections_minimal(self) -> list[SparsePauli]:
+        """Decompose into a *minimal* ordered product of Clifford transvections (pi/4 Pauli exponents).
+
+        Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+        ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces this Clifford's symplectic
+        (conjugation) action, with ``k`` the minimal transvection count (``r`` or ``r + 1``, where
+        ``r`` is the rank of the residue matrix). Pauli-image signs and the global phase are not
+        reproduced; :meth:`to_transvections` is the linear-time greedy variant, which may use more
+        factors.
+        """
+        ...
+
+    def centralizer(self) -> list[SparsePauli]:
+        """Generators of the centralizer: Paulis fixed up to sign under conjugation."""
+        ...
+
     def qubits(self) -> slice:
         """Return a slice representing the qubit indices."""
         ...
diff --git a/paulimer/bindings/python/src/py_clifford.rs b/paulimer/bindings/python/src/py_clifford.rs
index bdf51d2a..35836c17 100644
--- a/paulimer/bindings/python/src/py_clifford.rs
+++ b/paulimer/bindings/python/src/py_clifford.rs
@@ -1,7 +1,7 @@
 use derive_more::{Deref, DerefMut, From, Into};
 use paulimer::clifford::{
-    group_encoding_clifford_of, split_phased_css, split_qubit_cliffords_and_css, Clifford, CliffordMutable,
-    CliffordUnitary, XOrZ,
+    clifford_centralizer, clifford_to_transvections, clifford_to_transvections_minimal, group_encoding_clifford_of,
+    split_phased_css, split_qubit_cliffords_and_css, Clifford, CliffordMutable, CliffordUnitary, XOrZ,
 };
 use paulimer::pauli::{as_sparse, DensePauli, SparsePauli};
 use pyo3::exceptions::PyValueError;
@@ -290,6 +290,45 @@ impl PyCliffordUnitary {
         self.inner.symplectic_matrix().into()
     }
 
+    /// Decomposes this Clifford into an ordered product of Clifford transvections (pi/4 Pauli
+    /// exponents), reproducing its symplectic action with a linear number of factors.
+    ///
+    /// Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+    /// ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces the conjugation action of this
+    /// Clifford. Pauli-image signs and the global phase are not reproduced; a sign-exact
+    /// decomposition into Pauli exponents would preserve them, at the cost of ``O(n^2)`` factors.
+    fn to_transvections(&self) -> Vec<PySparsePauli> {
+        clifford_to_transvections(&self.inner)
+            .into_iter()
+            .map(PySparsePauli::from)
+            .collect()
+    }
+
+    /// Decomposes this Clifford into a *minimal* ordered product of Clifford transvections (pi/4
+    /// Pauli exponents), reproducing its symplectic action with the fewest possible factors.
+    ///
+    /// Returns Pauli operators ``[P_1, ..., P_k]`` such that applying ``exp(i pi/4 P_1)``, then
+    /// ``exp(i pi/4 P_2)``, ..., then ``exp(i pi/4 P_k)`` reproduces the conjugation action of this
+    /// Clifford, with ``k`` equal to the minimal transvection count (``r`` or ``r + 1``, where ``r``
+    /// is the rank of the residue matrix). Pauli-image signs and the global phase are not
+    /// reproduced; see :meth:`to_transvections` for the linear-time greedy decomposition, which may
+    /// use more factors.
+    fn to_transvections_minimal(&self) -> Vec<PySparsePauli> {
+        clifford_to_transvections_minimal(&self.inner)
+            .into_iter()
+            .map(PySparsePauli::from)
+            .collect()
+    }
+
+    /// Returns generators of this Clifford's centralizer: the Pauli operators fixed up to sign under
+    /// conjugation (``clifford * P * clifford_dagger == +/- P``).
+    fn centralizer(&self) -> Vec<PySparsePauli> {
+        clifford_centralizer(&self.inner)
+            .into_iter()
+            .map(PySparsePauli::from)
+            .collect()
+    }
+
     #[allow(clippy::needless_pass_by_value)]
     fn left_mul(&mut self, unitary_op: PyUnitaryOp, support: Vec<usize>) {
         self.inner.left_mul(unitary_op.into(), &support);
diff --git a/paulimer/bindings/python/tests/transvection_test.py b/paulimer/bindings/python/tests/transvection_test.py
new file mode 100644
index 00000000..829ab9af
--- /dev/null
+++ b/paulimer/bindings/python/tests/transvection_test.py
@@ -0,0 +1,149 @@
+"""Tests for the Clifford -> transvection decomposition bindings (arXiv:2102.11380).
+
+The decomposition reproduces a Clifford's symplectic (conjugation) action with a linear number of
+pi/4 Pauli exponents, ignoring Pauli-image signs and the global phase.
+"""
+
+from hypothesis import given, settings
+from hypothesis import strategies as st
+
+from paulimer import CliffordUnitary, SparsePauli, UnitaryOpcode
+
+
+def _rebuild_from_transvections(transvections, qubit_count):
+    rebuilt = CliffordUnitary.identity(qubit_count)
+    for pauli in transvections:
+        rebuilt.left_mul_pauli_exp(pauli)
+    return rebuilt
+
+
+def _residue_rank(clifford):
+    return 2 * clifford.qubit_count - len(clifford.centralizer())
+
+
+def _is_conjugation_fixed(clifford, pauli):
+    image = SparsePauli.from_dense(clifford.image_of(pauli))
+    return (image * pauli).weight == 0
+
+
+def _assert_valid_decomposition(clifford):
+    qubit_count = clifford.qubit_count
+    transvections = clifford.to_transvections()
+
+    rebuilt = _rebuild_from_transvections(transvections, qubit_count)
+    assert rebuilt.symplectic_matrix == clifford.symplectic_matrix
+
+    for pauli in transvections:
+        assert pauli.weight > 0
+
+    minimum = _residue_rank(clifford)
+    assert len(transvections) >= minimum
+    assert len(transvections) <= 4 * qubit_count + 2
+
+
+def _assert_valid_minimal_decomposition(clifford):
+    qubit_count = clifford.qubit_count
+    transvections = clifford.to_transvections_minimal()
+
+    rebuilt = _rebuild_from_transvections(transvections, qubit_count)
+    assert rebuilt.symplectic_matrix == clifford.symplectic_matrix
+
+    for pauli in transvections:
+        assert pauli.weight > 0
+
+    rank = _residue_rank(clifford)
+    assert len(transvections) in (rank, rank + 1)
+    assert len(transvections) <= len(clifford.to_transvections())
+
+
+def test_identity_has_no_transvections():
+    for qubit_count in range(5):
+        identity = CliffordUnitary.identity(qubit_count)
+        assert identity.to_transvections() == []
+        assert identity.to_transvections_minimal() == []
+        assert len(identity.centralizer()) == 2 * qubit_count
+
+
+def test_single_qubit_gate_lengths():
+    s_gate = CliffordUnitary.from_name("SqrtZ", [0], 1)
+    _assert_valid_decomposition(s_gate)
+    _assert_valid_minimal_decomposition(s_gate)
+    assert len(s_gate.to_transvections()) == 1
+    assert len(s_gate.to_transvections_minimal()) == 1
+
+    hadamard = CliffordUnitary.from_name("Hadamard", [0], 1)
+    _assert_valid_decomposition(hadamard)
+    _assert_valid_minimal_decomposition(hadamard)
+    assert len(hadamard.to_transvections()) == 1
+    assert len(hadamard.to_transvections_minimal()) == 1
+
+
+def test_swap_hyperbolic_branch():
+    swap = CliffordUnitary.from_name("Swap", [0, 1], 2)
+    _assert_valid_decomposition(swap)
+    _assert_valid_minimal_decomposition(swap)
+    assert _residue_rank(swap) == 2
+    assert len(swap.to_transvections()) == 3
+    assert len(swap.to_transvections_minimal()) == 3
+    assert len(swap.centralizer()) == 2
+
+
+def test_two_qubit_gates():
+    for name in ("ControlledX", "ControlledZ"):
+        clifford = CliffordUnitary.from_name(name, [0, 1], 2)
+        _assert_valid_decomposition(clifford)
+        _assert_valid_minimal_decomposition(clifford)
+
+
+def test_centralizer_generators_are_conjugation_fixed():
+    clifford = CliffordUnitary.identity(3)
+    clifford.left_mul(UnitaryOpcode.Hadamard, [0])
+    clifford.left_mul(UnitaryOpcode.ControlledX, [0, 1])
+    clifford.left_mul(UnitaryOpcode.SqrtZ, [2])
+    centralizer = clifford.centralizer()
+    assert all(_is_conjugation_fixed(clifford, pauli) for pauli in centralizer)
+    assert all(pauli.weight > 0 for pauli in centralizer)
+
+
+_SINGLE_QUBIT_GATES = ["Hadamard", "SqrtZ", "SqrtX", "X", "Y", "Z"]
+_TWO_QUBIT_GATES = ["ControlledX", "ControlledZ", "Swap"]
+
+
+@st.composite
+def _random_clifford(draw, max_qubits=5):
+    qubit_count = draw(st.integers(min_value=1, max_value=max_qubits))
+    gate_count = draw(st.integers(min_value=0, max_value=3 * qubit_count))
+    clifford = CliffordUnitary.identity(qubit_count)
+    for _ in range(gate_count):
+        if qubit_count >= 2 and draw(st.booleans()):
+            name = draw(st.sampled_from(_TWO_QUBIT_GATES))
+            first = draw(st.integers(min_value=0, max_value=qubit_count - 1))
+            second = draw(
+                st.integers(min_value=0, max_value=qubit_count - 1).filter(lambda q: q != first)
+            )
+            clifford.left_mul(getattr(UnitaryOpcode, name), [first, second])
+        else:
+            name = draw(st.sampled_from(_SINGLE_QUBIT_GATES))
+            qubit = draw(st.integers(min_value=0, max_value=qubit_count - 1))
+            clifford.left_mul(getattr(UnitaryOpcode, name), [qubit])
+    return clifford
+
+
+@settings(max_examples=200)
+@given(_random_clifford())
+def test_random_cliffords_reproduce_symplectic_action(clifford):
+    _assert_valid_decomposition(clifford)
+
+
+@settings(max_examples=200)
+@given(_random_clifford())
+def test_random_cliffords_minimal_reproduce_symplectic_action(clifford):
+    _assert_valid_minimal_decomposition(clifford)
+
+
+@settings(max_examples=200)
+@given(_random_clifford())
+def test_random_centralizers_are_conjugation_fixed(clifford):
+    for generator in clifford.centralizer():
+        assert _is_conjugation_fixed(clifford, generator)
+        assert generator.weight > 0
diff --git a/paulimer/docs/transvection-minimality-correction.md b/paulimer/docs/transvection-minimality-correction.md
new file mode 100644
index 00000000..7f12fb9a
--- /dev/null
+++ b/paulimer/docs/transvection-minimality-correction.md
@@ -0,0 +1,339 @@
+# Reassessment of the minimal transvection decomposition in arXiv:2102.11380
+
+This note reassesses the minimal-transvection construction in
+
+> T. Pllaha, K. Volanto, and O. Tirkkonen,
+> *Decomposition of Clifford Gates*,
+> 2021 IEEE Global Communications Conference (GLOBECOM), 2021, pp. 1–6.
+> DOI [10.1109/GLOBECOM46510.2021.9685501](https://doi.org/10.1109/GLOBECOM46510.2021.9685501),
+> arXiv:[2102.11380](https://arxiv.org/abs/2102.11380).
+
+against its cited primary source, Callan (1976), and an exhaustive two-qubit
+calculation.
+
+The paper's structural matrix identities are correct, and its top-level
+existence claim -- that a minimal decomposition algorithm exists -- is true.
+However, the construction given in and immediately before **Theorem 3** assumes
+that every non-hyperbolic binary symplectic map has a length-$r$
+decomposition. Callan explicitly classifies non-hyperbolic exceptions, and the
+assumption fails on two qubits. The replacement proof was checked independently
+with an interactive theorem prover. It repairs this step; it does not prove the
+paper's non-hyperbolic criterion.
+
+All matrices below are over $\mathbb{F}_2$.
+
+## 1. Summary
+
+The paper decomposes a symplectic map $\mathbf{F}\in\mathrm{Sp}(2m;2)$ into
+*symplectic transvections* and claims in the discussion preceding Theorem 3
+that the number of factors equals the **residue rank**
+
+$$
+r \;=\; \dim\operatorname{Res}(\mathbf F) \;=\; 2m-\dim\operatorname{Fix}(\mathbf F)
+$$
+
+whenever $\mathbf F$ is **non-hyperbolic**, and $r+1$ when $\mathbf F$
+is hyperbolic. Concretely, line 411 of the arXiv v1 source asserts that the
+residue core **can always be triangularized by congruence for a
+non-hyperbolic map**.
+
+That assertion is false over $\mathbb{F}_2$. Callan defines an element as
+*exceptional* precisely when it cannot be expressed using $r$ transvections.
+His §1.1 proves the universal bounds
+
+$$
+r\leq\ell(\mathbf F)\leq r+1,
+$$
+
+and Theorem 5.1 classifies the binary exceptions. They include
+non-hyperbolic maps. Thus the part that fails is the
+hyperbolic/non-hyperbolic classification, not the universal $r$/$r+1$ range.
+The smallest non-hyperbolic example already occurs on two qubits, with
+$r=3$ and $\ell(\mathbf F)=4$.
+
+The **correct criterion**, which we adopt in the implementation, is:
+
+> The minimal length is $r$ **iff** the invertible residue core $\mathbf E$ is
+> congruence-lower-triangularizable over $\mathbb{F}_2$; otherwise our construction returns a
+> decomposition of length $r+1$. Hyperbolicity ($\mathbf E$ *alternating*) is a special
+> $r+1$ sub-case, but it is **not** the only one: non-alternating cores can fail to be
+> triangularizable too.
+
+That independent check covers the bound, the criterion above, strict
+minimality, and the one-fix theorem used by the implementation.
+
+The Rust implementation should therefore retain its complete congruence search
+and $r+1$ fallback. No semantic rollback to the paper's non-hyperbolic branch is
+warranted.
+
+## 2. Setup and notation
+
+We follow the paper's conventions. Pauli operators on $m$ qubits are represented by row vectors
+$\mathbf v\in\mathbb{F}_2^{2m}$; a Clifford acts on them by a symplectic matrix
+$\mathbf F\in\mathrm{Sp}(2m;2)$ via the right action $\mathbf x\mapsto\mathbf x\mathbf F$. With
+$\boldsymbol\Omega=\left(\begin{smallmatrix}\mathbf 0&\mathbf I\\\mathbf I&\mathbf 0\end{smallmatrix}\right)$
+the symplectic form is $\langle\mathbf u,\mathbf v\rangle=\mathbf u\,\boldsymbol\Omega\,\mathbf v^{\mathsf T}$.
+
+A **symplectic transvection** is the map
+
+$$
+\mathbf T_{\mathbf v}\;=\;\mathbf I+\boldsymbol\Omega\,\mathbf v^{\mathsf T}\mathbf v,
+\qquad\text{i.e.}\qquad
+\mathbf x\,\mathbf T_{\mathbf v}=\mathbf x+\langle\mathbf x,\mathbf v\rangle\,\mathbf v ,
+$$
+
+the conjugation action of the Clifford transvection $\exp(i\tfrac{\pi}{4}P_{\mathbf v})$. It is
+classical that $\mathrm{Sp}(2m;2)$ is generated by transvections. The **fixed** and **residue**
+spaces are
+
+$$
+\operatorname{Fix}(\mathbf F)=\ker(\mathbf I+\mathbf F),\qquad
+\operatorname{Res}(\mathbf F)=\operatorname{rowsp}(\mathbf I+\mathbf F),\qquad
+r:=\dim\operatorname{Res}(\mathbf F)=\operatorname{rank}(\mathbf I+\mathbf F).
+$$
+
+$\mathbf F$ is **hyperbolic** iff $\langle\mathbf v,\mathbf v\mathbf F\rangle=0$ for all $\mathbf v$.
+The **residue matrix** is
+
+$$
+\widehat{\mathbf F}:=\boldsymbol\Omega(\mathbf I+\mathbf F),\qquad
+\operatorname{rowsp}(\widehat{\mathbf F})=\operatorname{Res}(\mathbf F),\qquad
+\operatorname{rank}(\widehat{\mathbf F})=r,
+$$
+
+The correct matrix test is that $\mathbf F$ is hyperbolic iff
+$\widehat{\mathbf F}$ is *alternating*, meaning symmetric with zero diagonal.
+Zero diagonal alone is necessary but not sufficient: a non-involution can have
+a nonsymmetric residue matrix with zero diagonal. Row-reducing
+$\widehat{\mathbf F}$ with a transform $\mathbf R$ yields the invertible
+**core**
+
+$$
+\mathbf R\,\widehat{\mathbf F}\,\mathbf R^{\mathsf T}=\begin{pmatrix}\mathbf E&\mathbf 0\\\mathbf 0&\mathbf 0\end{pmatrix},
+\qquad \mathbf E\in\mathrm{GL}(r;2).
+$$
+
+Write $\psi_{\mathbf E}(\mathbf x)=\mathbf x\,\mathbf E\,\mathbf x^{\mathsf T}$ for the associated
+quadratic form; $\mathbf E$ is *alternating* iff $\psi_{\mathbf E}\equiv 0$ iff $\mathbf F$ is
+hyperbolic.
+
+## 3. The paper's claim
+
+The paper's two key lemmas are correct and we use them:
+
+- **Lemma 2 (l-TET).** A length-$r$ basis $\mathbf Q\mathbf V$ of $\operatorname{Res}(\mathbf F)$
+  (with $\mathbf Q\in\mathrm{GL}(r;2)$) constitutes a transvection decomposition of $\mathbf F$ **iff**
+  $\mathbf Q\mathbf E\mathbf Q^{\mathsf T}=\mathbf B^{-\mathsf T}$, where $\mathbf B$ is the paper's
+  upper-triangular, unit-diagonal path-counting matrix.
+- **Lemma 3 (l-QEQ).** If $\mathbf Q\mathbf E\mathbf Q^{\mathsf T}$ is lower-triangular, then it
+  automatically equals $\mathbf B^{-\mathsf T}$. (Over $\mathbb{F}_2$ an invertible triangular matrix
+  necessarily has unit diagonal, so no separate diagonal condition is needed.)
+
+Together these give the correct reduction: **a length-$r$ transvection decomposition of $\mathbf F$
+exists iff there is $\mathbf Q\in\mathrm{GL}(r;2)$ making $\mathbf Q\mathbf E\mathbf Q^{\mathsf T}$
+lower-triangular** — i.e. iff $\mathbf E$ is congruence-triangularizable. So far, so good.
+
+An equivalent faithful residue presentation writes
+$\widehat{\mathbf F}=\mathbf V^{\mathsf T}\mathbf D\mathbf V$ and defines
+its core as $\mathbf D^{-\mathsf T}$. From the paper's row-reduction
+identity, $\mathbf D=\mathbf E^{-\mathsf T}$, so this core is exactly
+the paper's $\mathbf E$. The difference is notation, not a transpose or
+action convention.
+
+The error is the very next sentence (the paragraph following Lemmas 2–3; line 411 of the arXiv v1
+source), which asserts existence unconditionally:
+
+> "It also follows … that $\widehat{\mathbf F}$ *can* be triangularized by congruence for any
+> non-hyperbolic $\mathbf F$ (since for this, one would only need a transvection decomposition of
+> $\mathbf F$, which we know it always exists)."
+
+and the earlier statement (in §III, following the transvection definition; line 174 of the source)
+attributed to O'Meara and Callan:
+
+> "a non-hyperbolic map $\mathbf F$ can be written as a product of $r$ *independent* transvections."
+
+Theorem 3 (T-main1) then instructs one to "let $\mathbf Q$ be such that
+$\mathbf Q\mathbf E\mathbf Q^{\mathsf T}$ is lower triangular", assuming such $\mathbf Q$ exists for
+every non-hyperbolic $\mathbf F$.
+
+**What the cited literature actually says.** O'Meara's Theorem 2.1.11
+gives the simple non-hyperbolic/hyperbolic dichotomy under the hypothesis
+$F\neq\mathbb F_2$. His 2.1.18 extends it to involutions in characteristic
+two, including $\mathbb F_2$, while 2.1.19 again excludes $\mathbb F_2$ for
+general maps. The §2.3 comment that the binary theorem is "considerably more
+complicated" refers to the exceptional-class classification, not to failure
+of the $r+1$ upper bound.
+
+Callan is decisive because the paper cites it directly. Callan §1.1 proves
+that every binary symplectic map has length at most $r+1$, §2.4 identifies
+the residue-$3$ class-A exceptions in dimension $4$, and Theorem 5.1 gives
+the complete exceptional list. The paper overlooks those exceptions.
+
+**The flaw in the paper's own justification.** Independently of the incorrect attribution, the line-411
+argument is circular: "a transvection decomposition always exists" is true (transvections generate
+the group), but it only guarantees *some* decomposition — possibly of length $r+1$. A
+length-$(r+1)$ decomposition does **not** correspond to any $\mathbf Q\in\mathrm{GL}(r;2)$
+triangularizing the $r\times r$ core (Lemma 2 is specifically about length-$r$ bases of
+$\operatorname{Res}(\mathbf F)$). The argument conflates *existence of a decomposition* with
+*existence of a minimal, length-$r$ one*. Concretely, the greedy strategy — repeatedly pick
+$\mathbf x$ with $\langle\mathbf x,\mathbf x\mathbf F\rangle=1$ and reduce to
+$\mathbf F\mathbf T_{\mathbf v}$ with $r(\mathbf F\mathbf T_{\mathbf v})=r-1$ — can make an
+intermediate map **hyperbolic** before the residue is exhausted, at which point Lemma 1 (l-hyp) must
+spend an *extra* transvection, yielding $r+1$ overall. Over $\mathbb{F}_2$, non-hyperbolicity of the
+*initial* map does not prevent this.
+
+## 4. The gap: non-hyperbolic does **not** imply triangularizable
+
+Take $m=2$ qubits, coordinates $(x_0,x_1,z_0,z_1)$. The symplectic matrix (acting on the right)
+
+$$
+\mathbf F=\begin{pmatrix}1&0&1&0\\0&1&0&0\\0&1&1&0\\1&0&1&1\end{pmatrix}
+$$
+
+satisfies:
+
+- **Symplectic and non-hyperbolic.** $\mathbf F\in\mathrm{Sp}(4;2)$, and
+  $\langle\mathbf e_0,\mathbf e_0\mathbf F\rangle=1$, so $\mathbf F$
+  is not hyperbolic. The paper would therefore predict minimal length
+  $r=3$.
+- **Residue rank $r=3$.** $\operatorname{rank}(\mathbf I+\mathbf F)=3$.
+- **Residue matrix and core.**
+
+$$
+\widehat{\mathbf F}=\boldsymbol\Omega(\mathbf I+\mathbf F)=
+\begin{pmatrix}0&1&0&0\\1&0&1&0\\0&0&1&0\\0&0&0&0\end{pmatrix},
+\qquad
+\mathbf E=\begin{pmatrix}0&1&0\\1&0&0\\1&0&1\end{pmatrix}\in\mathrm{GL}(3;2).
+$$
+
+  The diagonal of $\widehat{\mathbf F}$ is nonzero, and
+  $\psi_{\mathbf E}(\mathbf x)=\mathbf x\mathbf E\mathbf x^{\mathsf T}\not\equiv 0$, confirming
+  $\mathbf E$ is **non-alternating** (again: non-hyperbolic).
+
+Two independent exhaustive computations refute the paper's claim for this $\mathbf F$:
+
+1. **Minimal transvection length is $4=r+1$.** Breadth-first search over the *entire* group
+   $\mathrm{Sp}(4;2)$ (720 elements) with the 15 transvections as generators gives Cayley distance
+   $\ell(\mathbf F)=4$. There is no product of three transvections equal to $\mathbf F$.
+2. **The core $\mathbf E$ is not congruence-triangularizable.** Enumerating all $168$ elements of
+   $\mathrm{GL}(3;2)$, **no** $\mathbf Q$ makes $\mathbf Q\mathbf E\mathbf Q^{\mathsf T}$
+   lower-triangular. So the length-$r$ criterion of Section 3 fails, consistent with (1).
+
+Because $\mathbf F$ is non-hyperbolic yet requires $r+1$ transvections, the sentence at line 411 —
+and Theorem 3's assumption that a triangularizing $\mathbf Q$ always exists — are false.
+
+This is a concrete instance of Callan's class A, not a new exceptional
+family. Callan §2.4 characterizes the residue-$3$ exceptions in
+$\mathrm{Sp}(4;2)$ as class A, and Theorem 5.1 includes class A in the
+complete binary exception list.
+
+**Cross-check against this repository's code.** Building the same map in `paulimer` (as the product
+of the four transvections $X_0,\,X_1,\,X_0X_1,\,Z_0$ found by the search) and calling the shipped
+decomposer reproduces exactly this behaviour:
+
+```text
+residue rank r         = 3          # 2*qubit_count - |centralizer|
+centralizer size       = 1          # = 2m - r
+to_transvections       len = 4
+to_transvections_minimal len = 4    # = r + 1, not r
+```
+
+**The phenomenon is common, not isolated.** An exhaustive census of
+$\mathrm{Sp}(4;2)$ (all $720$ elements) shows that $225$ of the $719$ non-identity maps require
+$r+1$ transvections — and of those $225$, **$210$ are non-hyperbolic** and only $15$ are hyperbolic.
+So over two qubits the paper's rule "non-hyperbolic $\Rightarrow$ length $r$" is violated by $210$
+distinct symplectic maps. The phenomenon first appears at $r=3$; the map above is one representative.
+(The same census confirms the minimal length never exceeds $r+1$ for $m=2$, so hyperbolicity is the
+*wrong* invariant, not the count $r+1$ itself — at least at this size.)
+
+## 5. Implementation note: triangularization requires a complete decision procedure
+
+The paper does not specify a greedy triangularization algorithm; it refers to
+Botha's work. Therefore a dead-ending greedy pivot choice is not a second
+error in the paper. It is an implementation pitfall.
+
+The Rust implementation explores every non-isotropic pivot and memoizes
+subspaces already proved unsolvable, so its decision procedure does not depend
+on a greedy choice. The earlier version of this note cited a specific
+proptest regression file that is no longer present; that historical claim is
+not needed for the paper counterexample or the correctness argument here.
+
+## 6. Corrected result and implementation
+
+Combining the correct Lemmas 2–3, Callan's $r+1$ bound, and the bordered
+construction gives:
+
+$$
+\ell(\mathbf F)=
+\begin{cases}
+r & \text{if } \mathbf E \text{ is congruence-lower-triangularizable over } \mathbb{F}_2,\\[2pt]
+r+1 & \text{otherwise (this includes, but is strictly larger than, the hyperbolic case).}
+\end{cases}
+$$
+
+Implementation ([`transvection.rs`](../src/clifford/transvection.rs)):
+
+- **Triangularization by complete search.** `congruence_triangularize` performs an exhaustive
+  backtracking search for $\mathbf Q$: at each node it enumerates all $\psi$-non-isotropic pivots,
+  recurses into the right-orthogonal complement, and **memoizes subspaces proven untriangularizable**
+  by a canonical row-reduced key.
+- **The $r+1$ fix vector.** When (and only when) no $\mathbf Q$ exists, `find_fix_vector` appends one
+  extra transvection $\mathbf T_{\mathbf w}$ chosen so that the residue-preserving update
+  $\mathbf F\mathbf T_{\mathbf w}$ *becomes* triangularizable at the same rank, and recurses. This is
+  the non-hyperbolic analogue of the paper's hyperbolic Lemma 1 patch — the case the paper's
+  construction omits.
+
+The test suite additionally checks the result against a brute-force BFS oracle on one and two
+qubits and against the $\{r,r+1\}$ range on up to six qubits. These computations are regression
+checks, not the justification for generality. The proof establishes for every finite $m$ that
+$r\leq\ell(\mathbf F)\leq r+1$, that length $r$ is equivalent to core triangularizability, and that
+otherwise a nonzero $\mathbf w\in\operatorname{Res}(\mathbf F)$ exists for which
+$\mathbf F\mathbf T_{\mathbf w}$ has the same residue rank and a triangularizable core. Thus the
+exhaustive `find_fix_vector` search is total on valid symplectic input.
+
+**Why "non-alternating" is not enough.** Alternating cores are
+untriangularizable, but the explicit core in Section 4 is non-alternating and
+still untriangularizable. The exact condition is congruence
+triangularizability of the full core, not merely its diagonal or associated
+quadratic form. Botha (1997), which the paper cites, studies this GF(2)
+congruence problem directly.
+
+## 7. Reproducing the verification
+
+The counterexample of Section 4 is fully finite and self-contained. Both checks — the
+$\mathrm{Sp}(4;2)$ Cayley-distance BFS (720 group elements) and the $\mathrm{GL}(3;2)$ congruence
+enumeration (168 candidates) — are small enough to run by hand or in a few lines of code, and the
+repository's own `clifford_to_transvections_minimal` reproduces $\ell(\mathbf F)=r+1$ on the same
+map. No floating point or randomness is involved. The symbolic proof is reproduced with
+`cd paulimer/formal && lake build`; it contains no admitted theorem or project-defined axiom.
+
+## 8. References
+
+Bibliographic details are taken from the paper's reference list and the cited
+primary sources.
+
+1. T. Pllaha, K. Volanto, O. Tirkkonen. *Decomposition of Clifford Gates.* 2021 IEEE Global
+   Communications Conference (GLOBECOM), 2021.
+   DOI [10.1109/GLOBECOM46510.2021.9685501](https://doi.org/10.1109/GLOBECOM46510.2021.9685501);
+   arXiv:[2102.11380](https://arxiv.org/abs/2102.11380). *(The paper corrected here.)*
+2. O. T. O'Meara. *Symplectic Groups.* Mathematical Surveys, vol. 16. American Mathematical Society,
+   Providence, R.I., 1978. *(Transvection generation and the greedy residue reduction. Theorem
+   2.1.11 gives the $r/(r+1)$ dichotomy for $F\neq\mathbb{F}_2$; results 2.1.17–2.1.19 and the §2.3
+   "Comments" treat, and explicitly except, the $\mathbb{F}_2$ case.)*
+3. J. Dieudonné. *Sur les générateurs des groupes classiques.* Summa Brasiliensis Mathematicae,
+   vol. 3, pp. 149–179, 1955. *(Original proof of the transvection-length theorem for
+   $F\neq\mathbb{F}_2$, as cited by O'Meara §2.3.)*
+4. D. Callan. *The generation of $\mathrm{Sp}(\mathbb{F}_2)$ by transvections.* Journal of Algebra,
+   vol. 42, no. 2, pp. 378–390, 1976. *(Section 1.1 proves the $r/r+1$ bound,
+   §2.4 identifies the class-A residue-$3$ exceptions, and Theorem 5.1 gives
+   the complete binary exception list.)*
+5. U. Spengler and H. Wolff. *Die Länge einer symplektischen Abbildung* ("The length of a symplectic
+   map"). Journal für die reine und angewandte Mathematik, vol. 274/275, pp. 150–157, 1975. *(The
+   transvection-length function itself, as cited by O'Meara §2.3.)*
+6. J. D. Botha. *Triangularizing matrices over GF(2) by congruence.* Linear and Multilinear Algebra,
+   vol. 42, no. 2, pp. 109–158, 1997.
+   DOI [10.1080/03081089708818495](https://doi.org/10.1080/03081089708818495).
+   *(GF(2) congruence triangularization.)*
+7. D. Maslov, M. Roetteler. *Shorter stabilizer circuits via Bruhat decomposition and quantum
+   circuit transformations.* IEEE Transactions on Information Theory, vol. 64, no. 7, pp. 4729–4738,
+   2018. *(Related symplectic/Bruhat structure.)*
diff --git a/paulimer/src/clifford.rs b/paulimer/src/clifford.rs
index 5fa75826..73a76b51 100644
--- a/paulimer/src/clifford.rs
+++ b/paulimer/src/clifford.rs
@@ -308,5 +308,8 @@ pub use clifford_impl::{
     z_images_partition_transform,
 };
 
+mod transvection;
+pub use transvection::{clifford_centralizer, clifford_to_transvections, clifford_to_transvections_minimal};
+
 #[derive(Debug, PartialEq, Eq, Default)]
 pub struct CliffordStringParsingError;
diff --git a/paulimer/src/clifford/transvection.rs b/paulimer/src/clifford/transvection.rs
new file mode 100644
index 00000000..52ef0ec0
--- /dev/null
+++ b/paulimer/src/clifford/transvection.rs
@@ -0,0 +1,671 @@
+//! Minimal decomposition of Clifford unitaries into Clifford transvections (`π/4` Pauli exponents).
+//!
+//! A *Clifford transvection* is the `π/4` Pauli exponent `exp(iπ/4·P_v)`, whose conjugation action
+//! on Pauli operators is the *symplectic transvection*
+//!
+//! ```text
+//! x ↦ x + ⟨x, v⟩ v,
+//! ```
+//!
+//! where `⟨·,·⟩` is the symplectic (commutation) form. This module follows the transvection
+//! framework of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380) (Pllaha, Volanto & Tirkkonen,
+//! *Decomposition of Clifford Gates*). The exact minimum is `r` or `r + 1`, where
+//! `r = 2n − dim Fix(F)`; congruence-triangularizability of the residue core, not hyperbolicity
+//! alone, decides which value occurs.
+//!
+//! Two decompositions are provided:
+//!
+//! * [`clifford_to_transvections`] uses a greedy O'Meara-style reduction: it always produces a
+//!   **linear number of factors** (`O(n)`), reproducing the symplectic action exactly, but it is
+//!   **not guaranteed to hit the strict `r`/`r + 1` minimum** — intermediate maps can become
+//!   hyperbolic, adding an occasional extra factor.
+//! * [`clifford_to_transvections_minimal`] produces the **strict minimum** number of factors
+//!   (`r` or `r + 1`) via a congruence-triangulation of the residue core.
+//!
+//! Unlike a sign-exact decomposition into Pauli exponents (which reproduces the full signed tableau,
+//! and hence an exact global phase when replayed on a phased operator, at `O(n²)` factors via
+//! Gaussian elimination), these decompositions reproduce only the **symplectic action** — they
+//! ignore Pauli-image signs and the global phase. Their advantage is the linear factor count `O(n)`.
+//!
+//! ## The minimum factor count
+//!
+//! [arXiv:2102.11380](https://arxiv.org/abs/2102.11380) states before and within Theorem 3 that the
+//! residue matrix `F̂` of any *non-hyperbolic* symplectic map can be triangularized by congruence,
+//! giving a decomposition into exactly `r = dim Res(F)` transvections. This is **not correct**:
+//! there exist non-hyperbolic maps whose residue core is *not* congruence-triangularizable and
+//! which therefore require `r + 1` transvections. The smallest examples occur on two qubits:
+//! `T_X₀ T_X₁ T_{X₀X₁} T_Z₀` has residue rank `3` and minimal length `4` despite being
+//! non-hyperbolic. The correct criterion, used here, is: the minimum is `r` when the residue core is
+//! congruence-triangularizable and `r + 1` otherwise (hyperbolicity is the special case where the
+//! core is *alternating*).
+
+use std::collections::HashSet;
+
+use binar::matrix::{AlignedBitMatrix, kernel_basis_matrix};
+use binar::{Bitwise, IndexSet};
+
+use crate::clifford::{Clifford, CliffordMutable, CliffordUnitary};
+use crate::pauli::DensePauli;
+use crate::{Pauli, PauliBinaryOps, PauliMutable, SparsePauli, anti_commutes_with};
+
+/// Decomposes `clifford` into an ordered product of Clifford transvections.
+///
+/// Returns a list of Pauli operators `[P₁, …, P_k]` such that left-multiplying the identity by the
+/// transvections `exp(iπ/4·P₁)`, then `exp(iπ/4·P₂)`, …, then `exp(iπ/4·P_k)` reproduces the
+/// **symplectic action** of `clifford` (its conjugation map on Pauli operators). The Pauli-image
+/// signs and the global phase are *not* reproduced; a sign-exact decomposition into Pauli
+/// exponents would preserve them, at the cost of `O(n²)` factors (see the module docs).
+///
+/// The number of factors is **linear** in the qubit count (`O(n)`). The strict minimum is either
+/// `r` or `r + 1`, where `r = 2n − dim Fix(clifford)`; the greedy reduction here can add an
+/// occasional extra factor when an intermediate map becomes hyperbolic. The count is always at
+/// least `r`.
+///
+/// Every factor is returned with phase exponent `0`; the sign of a transvection does not affect its
+/// symplectic action, so `exp(iπ/4·P)` and `exp(−iπ/4·P)` are interchangeable here.
+///
+/// The exact congruence search has exponential worst-case running time and memoization space in the
+/// residue rank. Candidates are generated lazily rather than materializing the full residue-space
+/// span up front. For large Cliffords where strict minimality is unnecessary, prefer
+/// [`clifford_to_transvections`].
+///
+/// # Examples
+///
+/// ```
+/// use paulimer::CliffordUnitary;
+/// use paulimer::clifford::{clifford_to_transvections, Clifford, CliffordMutable};
+///
+/// let mut clifford = CliffordUnitary::identity(2);
+/// clifford.left_mul_hadamard(0);
+/// clifford.left_mul_cx(0, 1);
+///
+/// let transvections = clifford_to_transvections(&clifford);
+///
+/// let mut rebuilt = CliffordUnitary::identity(2);
+/// for pauli in &transvections {
+///     rebuilt.left_mul_pauli_exp(pauli);
+/// }
+/// // The symplectic actions agree (signs and global phase may differ).
+/// assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+/// ```
+#[must_use]
+pub fn clifford_to_transvections(clifford: &CliffordUnitary) -> Vec<SparsePauli> {
+    let qubit_count = clifford.num_qubits();
+    let mut working = clifford.clone();
+    let mut recorded = Vec::new();
+    // Reduce the symplectic action to the identity by left-multiplying transvections `T_{v₁}, …,
+    // T_{v_k}`, so that `T_{v_k} ⋯ T_{v₁} · F = I` and hence `F = T_{v₁} ⋯ T_{v_k}`. Replaying the
+    // factors in reverse order rebuilds `F` from the identity.
+    while let Some(transvection) = next_transvection(&working) {
+        working.left_mul_pauli_exp(&transvection);
+        recorded.push(transvection);
+        debug_assert!(
+            recorded.len() <= 4 * qubit_count + 2,
+            "transvection reduction exceeded its linear termination bound"
+        );
+    }
+    recorded.reverse();
+    recorded
+}
+
+/// Returns generators of the Clifford's centralizer: the Pauli operators fixed (up to sign) by
+/// conjugation, i.e. the `P` with `clifford · P · clifford† = ±P`.
+///
+/// This is `Fix(F)`, the kernel of the residue map `P ↦ conj(P) · P`, computed as the left null
+/// space of the residue matrix over GF(2). The returned Paulis are independent generators (with
+/// phase exponent `0`); the centralizer they span has dimension `dim Fix(F) = 2n − r`, where `r` is
+/// the residue rank and a lower bound on every transvection decomposition.
+///
+/// # Examples
+///
+/// ```
+/// use paulimer::{CliffordUnitary, Pauli};
+/// use paulimer::clifford::{clifford_centralizer, Clifford, CliffordMutable};
+///
+/// let mut clifford = CliffordUnitary::identity(1);
+/// clifford.left_mul_root_z(0); // S fixes Z, sends X -> Y
+///
+/// let generators = clifford_centralizer(&clifford);
+/// // Every generator is fixed (up to sign) under conjugation.
+/// assert!(generators.iter().all(|pauli| {
+///     let image = clifford.image(pauli);
+///     image.x_bits() == pauli.x_bits() && image.z_bits() == pauli.z_bits()
+/// }));
+/// ```
+#[must_use]
+pub fn clifford_centralizer(clifford: &CliffordUnitary) -> Vec<SparsePauli> {
+    let qubit_count = clifford.num_qubits();
+    let dimension = 2 * qubit_count;
+    let mut residue = AlignedBitMatrix::zeros(dimension, dimension);
+    for (row, basis) in symplectic_basis(qubit_count).enumerate() {
+        let vector = residue_vector(&basis, &clifford.image(&basis));
+        for qubit in 0..qubit_count {
+            if vector.x_bits().index(qubit) {
+                residue.set((row, qubit), true);
+            }
+            if vector.z_bits().index(qubit) {
+                residue.set((row, qubit_count + qubit), true);
+            }
+        }
+    }
+    let kernel = kernel_basis_matrix(&residue.transposed());
+    (0..kernel.row_count())
+        .map(|row| {
+            let x_bits: IndexSet = (0..qubit_count).filter(|&qubit| kernel[(row, qubit)]).collect();
+            let z_bits: IndexSet = (0..qubit_count)
+                .filter(|&qubit| kernel[(row, qubit_count + qubit)])
+                .collect();
+            SparsePauli::from_bits(x_bits, z_bits, 0)
+        })
+        .collect()
+}
+
+/// The `2n` standard basis Pauli operators `X₀, …, X_{n−1}, Z₀, …, Z_{n−1}`.
+fn symplectic_basis(qubit_count: usize) -> impl Iterator<Item = SparsePauli> {
+    (0..qubit_count)
+        .map(move |qubit| SparsePauli::x(qubit, qubit_count))
+        .chain((0..qubit_count).map(move |qubit| SparsePauli::z(qubit, qubit_count)))
+}
+
+/// The next transvection `T_v` reducing the residue of `working`, or `None` if `working` already
+/// acts as the identity on Pauli operators (up to sign).
+///
+/// Following the O'Meara strategy of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380): find a
+/// vector `x` with `⟨x, conj(x)⟩ = 1` (`x` anticommutes with its own image) and set `v = x + conj(x)`
+/// — a residue vector — which lowers the residue rank by one. If no such `x` exists but `working`
+/// is non-trivial (the hyperbolic case), any nonzero residue vector `v` makes the action
+/// non-hyperbolic while preserving the residue space, costing one extra transvection.
+fn next_transvection(working: &CliffordUnitary) -> Option<SparsePauli> {
+    let qubit_count = working.num_qubits();
+    let basis: Vec<SparsePauli> = symplectic_basis(qubit_count).collect();
+    let images: Vec<DensePauli> = basis.iter().map(|pauli| working.image(pauli)).collect();
+
+    for (pauli, image) in basis.iter().zip(&images) {
+        if anti_commutes_with(pauli, image) {
+            return Some(residue_vector(pauli, image));
+        }
+    }
+
+    let dimension = basis.len();
+    for first in 0..dimension {
+        for second in (first + 1)..dimension {
+            let anticommuting =
+                anti_commutes_with(&basis[first], &images[second]) ^ anti_commutes_with(&basis[second], &images[first]);
+            if anticommuting {
+                let mut sum = basis[first].clone();
+                sum.mul_assign_left(&basis[second]);
+                let mut image = images[first].clone();
+                image.mul_assign_left(&images[second]);
+                return Some(residue_vector(&sum, &image));
+            }
+        }
+    }
+
+    basis
+        .iter()
+        .zip(&images)
+        .find(|(pauli, image)| !acts_trivially_on(pauli, image))
+        .map(|(pauli, image)| residue_vector(pauli, image))
+}
+
+/// The residue vector `v = x + conj(x)` as a phaseless Pauli (its symplectic vector is the product
+/// `x · conj(x)`).
+fn residue_vector(pauli: &SparsePauli, image: &DensePauli) -> SparsePauli {
+    let mut vector: SparsePauli = image.clone().into();
+    vector.mul_assign_left(pauli);
+    vector.assign_phase_exp(0);
+    vector
+}
+
+/// Whether `image` equals `pauli` as a symplectic vector (i.e. conjugation fixes `pauli` up to sign).
+fn acts_trivially_on(pauli: &SparsePauli, image: &DensePauli) -> bool {
+    let mut difference: SparsePauli = image.clone().into();
+    difference.mul_assign_left(pauli);
+    difference.x_bits().is_zero() && difference.z_bits().is_zero()
+}
+
+/// Decomposes `clifford` into a **minimal** ordered product of Clifford transvections.
+///
+/// Returns a list of Pauli operators `[P₁, …, P_k]` such that left-multiplying the identity by the
+/// transvections `exp(iπ/4·P₁)`, then `exp(iπ/4·P₂)`, …, then `exp(iπ/4·P_k)` reproduces the
+/// **symplectic action** of `clifford` (its conjugation map on Pauli operators). The Pauli-image
+/// signs and the global phase are *not* reproduced; a sign-exact decomposition into Pauli
+/// exponents would preserve them, at the cost of `O(n²)` factors (see the module docs).
+///
+/// The number of factors `k` is the strict minimum: `k = r` when the residue core is
+/// congruence-triangularizable and `k = r + 1` otherwise, where `r = 2n − dim Fix(clifford)` is the
+/// dimension of the residue space (see [`clifford_centralizer`] for `Fix`). This corrects the
+/// minimality criterion of [arXiv:2102.11380](https://arxiv.org/abs/2102.11380) (see the module
+/// docs). Contrast with [`clifford_to_transvections`], which is only near-minimal.
+///
+/// Every factor is returned with phase exponent `0`; the sign of a transvection does not affect its
+/// symplectic action, so `exp(iπ/4·P)` and `exp(−iπ/4·P)` are interchangeable here.
+///
+/// # Examples
+///
+/// ```
+/// use paulimer::CliffordUnitary;
+/// use paulimer::clifford::{clifford_to_transvections_minimal, Clifford, CliffordMutable};
+///
+/// let mut clifford = CliffordUnitary::identity(2);
+/// clifford.left_mul_hadamard(0);
+/// clifford.left_mul_cx(0, 1);
+///
+/// let transvections = clifford_to_transvections_minimal(&clifford);
+///
+/// let mut rebuilt = CliffordUnitary::identity(2);
+/// for pauli in &transvections {
+///     rebuilt.left_mul_pauli_exp(pauli);
+/// }
+/// // The symplectic actions agree (signs and global phase may differ).
+/// assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+/// ```
+#[must_use]
+pub fn clifford_to_transvections_minimal(clifford: &CliffordUnitary) -> Vec<SparsePauli> {
+    let qubit_count = clifford.num_qubits();
+    let action = action_matrix(clifford);
+    let vectors = minimal_decomposition(&action, qubit_count);
+    vectors
+        .iter()
+        .map(|vector| vector_to_pauli(vector, qubit_count))
+        .collect()
+}
+
+/// The `2n × 2n` symplectic action matrix of `clifford`, in the "image" convention: row `k` is the
+/// symplectic vector of the image of the `k`-th standard basis Pauli (`X₀, …, X_{n−1}, Z₀, …,
+/// Z_{n−1}`), with `x`-bits in columns `[0, n)` and `z`-bits in columns `[n, 2n)`.
+fn action_matrix(clifford: &CliffordUnitary) -> AlignedBitMatrix {
+    let qubit_count = clifford.num_qubits();
+    let dimension = 2 * qubit_count;
+    let mut matrix = AlignedBitMatrix::zeros(dimension, dimension);
+    for (row, basis) in symplectic_basis(qubit_count).enumerate() {
+        let image = clifford.image(&basis);
+        for qubit in 0..qubit_count {
+            if image.x_bits().index(qubit) {
+                matrix.set((row, qubit), true);
+            }
+            if image.z_bits().index(qubit) {
+                matrix.set((row, qubit_count + qubit), true);
+            }
+        }
+    }
+    matrix
+}
+
+/// The symplectic transvection matrix `T_v` (row `k` = `e_k + ⟨e_k, v⟩·v`), whose row-vector action
+/// `x ↦ x·T_v` equals `x + ⟨x, v⟩·v`.
+fn transvection_matrix(vector: &[bool], qubit_count: usize) -> AlignedBitMatrix {
+    let dimension = 2 * qubit_count;
+    let mut matrix = AlignedBitMatrix::identity(dimension);
+    for row in 0..dimension {
+        let coupling = if row < qubit_count {
+            vector[qubit_count + row]
+        } else {
+            vector[row - qubit_count]
+        };
+        if coupling {
+            for (column, &bit) in vector.iter().enumerate() {
+                if bit {
+                    matrix.negate((row, column));
+                }
+            }
+        }
+    }
+    matrix
+}
+
+/// The residue matrix `F̂ = Ω·(I + F)`, where `Ω` swaps the `x` and `z` halves of the rows. Its row
+/// space is the residue space `Res(F)`.
+fn residue_matrix(action: &AlignedBitMatrix, qubit_count: usize) -> AlignedBitMatrix {
+    let dimension = 2 * qubit_count;
+    let mut residue = AlignedBitMatrix::zeros(dimension, dimension);
+    for row in 0..dimension {
+        let swapped = if row < qubit_count {
+            row + qubit_count
+        } else {
+            row - qubit_count
+        };
+        for column in 0..dimension {
+            let mut bit = action.get((swapped, column));
+            if swapped == column {
+                bit ^= true;
+            }
+            if bit {
+                residue.set((row, column), true);
+            }
+        }
+    }
+    residue
+}
+
+/// Row-reduces `matrix` to reduced echelon form while tracking the transform.
+///
+/// Returns `(basis, transform)` where `basis` holds the `r` nonzero echelon rows (a basis of the row
+/// space) and `transform` is `r × rows` with `basis = transform · matrix`. Pivoting is over the
+/// columns of `matrix` only.
+fn row_reduce_with_transform(matrix: &AlignedBitMatrix) -> (AlignedBitMatrix, AlignedBitMatrix) {
+    let rows = matrix.row_count();
+    let columns = matrix.column_count();
+    let mut augmented = AlignedBitMatrix::zeros(rows, columns + rows);
+    for row in 0..rows {
+        for column in 0..columns {
+            if matrix.get((row, column)) {
+                augmented.set((row, column), true);
+            }
+        }
+        augmented.set((row, columns + row), true);
+    }
+    let mut pivot_row = 0;
+    for column in 0..columns {
+        let Some(selected) = (pivot_row..rows).find(|&row| augmented.get((row, column))) else {
+            continue;
+        };
+        augmented.swap_rows(pivot_row, selected);
+        for row in 0..rows {
+            if row != pivot_row && augmented.get((row, column)) {
+                augmented.add_into_row(row, pivot_row);
+            }
+        }
+        pivot_row += 1;
+    }
+    let rank = pivot_row;
+    let mut basis = AlignedBitMatrix::zeros(rank, columns);
+    let mut transform = AlignedBitMatrix::zeros(rank, rows);
+    for row in 0..rank {
+        for column in 0..columns {
+            if augmented.get((row, column)) {
+                basis.set((row, column), true);
+            }
+        }
+        for column in 0..rows {
+            if augmented.get((row, columns + column)) {
+                transform.set((row, column), true);
+            }
+        }
+    }
+    (basis, transform)
+}
+
+/// Extracts row `index` of `matrix` as a boolean vector of length `length`.
+fn matrix_row(matrix: &AlignedBitMatrix, index: usize, length: usize) -> Vec<bool> {
+    (0..length).map(|column| matrix.get((index, column))).collect()
+}
+
+/// The bitwise XOR of two equal-length boolean vectors.
+fn xor_vectors(left: &[bool], right: &[bool]) -> Vec<bool> {
+    left.iter().zip(right).map(|(&a, &b)| a ^ b).collect()
+}
+
+/// The value `x·E·yᵀ` of the bilinear form given by the square matrix `core`.
+fn bilinear(core: &AlignedBitMatrix, left: &[bool], right: &[bool]) -> bool {
+    let dimension = core.row_count();
+    (0..dimension).fold(false, |acc, i| {
+        let row = (0..dimension).fold(false, |inner, j| inner ^ (core.get((i, j)) & right[j]));
+        acc ^ (left[i] & row)
+    })
+}
+
+/// Packs `vectors` (each of length `columns`) into an `AlignedBitMatrix`.
+fn vectors_to_matrix(vectors: &[Vec<bool>], columns: usize) -> AlignedBitMatrix {
+    let mut matrix = AlignedBitMatrix::zeros(vectors.len(), columns);
+    for (row, vector) in vectors.iter().enumerate() {
+        for (column, &bit) in vector.iter().enumerate() {
+            if bit {
+                matrix.set((row, column), true);
+            }
+        }
+    }
+    matrix
+}
+
+/// Attempts to triangularize the `r × r` matrix `core` by congruence.
+///
+/// On success returns `Ok(q)` with `q ∈ GL(r, 2)` such that `q·core·qᵀ` is lower triangular; the
+/// rows of `q` are an ordered basis in which each vector is right-orthogonal (under the form
+/// `x·core·yᵀ`) to all later ones and non-isotropic (`x·core·xᵀ = 1`). Since `core` is invertible,
+/// a lower-triangular `q·core·qᵀ` automatically has an all-ones diagonal.
+///
+/// A triangularization exists exactly when the associated symplectic map is a product of `r`
+/// transvections. It is found by a backtracking search over the choice of each successive basis
+/// vector: after picking a non-isotropic `pick`, the search recurses into its right-orthogonal
+/// complement. A greedy (first-choice) search can dead-end even when a triangularization exists, so
+/// the choices are explored exhaustively, with subspaces proven unsolvable memoized to prune the
+/// search. On failure returns `Err(())`.
+fn congruence_triangularize(core: &AlignedBitMatrix) -> Result<AlignedBitMatrix, ()> {
+    let dimension = core.row_count();
+    if dimension == 0 {
+        return Ok(AlignedBitMatrix::zeros(0, 0));
+    }
+    let standard: Vec<Vec<bool>> = (0..dimension)
+        .map(|index| (0..dimension).map(|column| column == index).collect())
+        .collect();
+    let mut unsolvable: HashSet<Vec<bool>> = HashSet::new();
+    triangularize_subspace(core, &standard, dimension, &mut unsolvable)
+        .map(|picks| vectors_to_matrix(&picks, dimension))
+        .ok_or(())
+}
+
+/// Backtracking core of [`congruence_triangularize`]: finds an ordered basis of `span(basis)` in
+/// which each vector is non-isotropic and right-orthogonal to all later ones, or `None` if none
+/// exists. Subspaces proven to have no such basis are recorded in `unsolvable` (keyed by their
+/// canonical row-reduced form) so that they are never re-explored.
+fn triangularize_subspace(
+    core: &AlignedBitMatrix,
+    basis: &[Vec<bool>],
+    dimension: usize,
+    unsolvable: &mut HashSet<Vec<bool>>,
+) -> Option<Vec<Vec<bool>>> {
+    if basis.is_empty() {
+        return Some(Vec::new());
+    }
+    let key = subspace_key(basis, dimension);
+    if unsolvable.contains(&key) {
+        return None;
+    }
+    let mut explored: HashSet<Vec<bool>> = HashSet::new();
+    for pick in span_vectors(basis) {
+        if !bilinear(core, &pick, &pick) {
+            continue;
+        }
+        let Some(complement) = right_orthogonal_complement(core, &pick, basis) else {
+            continue;
+        };
+        let complement_key = subspace_key(&complement, dimension);
+        if !explored.insert(complement_key) {
+            continue;
+        }
+        if let Some(mut rest) = triangularize_subspace(core, &complement, dimension, unsolvable) {
+            let mut picks = Vec::with_capacity(rest.len() + 1);
+            picks.push(pick);
+            picks.append(&mut rest);
+            return Some(picks);
+        }
+    }
+    unsolvable.insert(key);
+    None
+}
+
+/// Lazily generates all `2ᵈ − 1` nonzero vectors in the span of a `d`-vector basis.
+struct SpanVectors<'a> {
+    basis: &'a [Vec<bool>],
+    coefficients: Vec<bool>,
+    current: Vec<bool>,
+    exhausted: bool,
+}
+
+impl<'a> SpanVectors<'a> {
+    fn new(basis: &'a [Vec<bool>]) -> Self {
+        Self {
+            basis,
+            coefficients: vec![false; basis.len()],
+            current: vec![false; basis.first().map_or(0, Vec::len)],
+            exhausted: basis.is_empty(),
+        }
+    }
+}
+
+impl Iterator for SpanVectors<'_> {
+    type Item = Vec<bool>;
+
+    fn next(&mut self) -> Option<Self::Item> {
+        if self.exhausted {
+            return None;
+        }
+        for (index, member) in self.basis.iter().enumerate() {
+            self.coefficients[index] ^= true;
+            for (slot, &bit) in self.current.iter_mut().zip(member) {
+                *slot ^= bit;
+            }
+            if self.coefficients[index] {
+                return Some(self.current.clone());
+            }
+        }
+        self.exhausted = true;
+        None
+    }
+}
+
+fn span_vectors(basis: &[Vec<bool>]) -> SpanVectors<'_> {
+    SpanVectors::new(basis)
+}
+
+/// A basis of `{y ∈ span(basis) : pick·core·yᵀ = 0}`, one dimension smaller than `basis`, or `None`
+/// if `pick` is right-orthogonal to the whole span (which cannot happen for a non-isotropic `pick`).
+fn right_orthogonal_complement(core: &AlignedBitMatrix, pick: &[bool], basis: &[Vec<bool>]) -> Option<Vec<Vec<bool>>> {
+    let couplings: Vec<bool> = basis.iter().map(|vector| bilinear(core, pick, vector)).collect();
+    let pivot = couplings.iter().position(|&bit| bit)?;
+    let mut complement = Vec::with_capacity(basis.len() - 1);
+    for (index, vector) in basis.iter().enumerate() {
+        if index == pivot {
+            continue;
+        }
+        if couplings[index] {
+            complement.push(xor_vectors(vector, &basis[pivot]));
+        } else {
+            complement.push(vector.clone());
+        }
+    }
+    Some(complement)
+}
+
+/// A canonical key for the subspace spanned by `basis`: its rows reduced to reduced row-echelon
+/// form and flattened, so that any two bases of the same subspace produce the same key.
+fn subspace_key(basis: &[Vec<bool>], dimension: usize) -> Vec<bool> {
+    let mut rows: Vec<Vec<bool>> = basis.to_vec();
+    let mut pivot = 0;
+    for column in 0..dimension {
+        let Some(selected) = (pivot..rows.len()).find(|&row| rows[row][column]) else {
+            continue;
+        };
+        rows.swap(pivot, selected);
+        for row in 0..rows.len() {
+            if row != pivot && rows[row][column] {
+                let reference = rows[pivot].clone();
+                for (slot, bit) in rows[row].iter_mut().zip(&reference) {
+                    *slot ^= *bit;
+                }
+            }
+        }
+        pivot += 1;
+    }
+    rows.truncate(pivot);
+    rows.into_iter().flatten().collect()
+}
+
+/// The residue core `E` and its residue-space basis `V` for the action matrix `action`.
+///
+/// Returns `(basis, rank, core)` where `basis` (`rank × 2n`) spans `Res(F)` and `core = V·Rᵀ` with
+/// `V = R·F̂` (`rank × rank`) is the matrix whose congruence-triangularizability governs minimality.
+fn residue_core(action: &AlignedBitMatrix, qubit_count: usize) -> (AlignedBitMatrix, usize, AlignedBitMatrix) {
+    let residue = residue_matrix(action, qubit_count);
+    let (basis, transform) = row_reduce_with_transform(&residue);
+    let rank = basis.row_count();
+    let core = basis.dot(&transform.transposed());
+    (basis, rank, core)
+}
+
+/// The minimal ordered transvection vectors for the symplectic action matrix `action`.
+fn minimal_decomposition(action: &AlignedBitMatrix, qubit_count: usize) -> Vec<Vec<bool>> {
+    let dimension = 2 * qubit_count;
+    let (basis, rank, core) = residue_core(action, qubit_count);
+    if rank == 0 {
+        return Vec::new();
+    }
+    let Ok(transform) = congruence_triangularize(&core) else {
+        let fix = find_fix_vector(action, qubit_count, &basis, rank);
+        let updated = action.dot(&transvection_matrix(&fix, qubit_count));
+        let mut vectors = minimal_decomposition(&updated, qubit_count);
+        vectors.push(fix);
+        return vectors;
+    };
+    let defining = transform.dot(&basis);
+    (0..rank).map(|row| matrix_row(&defining, row, dimension)).collect()
+}
+
+/// Finds a residue vector `v` such that `F·T_v` has a congruence-triangularizable residue core of
+/// the same rank, so that `F` decomposes into `rank + 1` transvections. Such a vector always exists
+/// in `Res(F)` (the map is a product of `rank + 1` transvections, and dropping the last factor
+/// leaves a product of `rank` transvections whose residue core is triangularizable).
+///
+/// Candidates are the nonzero residue vectors in ascending binary-coordinate order. The search is
+/// exhaustive over `Res(F)` and therefore always succeeds.
+fn find_fix_vector(action: &AlignedBitMatrix, qubit_count: usize, basis: &AlignedBitMatrix, rank: usize) -> Vec<bool> {
+    let dimension = 2 * qubit_count;
+    let lift = |coordinates: &[bool]| -> Vec<bool> {
+        let mut vector = vec![false; dimension];
+        for (row, &selected) in coordinates.iter().enumerate() {
+            if selected {
+                for (column, slot) in vector.iter_mut().enumerate() {
+                    *slot ^= basis.get((row, column));
+                }
+            }
+        }
+        vector
+    };
+    let candidate_accepts = |vector: &[bool]| -> bool {
+        if vector.iter().all(|&bit| !bit) {
+            return false;
+        }
+        let updated = action.dot(&transvection_matrix(vector, qubit_count));
+        let (_, updated_rank, updated_core) = residue_core(&updated, qubit_count);
+        updated_rank == rank && congruence_triangularize(&updated_core).is_ok()
+    };
+    let coordinate_basis: Vec<Vec<bool>> = (0..rank)
+        .map(|selected| (0..rank).map(|index| index == selected).collect())
+        .collect();
+    for coordinates in span_vectors(&coordinate_basis) {
+        let vector = lift(&coordinates);
+        if candidate_accepts(&vector) {
+            return vector;
+        }
+    }
+    unreachable!("a residue fix vector always exists for a non-triangularizable core")
+}
+
+/// Converts a `2n`-bit symplectic vector into a phaseless Pauli (`x`-bits in `[0, n)`, `z`-bits in
+/// `[n, 2n)`).
+fn vector_to_pauli(vector: &[bool], qubit_count: usize) -> SparsePauli {
+    let x_bits: IndexSet = (0..qubit_count).filter(|&qubit| vector[qubit]).collect();
+    let z_bits: IndexSet = (0..qubit_count).filter(|&qubit| vector[qubit_count + qubit]).collect();
+    SparsePauli::from_bits(x_bits, z_bits, 0)
+}
+
+#[cfg(test)]
+mod tests {
+    use super::span_vectors;
+
+    #[test]
+    fn span_vectors_supports_more_than_u64_bits_lazily() {
+        let dimension = 65;
+        let basis: Vec<Vec<bool>> = (0..dimension)
+            .map(|row| (0..dimension).map(|column| row == column).collect())
+            .collect();
+        let vectors: Vec<Vec<bool>> = span_vectors(&basis).take(4).collect();
+
+        assert_eq!(vectors.len(), 4);
+        assert_eq!(vectors[0].iter().filter(|&&bit| bit).count(), 1);
+        assert_eq!(vectors[1].iter().filter(|&&bit| bit).count(), 1);
+        assert_eq!(vectors[2].iter().filter(|&&bit| bit).count(), 2);
+        assert_eq!(vectors[3].iter().filter(|&&bit| bit).count(), 1);
+    }
+}
diff --git a/paulimer/tests/transvection_test.rs b/paulimer/tests/transvection_test.rs
new file mode 100644
index 00000000..c7299737
--- /dev/null
+++ b/paulimer/tests/transvection_test.rs
@@ -0,0 +1,540 @@
+//! Tests for the Clifford -> transvection decomposition (arXiv:2102.11380).
+//!
+//! The greedy decomposition reproduces the *symplectic action* (ignoring Pauli-image signs and the
+//! global phase) with a linear number of factors. It is not guaranteed to hit the strict minimum,
+//! so its tests validate the symplectic-action round trip, the linear factor bound, and the
+//! centralizer contract. Separate tests cover the minimal decomposition.
+
+use binar::{Bitwise, IndexSet};
+use paulimer::UnitaryOp;
+use paulimer::clifford::{Clifford, CliffordMutable, CliffordUnitary, clifford_centralizer, clifford_to_transvections};
+use paulimer::pauli::{Pauli, SparsePauli};
+use proptest::collection::vec;
+use proptest::prelude::*;
+use rand::SeedableRng;
+use rand::rngs::StdRng;
+
+/// Rebuilds a Clifford's symplectic action by replaying transvections on the identity.
+fn symplectic_action_from_transvections(transvections: &[SparsePauli], qubit_count: usize) -> CliffordUnitary {
+    let mut rebuilt = CliffordUnitary::identity(qubit_count);
+    for transvection in transvections {
+        rebuilt.left_mul_pauli_exp(transvection);
+    }
+    rebuilt
+}
+
+/// Whether conjugation by `clifford` fixes `pauli` as a symplectic vector (ignoring sign).
+fn is_conjugation_fixed(clifford: &CliffordUnitary, pauli: &SparsePauli) -> bool {
+    let image = clifford.image(pauli);
+    image.x_bits() == pauli.x_bits() && image.z_bits() == pauli.z_bits()
+}
+
+fn is_non_identity(pauli: &SparsePauli) -> bool {
+    !(pauli.x_bits().is_zero() && pauli.z_bits().is_zero())
+}
+
+/// The residue rank `r = 2n - dim Fix(F)`, a lower bound on every decomposition.
+fn residue_rank(clifford: &CliffordUnitary) -> usize {
+    2 * clifford.num_qubits() - clifford_centralizer(clifford).len()
+}
+
+fn assert_valid_decomposition(clifford: &CliffordUnitary) {
+    let qubit_count = clifford.num_qubits();
+    let transvections = clifford_to_transvections(clifford);
+
+    let rebuilt = symplectic_action_from_transvections(&transvections, qubit_count);
+    assert_eq!(
+        rebuilt.symplectic_matrix(),
+        clifford.symplectic_matrix(),
+        "replayed transvections must reproduce the symplectic action"
+    );
+
+    for transvection in &transvections {
+        assert_eq!(transvection.xz_phase_exponent(), 0, "factors carry no phase");
+        assert!(is_non_identity(transvection), "factors are non-identity Paulis");
+    }
+
+    let lower_bound = residue_rank(clifford);
+    assert!(
+        transvections.len() >= lower_bound,
+        "a decomposition cannot be shorter than the residue rank {lower_bound}, got {}",
+        transvections.len()
+    );
+    assert!(
+        transvections.len() <= 4 * qubit_count + 2,
+        "the decomposition must be linear in the qubit count, got {}",
+        transvections.len()
+    );
+}
+
+#[test]
+fn identity_decomposes_to_no_transvections() {
+    for qubit_count in 0..5 {
+        let identity = CliffordUnitary::identity(qubit_count);
+        let transvections = clifford_to_transvections(&identity);
+        assert!(
+            transvections.is_empty(),
+            "identity has no transvections (qubit_count {qubit_count})"
+        );
+        let centralizer = clifford_centralizer(&identity);
+        assert_eq!(
+            centralizer.len(),
+            2 * qubit_count,
+            "identity commutes with all {qubit_count} Pauli generators"
+        );
+    }
+}
+
+#[test]
+fn single_qubit_gates_reproduce_symplectic_action() {
+    let mut s_gate = CliffordUnitary::identity(1);
+    s_gate.left_mul_root_z(0);
+    assert_valid_decomposition(&s_gate);
+    assert_eq!(clifford_to_transvections(&s_gate).len(), 1, "S is one transvection T_Z");
+
+    let mut hadamard = CliffordUnitary::identity(1);
+    hadamard.left_mul_hadamard(0);
+    assert_valid_decomposition(&hadamard);
+    assert_eq!(
+        clifford_to_transvections(&hadamard).len(),
+        1,
+        "H is the transvection T_Y"
+    );
+}
+
+#[test]
+fn pauli_gates_are_conjugation_trivial() {
+    // Pauli operators act trivially by conjugation (sign-only), so their symplectic action is the
+    // identity and no transvections are needed.
+    for axis in 0..3 {
+        let mut clifford = CliffordUnitary::identity(1);
+        match axis {
+            0 => clifford.left_mul_pauli(&SparsePauli::x(0, 1)),
+            1 => clifford.left_mul_pauli(&SparsePauli::z(0, 1)),
+            _ => clifford.left_mul_pauli(&SparsePauli::y(0, 1)),
+        }
+        assert!(
+            clifford_to_transvections(&clifford).is_empty(),
+            "Pauli axis {axis} needs no factor"
+        );
+        assert_eq!(
+            clifford_centralizer(&clifford).len(),
+            2,
+            "a Pauli commutes with all generators"
+        );
+    }
+}
+
+#[test]
+fn swap_exercises_the_hyperbolic_branch() {
+    // SWAP is hyperbolic (its residue space is totally isotropic), so the greedy reduction returns
+    // r + 1 = 3 transvections, where r = 2n - dim Fix = 4 - 2 = 2.
+    let mut swap = CliffordUnitary::identity(2);
+    swap.left_mul_swap(0, 1);
+    assert_valid_decomposition(&swap);
+    assert_eq!(residue_rank(&swap), 2);
+    assert_eq!(clifford_to_transvections(&swap).len(), 3);
+    assert_eq!(clifford_centralizer(&swap).len(), 2);
+}
+
+#[test]
+fn two_qubit_gates_reproduce_symplectic_action() {
+    let mut cx = CliffordUnitary::identity(2);
+    cx.left_mul_cx(0, 1);
+    assert_valid_decomposition(&cx);
+
+    let mut cz = CliffordUnitary::identity(2);
+    cz.left_mul_cz(0, 1);
+    assert_valid_decomposition(&cz);
+}
+
+#[test]
+fn composite_circuit_reproduces_symplectic_action() {
+    let mut clifford = CliffordUnitary::identity(4);
+    clifford.left_mul_hadamard(0);
+    clifford.left_mul_cx(0, 1);
+    clifford.left_mul_root_z(2);
+    clifford.left_mul_cz(1, 3);
+    clifford.left_mul_swap(2, 3);
+    clifford.left_mul_hadamard(3);
+    assert_valid_decomposition(&clifford);
+}
+
+#[test]
+fn centralizer_generators_are_conjugation_fixed_and_independent() {
+    let mut clifford = CliffordUnitary::identity(3);
+    clifford.left_mul_hadamard(0);
+    clifford.left_mul_cx(0, 1);
+    clifford.left_mul_root_z(2);
+
+    let centralizer = clifford_centralizer(&clifford);
+    assert!(centralizer.iter().all(|pauli| is_conjugation_fixed(&clifford, pauli)));
+    assert!(centralizer.iter().all(is_non_identity));
+    assert_eq!(
+        centralizer.len(),
+        2 * clifford.num_qubits() - residue_rank(&clifford),
+        "the centralizer dimension is 2n - r"
+    );
+}
+
+fn random_clifford(qubit_count: usize, seed: u64) -> CliffordUnitary {
+    let mut random_number_generator = StdRng::seed_from_u64(seed);
+    CliffordUnitary::random(qubit_count, &mut random_number_generator)
+}
+
+#[test]
+fn many_random_cliffords_reproduce_symplectic_action() {
+    // A deterministic sweep giving broad coverage independent of the proptest shrink budget.
+    for qubit_count in 0..7 {
+        for seed in 0..200 {
+            assert_valid_decomposition(&random_clifford(qubit_count, seed));
+        }
+    }
+}
+
+/// A single Clifford generator, modeled as an operation so proptest can shrink a failing input down
+/// to a minimal gate sequence (unlike an opaque RNG seed).
+#[derive(Clone, Debug)]
+enum Gate {
+    Single { op: UnitaryOp, qubit: usize },
+    Two { op: UnitaryOp, first: usize, second: usize },
+}
+
+fn distinct_pair(qubit_count: usize) -> impl Strategy<Value = (usize, usize)> {
+    (0..qubit_count, 0..qubit_count - 1)
+        .prop_map(|(first, second)| (first, if second < first { second } else { second + 1 }))
+}
+
+fn gate_strategy(qubit_count: usize) -> BoxedStrategy<Gate> {
+    use UnitaryOp::{ControlledX, ControlledZ, Hadamard, SqrtX, SqrtZ, Swap, X, Y, Z};
+    let single = (
+        prop::sample::select(vec![Hadamard, SqrtZ, SqrtX, X, Y, Z]),
+        0..qubit_count,
+    )
+        .prop_map(|(op, qubit)| Gate::Single { op, qubit });
+    if qubit_count < 2 {
+        return single.boxed();
+    }
+    let two = (
+        prop::sample::select(vec![ControlledX, ControlledZ, Swap]),
+        distinct_pair(qubit_count),
+    )
+        .prop_map(|(op, (first, second))| Gate::Two { op, first, second });
+    prop_oneof![3 => single, 1 => two].boxed()
+}
+
+fn clifford_from_gates(qubit_count: usize, gates: &[Gate]) -> CliffordUnitary {
+    let mut clifford = CliffordUnitary::identity(qubit_count);
+    for gate in gates {
+        match *gate {
+            Gate::Single { op, qubit } => clifford.left_mul(op, &[qubit]),
+            Gate::Two { op, first, second } => clifford.left_mul(op, &[first, second]),
+        }
+    }
+    clifford
+}
+
+/// A qubit count paired with a random gate sequence acting on it.
+fn scenario() -> impl Strategy<Value = (usize, Vec<Gate>)> {
+    (1usize..7).prop_flat_map(|qubit_count| {
+        vec(gate_strategy(qubit_count), 0..=3 * qubit_count).prop_map(move |gates| (qubit_count, gates))
+    })
+}
+
+proptest! {
+    #![proptest_config(ProptestConfig::with_cases(512))]
+
+    #[test]
+    fn reproduces_symplectic_action((qubit_count, gates) in scenario()) {
+        let clifford = clifford_from_gates(qubit_count, &gates);
+        let transvections = clifford_to_transvections(&clifford);
+        let rebuilt = symplectic_action_from_transvections(&transvections, qubit_count);
+        prop_assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+    }
+
+    #[test]
+    fn decomposition_is_linear_and_no_shorter_than_minimum((qubit_count, gates) in scenario()) {
+        let clifford = clifford_from_gates(qubit_count, &gates);
+        let transvection_count = clifford_to_transvections(&clifford).len();
+        let lower_bound = residue_rank(&clifford);
+        prop_assert!(
+            transvection_count >= lower_bound,
+            "got {transvection_count} factors, below the residue rank {lower_bound}"
+        );
+        prop_assert!(
+            transvection_count <= 4 * qubit_count + 2,
+            "got {transvection_count} factors, above the linear bound"
+        );
+    }
+
+    #[test]
+    fn centralizer_is_conjugation_fixed((qubit_count, gates) in scenario()) {
+        let clifford = clifford_from_gates(qubit_count, &gates);
+        for generator in clifford_centralizer(&clifford) {
+            prop_assert!(is_conjugation_fixed(&clifford, &generator));
+            prop_assert!(is_non_identity(&generator), "centralizer generators must be non-identity");
+        }
+    }
+}
+
+use paulimer::clifford::clifford_to_transvections_minimal;
+use std::collections::HashMap;
+
+/// A symplectic action matrix over GF(2) as a row-major boolean grid (test-local, used only by the
+/// brute-force minimality oracle).
+type ActionMatrix = Vec<Vec<bool>>;
+
+/// The image-convention symplectic action of `clifford`: row `k` is the image of the `k`-th standard
+/// basis Pauli. The minimal transvection length is a conjugation invariant, so any faithful matrix
+/// realization yields the same brute-force minimum.
+fn action_of(clifford: &CliffordUnitary) -> ActionMatrix {
+    let qubit_count = clifford.num_qubits();
+    let dimension = 2 * qubit_count;
+    let basis: Vec<SparsePauli> = (0..qubit_count)
+        .map(|qubit| SparsePauli::x(qubit, qubit_count))
+        .chain((0..qubit_count).map(|qubit| SparsePauli::z(qubit, qubit_count)))
+        .collect();
+    let mut matrix = vec![vec![false; dimension]; dimension];
+    for (row, pauli) in basis.iter().enumerate() {
+        let image = clifford.image(pauli);
+        for qubit in 0..qubit_count {
+            matrix[row][qubit] = image.x_bits().index(qubit);
+            matrix[row][qubit_count + qubit] = image.z_bits().index(qubit);
+        }
+    }
+    matrix
+}
+
+fn multiply(left: &ActionMatrix, right: &ActionMatrix) -> ActionMatrix {
+    let dimension = left.len();
+    let mut product = vec![vec![false; dimension]; dimension];
+    for i in 0..dimension {
+        for k in 0..dimension {
+            if left[i][k] {
+                for j in 0..dimension {
+                    product[i][j] ^= right[k][j];
+                }
+            }
+        }
+    }
+    product
+}
+
+fn transvection(vector: &[bool], qubit_count: usize) -> ActionMatrix {
+    let dimension = 2 * qubit_count;
+    let mut matrix = vec![vec![false; dimension]; dimension];
+    for (row, output) in matrix.iter_mut().enumerate() {
+        output[row] = true;
+        let coupling = if row < qubit_count {
+            vector[qubit_count + row]
+        } else {
+            vector[row - qubit_count]
+        };
+        if coupling {
+            for (column, slot) in output.iter_mut().enumerate() {
+                *slot ^= vector[column];
+            }
+        }
+    }
+    matrix
+}
+
+fn encode(matrix: &ActionMatrix) -> u32 {
+    let mut key = 0u32;
+    let mut bit = 0;
+    for row in matrix {
+        for &value in row {
+            if value {
+                key |= 1 << bit;
+            }
+            bit += 1;
+        }
+    }
+    key
+}
+
+/// The exact minimal transvection length of `clifford`'s symplectic action, by breadth-first search
+/// over the symplectic group. Only tractable for small qubit counts (`n <= 2`).
+fn minimal_length_oracle(clifford: &CliffordUnitary) -> usize {
+    let qubit_count = clifford.num_qubits();
+    let dimension = 2 * qubit_count;
+    let identity: ActionMatrix = (0..dimension)
+        .map(|i| (0..dimension).map(|j| i == j).collect())
+        .collect();
+    let target = encode(&action_of(clifford));
+    let generators: Vec<ActionMatrix> = (1..(1u32 << dimension))
+        .map(|mask| {
+            let vector: Vec<bool> = (0..dimension).map(|bit| mask & (1 << bit) != 0).collect();
+            transvection(&vector, qubit_count)
+        })
+        .collect();
+    let mut distances: HashMap<u32, usize> = HashMap::new();
+    distances.insert(encode(&identity), 0);
+    let mut frontier = vec![identity];
+    let mut distance = 0;
+    while !frontier.is_empty() {
+        if distances.contains_key(&target) {
+            break;
+        }
+        let mut next = Vec::new();
+        for current in &frontier {
+            for generator in &generators {
+                let product = multiply(current, generator);
+                let key = encode(&product);
+                if let std::collections::hash_map::Entry::Vacant(entry) = distances.entry(key) {
+                    entry.insert(distance + 1);
+                    next.push(product);
+                }
+            }
+        }
+        frontier = next;
+        distance += 1;
+    }
+    distances[&target]
+}
+
+fn assert_valid_minimal_decomposition(clifford: &CliffordUnitary) {
+    let qubit_count = clifford.num_qubits();
+    let transvections = clifford_to_transvections_minimal(clifford);
+
+    let rebuilt = symplectic_action_from_transvections(&transvections, qubit_count);
+    assert_eq!(
+        rebuilt.symplectic_matrix(),
+        clifford.symplectic_matrix(),
+        "replayed transvections must reproduce the symplectic action"
+    );
+
+    for transvection in &transvections {
+        assert_eq!(transvection.xz_phase_exponent(), 0, "factors carry no phase");
+        assert!(is_non_identity(transvection), "factors are non-identity Paulis");
+    }
+
+    let rank = residue_rank(clifford);
+    assert!(
+        transvections.len() == rank || transvections.len() == rank + 1,
+        "the minimal count is r or r + 1 (r = {rank}), got {}",
+        transvections.len()
+    );
+    assert!(
+        transvections.len() <= clifford_to_transvections(clifford).len(),
+        "the minimal decomposition cannot exceed the greedy one"
+    );
+}
+
+#[test]
+fn minimal_identity_decomposes_to_no_transvections() {
+    for qubit_count in 0..5 {
+        assert!(clifford_to_transvections_minimal(&CliffordUnitary::identity(qubit_count)).is_empty());
+    }
+}
+
+#[test]
+fn minimal_single_qubit_gates() {
+    let mut s_gate = CliffordUnitary::identity(1);
+    s_gate.left_mul_root_z(0);
+    assert_valid_minimal_decomposition(&s_gate);
+    assert_eq!(clifford_to_transvections_minimal(&s_gate).len(), 1);
+
+    let mut hadamard = CliffordUnitary::identity(1);
+    hadamard.left_mul_hadamard(0);
+    assert_valid_minimal_decomposition(&hadamard);
+    assert_eq!(clifford_to_transvections_minimal(&hadamard).len(), 1);
+}
+
+#[test]
+fn minimal_swap_needs_r_plus_one() {
+    let mut swap = CliffordUnitary::identity(2);
+    swap.left_mul_swap(0, 1);
+    assert_valid_minimal_decomposition(&swap);
+    assert_eq!(residue_rank(&swap), 2);
+    assert_eq!(clifford_to_transvections_minimal(&swap).len(), 3);
+}
+
+#[test]
+fn minimal_callan_class_a_needs_r_plus_one() {
+    let centers = [
+        SparsePauli::x(0, 2),
+        SparsePauli::x(1, 2),
+        SparsePauli::from_bits([0, 1].into_iter().collect(), IndexSet::new(), 0),
+        SparsePauli::z(0, 2),
+    ];
+    let clifford = symplectic_action_from_transvections(&centers, 2);
+    let action = action_of(&clifford);
+
+    assert_eq!(
+        action,
+        vec![
+            vec![true, false, true, false],
+            vec![false, true, false, false],
+            vec![false, true, true, false],
+            vec![true, false, true, true],
+        ]
+    );
+    assert!(action[0][2], "⟨X₀, X₀F⟩ = 1, so F is non-hyperbolic");
+    assert_eq!(residue_rank(&clifford), 3);
+    assert_eq!(minimal_length_oracle(&clifford), 4);
+    assert_eq!(clifford_to_transvections_minimal(&clifford).len(), 4);
+}
+
+#[test]
+fn minimal_two_qubit_gates() {
+    let mut cx = CliffordUnitary::identity(2);
+    cx.left_mul_cx(0, 1);
+    assert_valid_minimal_decomposition(&cx);
+
+    let mut cz = CliffordUnitary::identity(2);
+    cz.left_mul_cz(0, 1);
+    assert_valid_minimal_decomposition(&cz);
+}
+
+#[test]
+fn minimal_composite_circuit() {
+    let mut clifford = CliffordUnitary::identity(4);
+    clifford.left_mul_hadamard(0);
+    clifford.left_mul_cx(0, 1);
+    clifford.left_mul_root_z(2);
+    clifford.left_mul_cz(1, 3);
+    clifford.left_mul_swap(2, 3);
+    clifford.left_mul_hadamard(3);
+    assert_valid_minimal_decomposition(&clifford);
+}
+
+#[test]
+fn minimal_matches_brute_force_oracle_on_one_and_two_qubits() {
+    // Exact minimality against an independent breadth-first search over the symplectic group.
+    for qubit_count in 0..=2 {
+        for seed in 0..400 {
+            let clifford = random_clifford(qubit_count, seed);
+            let decomposed = clifford_to_transvections_minimal(&clifford);
+            let rebuilt = symplectic_action_from_transvections(&decomposed, qubit_count);
+            assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+            assert_eq!(
+                decomposed.len(),
+                minimal_length_oracle(&clifford),
+                "decomposition length must equal the brute-force minimum (n={qubit_count}, seed={seed})"
+            );
+        }
+    }
+}
+
+proptest! {
+    #![proptest_config(ProptestConfig::with_cases(512))]
+
+    #[test]
+    fn minimal_reproduces_symplectic_action((qubit_count, gates) in scenario()) {
+        let clifford = clifford_from_gates(qubit_count, &gates);
+        let transvections = clifford_to_transvections_minimal(&clifford);
+        let rebuilt = symplectic_action_from_transvections(&transvections, qubit_count);
+        prop_assert_eq!(rebuilt.symplectic_matrix(), clifford.symplectic_matrix());
+    }
+
+    #[test]
+    fn minimal_is_r_or_r_plus_one_and_at_most_greedy((qubit_count, gates) in scenario()) {
+        let clifford = clifford_from_gates(qubit_count, &gates);
+        let minimal = clifford_to_transvections_minimal(&clifford).len();
+        let greedy = clifford_to_transvections(&clifford).len();
+        let residue = residue_rank(&clifford);
+        prop_assert!(minimal == residue || minimal == residue + 1, "got {minimal}, r = {residue}");
+        prop_assert!(minimal <= greedy, "minimal {minimal} exceeded greedy {greedy}");
+    }
+}