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 f072dac0..c2185921 100644
--- a/paulimer/bindings/python/paulimer.pyi
+++ b/paulimer/bindings/python/paulimer.pyi
@@ -543,6 +543,18 @@ class CliffordUnitary:
         """
         ...
 
+    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."""
         ...
diff --git a/paulimer/bindings/python/src/py_clifford.rs b/paulimer/bindings/python/src/py_clifford.rs
index 4ef6c88e..e055e5b3 100644
--- a/paulimer/bindings/python/src/py_clifford.rs
+++ b/paulimer/bindings/python/src/py_clifford.rs
@@ -1,6 +1,7 @@
 use derive_more::{Deref, DerefMut, From, Into};
 use paulimer::clifford::{
-    clifford_centralizer, clifford_to_transvections, group_encoding_clifford_of, split_phased_css,
+    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};
@@ -280,12 +281,28 @@ impl PyCliffordUnitary {
     ///
     /// 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; see
-    /// :meth:`to_pauli_exponents` for the sign-exact (but ``O(n^2)``) decomposition.
+    /// 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> {
diff --git a/paulimer/bindings/python/tests/transvection_test.py b/paulimer/bindings/python/tests/transvection_test.py
index a72ece10..829ab9af 100644
--- a/paulimer/bindings/python/tests/transvection_test.py
+++ b/paulimer/bindings/python/tests/transvection_test.py
@@ -41,28 +41,50 @@ def _assert_valid_decomposition(clifford):
     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
 
 
@@ -70,6 +92,7 @@ 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():
@@ -112,6 +135,12 @@ 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):
diff --git a/paulimer/docs/transvection-minimality-correction.md b/paulimer/docs/transvection-minimality-correction.md
new file mode 100644
index 00000000..819b2ba0
--- /dev/null
+++ b/paulimer/docs/transvection-minimality-correction.md
@@ -0,0 +1,328 @@
+# A correction to the minimal transvection decomposition of Clifford gates (arXiv:2102.11380)
+
+This note documents two issues we uncovered 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).
+
+while implementing [`clifford_to_transvections_minimal`](../src/clifford/transvection.rs)
+(the minimal-length decomposition of a Clifford's symplectic action into `π/4` Pauli exponents).
+It states the correct minimality criterion we adopted, and backs every claim with a finite,
+machine-checkable computation. 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 (Theorem 6, "Transvection Decomposition of Symplectic Matrices", and the
+paragraph preceding it) 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,
+the paper argues (in the paragraph immediately following Lemmas 2–3; line 411 of the arXiv v1
+source) that the *residue matrix* $\widehat{\mathbf F}$ **can always be
+triangularized by congruence when $\mathbf F$ is non-hyperbolic**, which is what would make the
+length-$r$ decomposition exist.
+
+**This is not correct over $\mathbb{F}_2$** — the field of interest for qubit Cliffords. The clean
+"$r$ if non-hyperbolic, $r+1$ if hyperbolic" dichotomy is a classical theorem of Dieudonné (see
+O'Meara, *Symplectic Groups*, Theorem 2.1.11), but that theorem is **stated only for fields
+$F\neq\mathbb{F}_2$**, and O'Meara explicitly warns that over $\mathbb{F}_2$ "the theorem fails … it
+is no longer possible to express every $\sigma$ … as a product of $\operatorname{res}\sigma$ or of
+$(\operatorname{res}\sigma)+1$ transvections" (§2.3, Comments). The paper invokes the dichotomy over
+exactly the one field the classical result excludes. Consequently there exist non-hyperbolic
+symplectic maps whose residue core is *not* congruence-triangularizable and whose minimal
+transvection length is therefore $r+1$, not $r$. The **smallest such example already occurs on two
+qubits** ($m=2$), with residue rank $r=3$ and minimal length $4$; we exhibit one explicitly and
+verify it two independent ways by exhaustive search.
+
+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.
+
+(The exact binary length function in full generality is the more intricate object studied by Callan
+and by Spengler–Wolff; see §6 and the references. Our criterion and the $r/(r+1)$ range are
+verified computationally for up to six qubits.)
+
+Triangularizability is *strictly stronger* than being non-alternating; it depends on the
+non-symmetric part of $\mathbf E$ and is not captured by any invariant of the associated quadratic
+form alone. This is exactly the subtle question studied in Botha's work on GF(2) congruence
+triangularization, which the paper itself cites but does not use to qualify the claim.
+
+## 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,
+$$
+
+and $\mathbf x\,\widehat{\mathbf F}^{\mathsf T}\mathbf x^{\mathsf T}=\langle\mathbf x,\mathbf x\mathbf F\rangle$,
+so $\widehat{\mathbf F}$ has all-zero diagonal iff $\mathbf F$ is hyperbolic (in which case
+$\widehat{\mathbf F}$ is *alternating*: symmetric 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.
+
+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 6 (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$.
+
+**Why the attribution does not carry over to $\mathbb{F}_2$.** The cited length statement is
+O'Meara's Theorem 2.1.11 (originally Dieudonné): *if $\sigma\neq 1$ is non-hyperbolic it is a product
+of $\operatorname{res}\sigma$ transvections; if hyperbolic, of $(\operatorname{res}\sigma)+1$ but not
+$\operatorname{res}\sigma$.* **Its hypothesis is $F\neq\mathbb{F}_2$.** O'Meara devotes separate
+results (2.1.17–2.1.19) to characteristic $2$ and proves the dichotomy for $\mathbb{F}_2$ only for
+*involutions* (2.1.18) and, in the general case, again only for $F\neq\mathbb{F}_2$ (2.1.19). His
+§2.3 "Comments" then states plainly: "If the underlying field is $\mathbb{F}_2$, then the theorem
+fails … There is a theorem for $\mathbb{F}_2$, but it is considerably more complicated," pointing to
+Callan (1976) and to Spengler–Wolff, *Die Länge einer symplektischen Abbildung* ("The length of a
+symplectic map"). The proof of 2.1.19 makes the gap explicit: it needs, at each step, a transvection
+$\mathbf T$ with $\operatorname{res}(\mathbf T\sigma)<\operatorname{res}\sigma$ **and $\mathbf T\sigma$
+still non-hyperbolic**, and finding such a $\mathbf T$ uses the extra field elements available only
+when $F\neq\mathbb{F}_2$.
+
+**The flaw in the paper's own justification.** Independently of the mis-cited hypothesis, 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. Issue 1: non-hyperbolic does **not** imply triangularizable (a machine-checked counterexample)
+
+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 v,\mathbf v\mathbf F\rangle=1$ for some $\mathbf v$, 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 6's assumption that a triangularizing $\mathbf Q$ always exists — are false.
+
+**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
+```
+
+**Non-hyperbolic maps of this kind are the rule, not the exception.** 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. Issue 2: even when a triangularization exists, the greedy search is incomplete
+
+The paper does not give an explicit triangularization procedure; it defers to Botha's algorithms.
+A natural but naive implementation extends O'Meara's idea directly: pick a vector $\mathbf u$ with
+$\psi_{\mathbf E}(\mathbf u)=1$ (a "unit-diagonal pivot"), use it as the first triangularization
+step, and recurse into its right-orthogonal complement. **This forward-greedy search is incomplete
+for $r\ge 5$:** a locally valid pivot choice can drive the remaining subspace to become entirely
+$\psi$-isotropic (alternating) and dead-end, *even when a triangularization of the whole core
+exists* via a different sequence of pivots. We encountered this on a 4-qubit instance (residue rank
+$r=7$) where the greedy triangularization returns an obstruction although the core is in fact
+triangularizable and a length-$r$ decomposition exists; a one-step look-ahead pivot rule does not
+fix it either. This case is preserved as a regression seed in
+[`tests/transvection_test.proptest-regressions`](../tests/transvection_test.proptest-regressions).
+
+This is a *practical* pitfall distinct from Issue 1: Issue 1 is a false mathematical claim (the
+target $\mathbf Q$ may not exist); Issue 2 is that *finding* $\mathbf Q$ when it does exist requires
+more than a greedy pivot walk.
+
+## 6. The correction we adopted
+
+Combining the (correct) Lemmas 2–3 with the two issues above, the length our algorithm produces is:
+
+$$
+\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. Memoization keeps the search tractable (a few thousand nodes even
+  at $r=9$) while guaranteeing completeness — fixing Issue 2.
+- **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
+  algorithm omits — and fixes Issue 1.
+
+We verify the result in the test suite against a brute-force BFS oracle on one and two qubits, and
+against the $\{r,\,r+1\}$ range on up to six qubits. Two honest caveats about *strict minimality* for
+large systems: (i) O'Meara's §2.3 warns that over $\mathbb{F}_2$ the exact length function is
+"considerably more complicated" than $r/(r+1)$, so we do **not** claim $\ell(\mathbf F)\in\{r,r+1\}$
+holds for *all* $m$ — only that it does in every case we have checked (through $m=6$), and that
+whenever `find_fix_vector` succeeds the returned length $r+1$ *is* minimal (since a non-triangularizable
+core rules out length $r$). The exact binary length is the object studied by Callan and by
+Spengler–Wolff. (ii) Our `find_fix_vector` restores triangularizability with a *single* extra
+transvection in all tested cases; a rigorous proof that one fix always suffices — or a construction
+handling the rare cases where it might not for large $m$ — is left as a follow-up.
+
+**Why "non-alternating" is not enough.** Triangularizability of $\mathbf E$ is not determined by any
+symmetric invariant of $\psi_{\mathbf E}$: neither the Arf invariant of $\psi_{\mathbf E}$ nor
+whether $\psi_{\mathbf E}$ is nonzero on the radical of its polar form
+$B(\mathbf x,\mathbf y)=\mathbf x(\mathbf E+\mathbf E^{\mathsf T})\mathbf y^{\mathsf T}$ decides it —
+solvability depends on the *non-symmetric* part of $\mathbf E$. (Empirically, an odd Arf invariant or
+a $\psi$ that is nonzero on the radical is always solvable; only the remaining regime is mixed,
+which is why no closed-form symmetric criterion exists and a search is required.) The precise
+characterization of GF(2) congruence triangularization is the subject of Botha (1997), which the
+paper cites but does not use to qualify Theorem 6.
+
+## 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.
+
+## 8. References (verified)
+
+Bibliographic details are taken from the corrected paper's reference list and from O'Meara's own
+bibliography and §2.3, and confirmed against the 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. *(The $\mathbb{F}_2$ case, which O'Meara notes is "considerably
+   more complicated" and which Callan shows Dieudonné's treatment handled incompletely.)*
+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). *(The precise
+   criterion and algorithms for GF(2) congruence triangularization — the subject the flawed claim
+   glosses over.)*
+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 bdc586ea..73a76b51 100644
--- a/paulimer/src/clifford.rs
+++ b/paulimer/src/clifford.rs
@@ -309,7 +309,7 @@ pub use clifford_impl::{
 };
 
 mod transvection;
-pub use transvection::{clifford_centralizer, clifford_to_transvections};
+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
index 60b00890..68cf5379 100644
--- a/paulimer/src/clifford/transvection.rs
+++ b/paulimer/src/clifford/transvection.rs
@@ -13,18 +13,33 @@
 //! *minimal* number of factors is `r = 2n − dim Fix(F)` (or `r + 1` when the symplectic action `F`
 //! is hyperbolic), where `Fix(F)` is the space of Pauli operators fixed by conjugation.
 //!
-//! The decomposition here 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. In practice it stays within a small additive constant of the
-//! minimum. The strict-minimum variant (via the paper's congruence-triangulation machinery) is
-//! tracked as a follow-up.
+//! Two decompositions are provided:
 //!
-//! Unlike [`clifford_to_pauli_exponents`](super::clifford_to_pauli_exponents), which reproduces the
-//! full signed tableau (and hence an exact global phase when replayed on a phased operator), this
-//! decomposition reproduces only the **symplectic action** — it ignores Pauli-image signs and the
-//! global phase. Its advantage is the linear factor count `O(n)`, versus `O(n²)` for the
-//! Gaussian-elimination decomposition.
+//! * [`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 matrix.
+//!
+//! 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 (Theorem 1) 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 already on two qubits; for instance
+//! the symplectic action with residue rank `3` fixed by the standard basis order
+//! `X₀, X₁, Z₀, Z₁` requires four transvections 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::{kernel_basis_matrix, AlignedBitMatrix};
 use binar::{Bitwise, IndexSet};
@@ -38,8 +53,8 @@ use crate::{anti_commutes_with, Pauli, PauliBinaryOps, PauliMutable, SparsePauli
 /// 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; see the module docs for the contrast with
-/// [`clifford_to_pauli_exponents`](super::clifford_to_pauli_exponents).
+/// 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)`). It is close to, but not
 /// guaranteed to equal, the strict minimum `r = 2n − dim Fix(clifford)` (`r + 1` when the
@@ -50,6 +65,11 @@ use crate::{anti_commutes_with, Pauli, PauliBinaryOps, PauliMutable, SparsePauli
 /// 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
 ///
 /// ```
@@ -202,3 +222,458 @@ fn acts_trivially_on(pauli: &SparsePauli, image: &DensePauli) -> bool {
     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.proptest-regressions b/paulimer/tests/transvection_test.proptest-regressions
new file mode 100644
index 00000000..2c160b84
--- /dev/null
+++ b/paulimer/tests/transvection_test.proptest-regressions
@@ -0,0 +1,7 @@
+# Seeds for failure cases proptest has generated in the past. It is
+# automatically read and these particular cases re-run before any
+# novel cases are generated.
+#
+# It is recommended to check this file in to source control so that
+# everyone who runs the test benefits from these saved cases.
+cc 7fdf3d7e2bef80dfea588254563fa2aa313b4e4b882bffeb71506d4f33678e5d # shrinks to qubit_count = 4, seed = 6167596993315164505
diff --git a/paulimer/tests/transvection_test.rs b/paulimer/tests/transvection_test.rs
index 0e38dba7..031dc432 100644
--- a/paulimer/tests/transvection_test.rs
+++ b/paulimer/tests/transvection_test.rs
@@ -220,3 +220,234 @@ proptest! {
         }
     }
 }
+
+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 minimum = residue_rank(clifford);
+    assert!(
+        transvections.len() == minimum || transvections.len() == minimum + 1,
+        "the minimal count is r or r + 1 (r = {minimum}), 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_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 in 0usize..7, seed in any::<u64>()) {
+        let clifford = random_clifford(qubit_count, seed);
+        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 in 0usize..7, seed in any::<u64>()) {
+        let clifford = random_clifford(qubit_count, seed);
+        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}");
+    }
+}