From 53ebf1b81d71acd0a649bb184cb9ddbca1d8fe95 Mon Sep 17 00:00:00 2001
From: Hood Chatham <roberthoodchatham@gmail.com>
Date: Thu, 9 Jul 2026 11:44:47 -0700
Subject: [PATCH 1/4] Add ext_py matrices and vectors

Adds a basic Python project scafold including basic vector bindings and tests,
and CI integration.
---
 .github/workflows/ext_py.yaml            |   43 +
 ext/crates/fp/src/field/dyn_field.rs     |  139 +++
 ext/crates/fp/src/field/mod.rs           |    2 +
 ext/crates/fp/src/matrix/matrix_inner.rs |    4 +
 ext_py/.gitignore                        |    2 +
 ext_py/Cargo.toml                        |   19 +
 ext_py/pyproject.toml                    |   24 +
 ext_py/python/ext/__init__.py            |    1 +
 ext_py/src/fp_mod.rs                     |  494 +++++++++
 ext_py/src/fp_mod/matrices.rs            | 1296 ++++++++++++++++++++++
 ext_py/src/fp_mod/vectors.rs             |  701 ++++++++++++
 ext_py/src/lib.rs                        |   12 +
 ext_py/tests/test_affine_subspace.py     |  181 +++
 ext_py/tests/test_augmented_matrix.py    |  249 +++++
 ext_py/tests/test_fp_vector_slices.py    |  329 ++++++
 ext_py/tests/test_matrix.py              |  386 +++++++
 ext_py/tests/test_matrix_slice_mut.py    |  196 ++++
 ext_py/tests/test_quasi_inverse.py       |  251 +++++
 ext_py/tests/test_subquotient.py         |  228 ++++
 ext_py/tests/test_subspace.py            |  173 +++
 20 files changed, 4730 insertions(+)
 create mode 100644 .github/workflows/ext_py.yaml
 create mode 100644 ext/crates/fp/src/field/dyn_field.rs
 create mode 100644 ext_py/.gitignore
 create mode 100644 ext_py/Cargo.toml
 create mode 100644 ext_py/pyproject.toml
 create mode 100644 ext_py/python/ext/__init__.py
 create mode 100644 ext_py/src/fp_mod.rs
 create mode 100644 ext_py/src/fp_mod/matrices.rs
 create mode 100644 ext_py/src/fp_mod/vectors.rs
 create mode 100644 ext_py/src/lib.rs
 create mode 100644 ext_py/tests/test_affine_subspace.py
 create mode 100644 ext_py/tests/test_augmented_matrix.py
 create mode 100644 ext_py/tests/test_fp_vector_slices.py
 create mode 100644 ext_py/tests/test_matrix.py
 create mode 100644 ext_py/tests/test_matrix_slice_mut.py
 create mode 100644 ext_py/tests/test_quasi_inverse.py
 create mode 100644 ext_py/tests/test_subquotient.py
 create mode 100644 ext_py/tests/test_subspace.py

diff --git a/.github/workflows/ext_py.yaml b/.github/workflows/ext_py.yaml
new file mode 100644
index 0000000000..11ef35878e
--- /dev/null
+++ b/.github/workflows/ext_py.yaml
@@ -0,0 +1,43 @@
+name: Test ext_py
+
+on: [push, pull_request]
+
+jobs:
+  test:
+    if: ${{ github.event_name == 'push' || github.event.pull_request.head.repo.full_name != github.repository }}
+    runs-on: ubuntu-latest
+    env:
+      RUST_BACKTRACE: 1
+      RUSTFLAGS: "-D warnings"
+
+    steps:
+      - uses: actions/checkout@v5
+
+      - name: Set up Python 3.14
+        uses: actions/setup-python@v6
+        with:
+          python-version: "3.14"
+
+      - name: Install rustup
+        id: rustup
+        uses: dtolnay/rust-toolchain@v1
+        with:
+          toolchain: stable
+
+      - name: Cache files
+        uses: actions/cache@v5
+        with:
+          path: |
+            ~/.cargo
+            /usr/share/rust/.cargo
+            **/target
+          key: ext-py-${{ steps.rustup.outputs.rustc_hash }}-${{ hashFiles('**/Cargo.*') }}
+
+      - name: Build wheel
+        run: cd ext_py && pip install --group dev && maturin build --out dist
+
+      - name: Install wheel and run tests
+        run: |
+          cd ext_py
+          pip install --force-reinstall dist/*.whl
+          pytest tests
diff --git a/ext/crates/fp/src/field/dyn_field.rs b/ext/crates/fp/src/field/dyn_field.rs
new file mode 100644
index 0000000000..0f8e5f1c0c
--- /dev/null
+++ b/ext/crates/fp/src/field/dyn_field.rs
@@ -0,0 +1,139 @@
+use std::{
+    fmt,
+    ops::{Add, Mul, Neg, Sub},
+};
+
+use super::{Fp, SmallFq, element::FieldElement};
+use crate::prime::ValidPrime;
+
+/// A field element whose field type has been erased to one of the two concrete kinds the public
+/// API exposes: a prime field [`Fp`] or a small extension field [`SmallFq`], each over a dynamic
+/// [`ValidPrime`].
+///
+/// This lets callers that cannot be generic over the [`Field`](super::Field) trait — most notably
+/// the Python bindings, where a `pyclass` cannot carry generic parameters — manipulate field
+/// elements of either kind through a single type, with the arithmetic operators implemented once
+/// here rather than re-matched at every call site.
+///
+/// Combining two elements that do not live in the same field (a different kind, or the same kind
+/// over a different prime/degree) is a genuine failure, not a value: the binary operators report it
+/// as [`None`]/[`Err`] rather than silently coercing one operand or returning a default element.
+#[derive(Clone, Copy, PartialEq, Eq, Hash, Debug)]
+pub enum DynFieldElement {
+    Fp(FieldElement<Fp<ValidPrime>>),
+    SmallFq(FieldElement<SmallFq<ValidPrime>>),
+}
+
+/// Why [`DynFieldElement::try_div`] could not divide.
+///
+/// Division has two distinct failure modes that callers may want to surface differently (the Python
+/// bindings map them to `ValueError` and `ZeroDivisionError` respectively), so unlike the additive
+/// operators this is reported with a dedicated enum rather than a bare [`Option`].
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum DivError {
+    /// The two operands do not live in the same field.
+    MismatchedField,
+    /// The divisor is the zero element.
+    DivisionByZero,
+}
+
+impl fmt::Display for DivError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Self::MismatchedField => f.write_str("cannot combine elements from different fields"),
+            Self::DivisionByZero => f.write_str("division by zero"),
+        }
+    }
+}
+
+impl std::error::Error for DivError {}
+
+/// Implement an additive/multiplicative operator once, returning `None` when the operands are not
+/// in the same field. The field-equality guard is necessary: [`FieldElement`]'s own operators use
+/// the left operand's field unconditionally, so without it a mismatch would silently compute a
+/// bogus result instead of signalling failure.
+macro_rules! impl_binop {
+    ($trait:ident, $method:ident) => {
+        impl $trait for DynFieldElement {
+            type Output = Option<Self>;
+
+            fn $method(self, rhs: Self) -> Self::Output {
+                match (self, rhs) {
+                    (Self::Fp(a), Self::Fp(b)) if a.field() == b.field() => {
+                        Some(Self::Fp(a.$method(b)))
+                    }
+                    (Self::SmallFq(a), Self::SmallFq(b)) if a.field() == b.field() => {
+                        Some(Self::SmallFq(a.$method(b)))
+                    }
+                    _ => None,
+                }
+            }
+        }
+    };
+}
+
+impl_binop!(Add, add);
+impl_binop!(Sub, sub);
+impl_binop!(Mul, mul);
+
+impl Neg for DynFieldElement {
+    type Output = Self;
+
+    fn neg(self) -> Self::Output {
+        match self {
+            Self::Fp(x) => Self::Fp(-x),
+            Self::SmallFq(x) => Self::SmallFq(-x),
+        }
+    }
+}
+
+impl DynFieldElement {
+    /// Divide, distinguishing a mismatched-field failure from division by zero. Returns
+    /// [`DivError::MismatchedField`] when the operands are not in the same field and
+    /// [`DivError::DivisionByZero`] when `rhs` is the zero element.
+    pub fn try_div(self, rhs: Self) -> Result<Self, DivError> {
+        match (self, rhs) {
+            (Self::Fp(a), Self::Fp(b)) if a.field() == b.field() => {
+                (a / b).map(Self::Fp).ok_or(DivError::DivisionByZero)
+            }
+            (Self::SmallFq(a), Self::SmallFq(b)) if a.field() == b.field() => {
+                (a / b).map(Self::SmallFq).ok_or(DivError::DivisionByZero)
+            }
+            _ => Err(DivError::MismatchedField),
+        }
+    }
+
+    /// The multiplicative inverse, or [`None`] for the zero element.
+    pub fn inv(self) -> Option<Self> {
+        match self {
+            Self::Fp(x) => x.inv().map(Self::Fp),
+            Self::SmallFq(x) => x.inv().map(Self::SmallFq),
+        }
+    }
+
+    /// The Frobenius endomorphism `x -> x^p`.
+    pub fn frobenius(self) -> Self {
+        match self {
+            Self::Fp(x) => Self::Fp(x.frobenius()),
+            Self::SmallFq(x) => Self::SmallFq(x.frobenius()),
+        }
+    }
+
+    /// The canonical `u32` value of a prime-field element, or [`None`] for a [`SmallFq`] element,
+    /// which has no canonical integer representative.
+    pub fn try_as_u32(self) -> Option<u32> {
+        match self {
+            Self::Fp(x) => Some(*x),
+            Self::SmallFq(_) => None,
+        }
+    }
+}
+
+impl fmt::Display for DynFieldElement {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Self::Fp(x) => write!(f, "{x}"),
+            Self::SmallFq(x) => write!(f, "{x}"),
+        }
+    }
+}
diff --git a/ext/crates/fp/src/field/mod.rs b/ext/crates/fp/src/field/mod.rs
index 714a9e735d..1d5d4e08c5 100644
--- a/ext/crates/fp/src/field/mod.rs
+++ b/ext/crates/fp/src/field/mod.rs
@@ -1,12 +1,14 @@
 use self::field_internal::FieldInternal;
 use crate::prime::Prime;
 
+pub mod dyn_field;
 pub mod element;
 pub(crate) mod field_internal;
 
 pub mod fp;
 pub mod smallfq;
 
+pub use dyn_field::{DivError, DynFieldElement};
 use element::FieldElement;
 pub use fp::Fp;
 pub use smallfq::SmallFq;
diff --git a/ext/crates/fp/src/matrix/matrix_inner.rs b/ext/crates/fp/src/matrix/matrix_inner.rs
index 54b90924ad..023894db48 100644
--- a/ext/crates/fp/src/matrix/matrix_inner.rs
+++ b/ext/crates/fp/src/matrix/matrix_inner.rs
@@ -1558,6 +1558,10 @@ impl<'a> MatrixSliceMut<'a> {
             l.add_masked(r, 1, mask);
         }
     }
+
+    pub fn to_vec(&self) -> Vec<Vec<u32>> {
+        self.iter().map(|row| row.iter().collect()).collect()
+    }
 }
 
 #[cfg(test)]
diff --git a/ext_py/.gitignore b/ext_py/.gitignore
new file mode 100644
index 0000000000..5ddc4ed4e8
--- /dev/null
+++ b/ext_py/.gitignore
@@ -0,0 +1,2 @@
+__pycache__
+*.so
diff --git a/ext_py/Cargo.toml b/ext_py/Cargo.toml
new file mode 100644
index 0000000000..5853af808d
--- /dev/null
+++ b/ext_py/Cargo.toml
@@ -0,0 +1,19 @@
+[package]
+name = "ext_py"
+version = "0.1.0"
+edition = "2021"
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+[lib]
+name = "ext_py"
+crate-type = ["cdylib"]
+
+[dependencies]
+anyhow = "1.0"
+fp = { path = "../ext/crates/fp" }
+
+pyo3 = "0.29.0"
+serde_json = "1.0.141"
+
+[dev-dependencies]
+tempfile = "3.20.0"
diff --git a/ext_py/pyproject.toml b/ext_py/pyproject.toml
new file mode 100644
index 0000000000..a4e3a72daa
--- /dev/null
+++ b/ext_py/pyproject.toml
@@ -0,0 +1,24 @@
+[build-system]
+requires = ["maturin>=1.8"]
+build-backend = "maturin"
+
+[project]
+name = "sseq_ext"
+requires-python = ">=3.10"
+classifiers = [
+    "Programming Language :: Rust",
+    "Programming Language :: Python :: Implementation :: CPython",
+    "Programming Language :: Python :: Implementation :: PyPy",
+]
+dynamic = ["version"]
+
+[tool.maturin]
+features = ["pyo3/extension-module"]
+module-name = "ext"
+python-source = "python"
+
+[dependency-groups]
+dev = [
+    "maturin>=1.8",
+    "pytest",
+]
diff --git a/ext_py/python/ext/__init__.py b/ext_py/python/ext/__init__.py
new file mode 100644
index 0000000000..ec66790593
--- /dev/null
+++ b/ext_py/python/ext/__init__.py
@@ -0,0 +1 @@
+from .ext import *
diff --git a/ext_py/src/fp_mod.rs b/ext_py/src/fp_mod.rs
new file mode 100644
index 0000000000..fb62c387f5
--- /dev/null
+++ b/ext_py/src/fp_mod.rs
@@ -0,0 +1,494 @@
+use std::hash::{DefaultHasher, Hash, Hasher};
+
+use fp::{
+    field::{DivError, DynFieldElement, Field, Fp as RustFp, SmallFq as RustSmallFq},
+    matrix::Matrix as RustMatrix,
+    prime::{self, Prime},
+};
+use pyo3::{
+    basic::CompareOp,
+    exceptions::{PyIndexError, PyRuntimeError, PyValueError, PyZeroDivisionError},
+    prelude::*,
+    types::PyBytes,
+    PyResult,
+};
+
+mod matrices;
+mod vectors;
+
+const MAX_VALID_PRIME: u32 = 1 << 31;
+
+type DynFp = RustFp<prime::ValidPrime>;
+type DynSmallFq = RustSmallFq<prime::ValidPrime>;
+
+fn valid_prime(p: u32) -> PyResult<prime::ValidPrime> {
+    if !(2..MAX_VALID_PRIME).contains(&p) {
+        return Err(PyValueError::new_err(format!("{p} is not prime")));
+    }
+    prime::ValidPrime::try_from(p).map_err(|_| PyValueError::new_err(format!("{p} is not prime")))
+}
+
+/// Build a `SmallFq` from a Python-supplied prime and degree.
+fn small_fq(p: u32, degree: u32) -> PyResult<DynSmallFq> {
+    let p = valid_prime(p)?;
+    Ok(DynSmallFq::new(p, degree))
+}
+
+fn py_hash<T: Hash>(value: &T) -> isize {
+    let mut hasher = DefaultHasher::new();
+    value.hash(&mut hasher);
+    match hasher.finish() as isize {
+        -1 => -2,
+        hash => hash,
+    }
+}
+
+fn checked_range(start: usize, end: usize, len: usize) -> PyResult<()> {
+    if start <= end && end <= len {
+        Ok(())
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "range {start}..{end} out of range for vector of length {len}"
+        )))
+    }
+}
+
+pub(crate) fn checked_row(row: usize, rows: usize) -> PyResult<usize> {
+    if row < rows {
+        Ok(row)
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "row {row} out of range for matrix with {rows} rows"
+        )))
+    }
+}
+
+fn borrow_error(err: impl ToString) -> PyErr {
+    PyRuntimeError::new_err(err.to_string())
+}
+
+/// Map any stringifiable error (e.g. `std::io::Error` from
+/// (de)serialization) into the `RuntimeError` used uniformly across the
+/// `to_bytes`/`from_bytes` methods.
+fn io_err(e: impl ToString) -> PyErr {
+    PyRuntimeError::new_err(e.to_string())
+}
+
+/// Run a `to_bytes`-style writer into a fresh buffer and wrap the result as
+/// `PyBytes`, mapping I/O errors through [`io_err`].
+fn serialize_to_pybytes<'py>(
+    py: Python<'py>,
+    f: impl FnOnce(&mut Vec<u8>) -> std::io::Result<()>,
+) -> PyResult<Bound<'py, PyBytes>> {
+    let mut buffer = Vec::new();
+    f(&mut buffer).map_err(io_err)?;
+    Ok(PyBytes::new(py, &buffer))
+}
+
+/// Uniform error for using a value that has been moved out (consumed) by a
+/// consuming method. Mirrors `borrow_error` for the move-and-invalidate
+/// pyclasses (e.g. the augmented matrices).
+fn consumed_error(label: &str) -> PyErr {
+    PyRuntimeError::new_err(format!("{label} has been consumed"))
+}
+
+/// A value that a consuming method can `take()` out, after which any further
+/// access raises `RuntimeError("<label> has been consumed")` instead of
+/// panicking or operating on stale data. Used to model upstream consuming
+/// semantics (`into_*`, `compute_quasi_inverses`) across the PyO3 boundary,
+/// where methods borrow the pyclass and cannot move out of `self` directly.
+pub(crate) struct Consumable<T> {
+    value: Option<T>,
+    label: &'static str,
+}
+
+impl<T> Consumable<T> {
+    pub(crate) fn new(label: &'static str, value: T) -> Self {
+        Self {
+            value: Some(value),
+            label,
+        }
+    }
+
+    pub(crate) fn get(&self) -> PyResult<&T> {
+        self.value
+            .as_ref()
+            .ok_or_else(|| consumed_error(self.label))
+    }
+
+    pub(crate) fn get_mut(&mut self) -> PyResult<&mut T> {
+        self.value
+            .as_mut()
+            .ok_or_else(|| consumed_error(self.label))
+    }
+
+    pub(crate) fn take(&mut self) -> PyResult<T> {
+        self.value.take().ok_or_else(|| consumed_error(self.label))
+    }
+}
+
+fn checked_equal_len(lhs: usize, rhs: usize) -> PyResult<()> {
+    if lhs == rhs {
+        Ok(())
+    } else {
+        Err(PyValueError::new_err(format!(
+            "length mismatch: {lhs} != {rhs}"
+        )))
+    }
+}
+
+fn checked_same_prime(lhs: u32, rhs: u32) -> PyResult<()> {
+    if lhs == rhs {
+        Ok(())
+    } else {
+        Err(PyValueError::new_err(format!(
+            "prime mismatch: {lhs} != {rhs}"
+        )))
+    }
+}
+
+/// Python `repr` of the field an erased element lives in (e.g. `Fp(2)` or
+/// `SmallFq(2, 3)`), matching the `__repr__` of the `Fp`/`SmallFq` classes.
+/// Used to build the mismatched-field error message.
+fn field_repr(x: DynFieldElement) -> String {
+    match x {
+        DynFieldElement::Fp(x) => format!("Fp({})", x.field().characteristic().as_u32()),
+        DynFieldElement::SmallFq(x) => {
+            let f = x.field();
+            format!("SmallFq({}, {})", f.characteristic().as_u32(), f.degree())
+        }
+    }
+}
+
+/// The `ValueError` raised when a binary field operation is given operands
+/// from two different fields (upstream signals this as `None`/`MismatchedField`).
+fn mismatched_field_error(lhs: DynFieldElement, rhs: DynFieldElement) -> PyErr {
+    PyValueError::new_err(format!(
+        "cannot combine elements from {} and {}",
+        field_repr(lhs),
+        field_repr(rhs)
+    ))
+}
+
+/// A matrix-like parent that can back a borrowed row or rectangle view.
+///
+/// A plain `Matrix` is held directly; an `AugmentedMatrix<N>` is held as its
+/// concrete pyclass and accessed through its `Deref<Target = Matrix>` so
+/// that segment rectangles and segment rows can revalidate against the inner
+/// matrix's current dimensions. We keep the parent Python object alive and
+/// reconstruct the underlying Rust matrix view on each call.
+pub(crate) enum MatrixParent {
+    Matrix(Py<matrices::PyMatrix>),
+    Augmented2(Py<matrices::PyAugmentedMatrix2>),
+    Augmented3(Py<matrices::PyAugmentedMatrix3>),
+}
+
+impl MatrixParent {
+    fn clone_ref(&self, py: Python<'_>) -> Self {
+        match self {
+            Self::Matrix(m) => Self::Matrix(m.clone_ref(py)),
+            Self::Augmented2(m) => Self::Augmented2(m.clone_ref(py)),
+            Self::Augmented3(m) => Self::Augmented3(m.clone_ref(py)),
+        }
+    }
+
+    /// Run `f` on the current inner `Matrix`, holding the borrow for the
+    /// duration of the call. Deref coercion turns an `&AugmentedMatrix<N>`
+    /// into the `&Matrix` expected by `f`.
+    fn with_matrix<R>(&self, py: Python<'_>, f: impl FnOnce(&RustMatrix) -> R) -> PyResult<R> {
+        match self {
+            Self::Matrix(m) => Ok(f(&m.try_borrow(py).map_err(borrow_error)?.0)),
+            Self::Augmented2(m) => Ok(f(m.try_borrow(py).map_err(borrow_error)?.0.get()?)),
+            Self::Augmented3(m) => Ok(f(m.try_borrow(py).map_err(borrow_error)?.0.get()?)),
+        }
+    }
+
+    /// Whether `self` and `other` are backed by the same Python object
+    /// (same `Matrix`/`AugmentedMatrix` instance). Different enum variants
+    /// are necessarily different objects. Used to decide whether a shared
+    /// borrow of one and a mutable borrow of the other would collide.
+    fn same_object(&self, py: Python<'_>, other: &MatrixParent) -> bool {
+        match (self, other) {
+            (Self::Matrix(a), Self::Matrix(b)) => a.bind(py).is(b.bind(py)),
+            (Self::Augmented2(a), Self::Augmented2(b)) => a.bind(py).is(b.bind(py)),
+            (Self::Augmented3(a), Self::Augmented3(b)) => a.bind(py).is(b.bind(py)),
+            _ => false,
+        }
+    }
+
+    /// Run `f` on the current inner `Matrix` mutably, holding the borrow for
+    /// the duration of the call.
+    fn with_matrix_mut<R>(
+        &self,
+        py: Python<'_>,
+        f: impl FnOnce(&mut RustMatrix) -> R,
+    ) -> PyResult<R> {
+        match self {
+            Self::Matrix(m) => Ok(f(&mut m.try_borrow_mut(py).map_err(borrow_error)?.0)),
+            Self::Augmented2(m) => Ok(f(m
+                .try_borrow_mut(py)
+                .map_err(borrow_error)?
+                .0
+                .get_mut()?)),
+            Self::Augmented3(m) => Ok(f(m
+                .try_borrow_mut(py)
+                .map_err(borrow_error)?
+                .0
+                .get_mut()?)),
+        }
+    }
+}
+
+/// The source backing a slice handle: either an owned vector, or a row of a
+/// matrix-like parent. In both cases we keep the parent Python object alive
+/// and store enough metadata to reconstruct the underlying Rust slice on
+/// each call, revalidating against the parent's current dimensions first.
+pub(crate) enum SliceParent {
+    Vector(Py<vectors::PyFpVector>),
+    MatrixRow { matrix: MatrixParent, row: usize },
+}
+
+impl SliceParent {
+    fn clone_ref(&self, py: Python<'_>) -> Self {
+        match self {
+            Self::Vector(v) => Self::Vector(v.clone_ref(py)),
+            Self::MatrixRow { matrix, row } => Self::MatrixRow {
+                matrix: matrix.clone_ref(py),
+                row: *row,
+            },
+        }
+    }
+
+    /// Whether `self` and `other` are backed by the same Python object, so that
+    /// taking a shared borrow of one while the other is mutably borrowed would
+    /// collide in PyO3.
+    ///
+    /// Two `MatrixRow`s of the same matrix object count as the same object
+    /// regardless of row index: the whole `Matrix` pyclass is borrowed as a
+    /// unit.
+    fn same_object(&self, py: Python<'_>, other: &SliceParent) -> bool {
+        match (self, other) {
+            (Self::Vector(a), Self::Vector(b)) => a.bind(py).is(b.bind(py)),
+            (Self::MatrixRow { matrix: a, .. }, Self::MatrixRow { matrix: b, .. }) => {
+                a.same_object(py, b)
+            }
+            _ => false,
+        }
+    }
+
+    /// Whether `self` is backed by the same Python object as the bound
+    /// `FpVector` `other`.
+    fn same_vector(&self, py: Python<'_>, other: &Bound<'_, vectors::PyFpVector>) -> bool {
+        match self {
+            Self::Vector(v) => v.bind(py).is(other),
+            _ => false,
+        }
+    }
+}
+
+#[pymodule]
+#[pyo3(name = "fp")]
+pub mod fp_py {
+    #[pymodule_export]
+    pub use super::matrices::{
+        PyAffineSubspace, PyMatrix, PyMatrixSliceMut, PyQuasiInverse, PySubquotient, PySubspace,
+        PySubspaceVectorIterator,
+    };
+    use super::matrices::{PyAugmentedMatrix2, PyAugmentedMatrix3};
+    #[pymodule_export]
+    pub use super::vectors::{PyFpSlice, PyFpSliceMut, PyFpVector, PyFpVectorIterator};
+    use super::*;
+
+    #[pyclass(name = "Fp", frozen, from_py_object)]
+    #[derive(Clone, Copy)]
+    struct PyFp(DynFp);
+
+    #[pyclass(name = "SmallFq", frozen, from_py_object)]
+    #[derive(Clone, Copy)]
+    struct PySmallFq(DynSmallFq);
+
+    #[pyclass(name = "FieldElement", frozen, from_py_object)]
+    #[derive(Clone, Copy)]
+    struct PyFieldElement(DynFieldElement);
+
+    /// The value-equality field types (`PyFp`, `PySmallFq`, `PyFieldElement`)
+    /// all expose byte-identical `__richcmp__`/`__hash__` methods: equality
+    /// only compares against another instance of the same class (`PyRef<Self>`)
+    /// on the wrapped value, `Eq`/`Ne` are the only supported operators, and
+    /// hashes go through the shared `py_hash` helper.
+    ///
+    /// This crate does not enable PyO3's `multiple-pymethods` feature, so each
+    /// class may have only one `#[pymethods]` block, and PyO3's `#[pymethods]`
+    /// proc-macro additionally rejects `macro_rules!` invocations *inside* the
+    /// block. So this macro emits the entire `#[pymethods]` block.
+    macro_rules! eq_hash_pymethods {
+        ($ty:ident, { $($extra:tt)* }) => {
+            #[pymethods]
+            impl $ty {
+                $($extra)*
+
+                pub fn __richcmp__(&self, other: &Bound<'_, PyAny>, op: CompareOp) -> bool {
+                    let eq = other
+                        .extract::<PyRef<Self>>()
+                        .is_ok_and(|other| self.0 == other.0);
+                    match op {
+                        CompareOp::Eq => eq,
+                        CompareOp::Ne => !eq,
+                        _ => false,
+                    }
+                }
+
+                pub fn __hash__(&self) -> isize {
+                    py_hash(&self.0)
+                }
+            }
+        };
+    }
+
+    eq_hash_pymethods!(PyFp, {
+        #[new]
+        pub fn new(p: u32) -> PyResult<Self> {
+            Ok(Self(DynFp::new(valid_prime(p)?)))
+        }
+
+        pub fn characteristic(&self) -> u32 {
+            self.0.characteristic().as_u32()
+        }
+
+        pub fn degree(&self) -> u32 {
+            self.0.degree()
+        }
+
+        pub fn zero(&self) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::Fp(self.0.zero()))
+        }
+
+        pub fn one(&self) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::Fp(self.0.one()))
+        }
+
+        pub fn element(&self, value: u32) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::Fp(self.0.element(value)))
+        }
+
+        pub fn __repr__(&self) -> String {
+            format!("Fp({})", self.characteristic())
+        }
+    });
+
+    eq_hash_pymethods!(PySmallFq, {
+        #[new]
+        pub fn new(p: u32, degree: u32) -> PyResult<Self> {
+            Ok(Self(small_fq(p, degree)?))
+        }
+
+        pub fn p(&self) -> u32 {
+            self.0.characteristic().as_u32()
+        }
+
+        pub fn degree(&self) -> u32 {
+            self.0.degree()
+        }
+
+        pub fn a(&self) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::SmallFq(self.0.a()))
+        }
+
+        pub fn q(&self) -> u32 {
+            self.0.q()
+        }
+
+        pub fn zero(&self) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::SmallFq(self.0.zero()))
+        }
+
+        pub fn one(&self) -> PyFieldElement {
+            PyFieldElement(DynFieldElement::SmallFq(self.0.one()))
+        }
+
+        pub fn __repr__(&self) -> String {
+            format!("SmallFq({}, {})", self.p(), self.degree())
+        }
+    });
+
+    eq_hash_pymethods!(PyFieldElement, {
+        pub fn inv(&self) -> Option<Self> {
+            self.0.inv().map(Self)
+        }
+
+        pub fn frobenius(&self) -> Self {
+            Self(self.0.frobenius())
+        }
+
+        pub fn field<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyAny>> {
+            match self.0 {
+                DynFieldElement::Fp(x) => {
+                    Py::new(py, PyFp(x.field())).map(|x| x.into_bound(py).into_any())
+                }
+                DynFieldElement::SmallFq(x) => {
+                    Py::new(py, PySmallFq(x.field())).map(|x| x.into_bound(py).into_any())
+                }
+            }
+        }
+
+        pub fn __add__(&self, rhs: Self) -> PyResult<Self> {
+            (self.0 + rhs.0)
+                .map(Self)
+                .ok_or_else(|| mismatched_field_error(self.0, rhs.0))
+        }
+
+        pub fn __sub__(&self, rhs: Self) -> PyResult<Self> {
+            (self.0 - rhs.0)
+                .map(Self)
+                .ok_or_else(|| mismatched_field_error(self.0, rhs.0))
+        }
+
+        pub fn __mul__(&self, rhs: Self) -> PyResult<Self> {
+            (self.0 * rhs.0)
+                .map(Self)
+                .ok_or_else(|| mismatched_field_error(self.0, rhs.0))
+        }
+
+        pub fn __truediv__(&self, rhs: Self) -> PyResult<Self> {
+            self.0.try_div(rhs.0).map(Self).map_err(|e| match e {
+                DivError::MismatchedField => mismatched_field_error(self.0, rhs.0),
+                DivError::DivisionByZero => PyZeroDivisionError::new_err("division by zero"),
+            })
+        }
+
+        pub fn __neg__(&self) -> Self {
+            Self(-self.0)
+        }
+
+        pub fn __int__(&self) -> PyResult<u32> {
+            self.0.try_as_u32().ok_or_else(|| {
+                PyValueError::new_err("SmallFq elements do not have a canonical integer value")
+            })
+        }
+
+        pub fn __repr__(&self) -> String {
+            match self.0 {
+                DynFieldElement::Fp(x) => {
+                    format!("FieldElement(Fp({}), {x})", x.field().characteristic())
+                }
+                DynFieldElement::SmallFq(x) => {
+                    let f = x.field();
+                    format!(
+                        "FieldElement(SmallFq({}, {}), {x})",
+                        f.characteristic(),
+                        f.degree()
+                    )
+                }
+            }
+        }
+    });
+
+    #[pymodule_init]
+    fn init(m: &Bound<'_, PyModule>) -> PyResult<()> {
+        // The `AugmentedMatrix2`/`AugmentedMatrix3` classes are produced by a
+        // macro, which `#[pymodule]` cannot auto-collect.
+        m.add_class::<PyAugmentedMatrix2>()?;
+        m.add_class::<PyAugmentedMatrix3>()?;
+        Ok(())
+    }
+}
diff --git a/ext_py/src/fp_mod/matrices.rs b/ext_py/src/fp_mod/matrices.rs
new file mode 100644
index 0000000000..2b0c01d7a4
--- /dev/null
+++ b/ext_py/src/fp_mod/matrices.rs
@@ -0,0 +1,1296 @@
+use std::io::Cursor;
+
+use fp::{
+    matrix::{
+        AffineSubspace as RustAffineSubspace, AugmentedMatrix as RustAugmentedMatrix,
+        Matrix as RustMatrix, MatrixSliceMut as RustMatrixSliceMut,
+        QuasiInverse as RustQuasiInverse, Subquotient as RustSubquotient, Subspace as RustSubspace,
+    },
+    prime::Prime,
+    vector::FpVector as RustFpVector,
+};
+use pyo3::{
+    exceptions::{PyIndexError, PyValueError},
+    types::PyBytes,
+};
+
+use super::{
+    vectors::{with_input_slice, with_target_slice_mut, PyFpSlice, PyFpSliceMut, PyFpVector},
+    *,
+};
+
+#[pyclass(name = "MatrixSliceMut")]
+pub struct PyMatrixSliceMut {
+    pub(crate) parent: MatrixParent,
+    row_start: usize,
+    row_end: usize,
+    col_start: usize,
+    col_end: usize,
+}
+
+#[pyclass(name = "Matrix")]
+pub struct PyMatrix(pub(crate) RustMatrix);
+
+#[pyclass(name = "Subspace")]
+pub struct PySubspace(RustSubspace);
+
+#[pyclass(name = "QuasiInverse")]
+pub struct PyQuasiInverse(RustQuasiInverse);
+
+#[pyclass(name = "Subquotient")]
+pub struct PySubquotient(RustSubquotient);
+
+#[pyclass(name = "AffineSubspace")]
+pub struct PyAffineSubspace(RustAffineSubspace);
+
+/// Lazy iterator over every vector in a subspace.
+///
+/// The upstream `Subspace::iter_all_vectors` iterator borrows the subspace,
+/// so it cannot be stored alongside an owned subspace in a `#[pyclass]`
+/// without a self-referential struct. Instead we own a clone of the
+/// subspace and an index counter, regenerating the i-th vector on each
+/// `__next__` from the base-`p` decomposition of the index.
+#[pyclass(name = "SubspaceVectorIterator")]
+pub struct PySubspaceVectorIterator {
+    subspace: RustSubspace,
+    index: u128,
+    total: u128,
+}
+
+/// Validate a `row_start..row_end` x `col_start..col_end` rectangle against
+/// a matrix's current `rows` x `columns`, raising `IndexError` otherwise.
+pub(crate) fn checked_rect(
+    row_start: usize,
+    row_end: usize,
+    col_start: usize,
+    col_end: usize,
+    rows: usize,
+    columns: usize,
+) -> PyResult<()> {
+    if row_start <= row_end && row_end <= rows && col_start <= col_end && col_end <= columns {
+        Ok(())
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "rectangle [{row_start}..{row_end}] x [{col_start}..{col_end}] out of range for \
+             matrix with {rows} rows and {columns} columns"
+        )))
+    }
+}
+
+/// Validate that `seg` is a segment index in `0..n`.
+fn checked_segment(seg: usize, n: usize) -> PyResult<()> {
+    if seg < n {
+        Ok(())
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "segment {seg} out of range for {n} segments"
+        )))
+    }
+}
+
+/// Validate a `[start, end]` segment-index range against an augmented
+/// matrix and return the width (column count) of the spanned rectangle.
+fn segment_cols<const N: usize>(
+    m: &RustAugmentedMatrix<N>,
+    start: usize,
+    end: usize,
+) -> PyResult<usize> {
+    checked_segment(start, N)?;
+    checked_segment(end, N)?;
+    let lo = m.start[start];
+    let hi = m.end[end];
+    if lo > hi {
+        return Err(PyValueError::new_err(format!(
+            "segment range [{start}, {end}] is empty or inverted"
+        )));
+    }
+    Ok(hi - lo)
+}
+
+impl PyMatrixSliceMut {
+    /// Number of rows spanned by the rectangle (cached; `with_slice_mut`
+    /// revalidates against the parent before any data access).
+    fn rows_span(&self) -> usize {
+        self.row_end - self.row_start
+    }
+
+    /// Number of columns spanned by the rectangle (cached; see `rows_span`).
+    fn cols_span(&self) -> usize {
+        self.col_end - self.col_start
+    }
+
+    /// Run `f` on the reconstructed `MatrixSliceMut`, after revalidating the
+    /// rectangle against the parent's current dimensions.
+    fn with_slice_mut<R>(
+        &self,
+        py: Python<'_>,
+        f: impl FnOnce(RustMatrixSliceMut<'_>) -> R,
+    ) -> PyResult<R> {
+        self.parent.with_matrix_mut(py, |m| {
+            checked_rect(
+                self.row_start,
+                self.row_end,
+                self.col_start,
+                self.col_end,
+                m.rows(),
+                m.columns(),
+            )?;
+            Ok(f(m.slice_mut(
+                self.row_start,
+                self.row_end,
+                self.col_start,
+                self.col_end,
+            )))
+        })?
+    }
+}
+
+#[pymethods]
+impl PyMatrixSliceMut {
+    #[getter]
+    pub fn prime(&self, py: Python<'_>) -> PyResult<u32> {
+        self.with_slice_mut(py, |s| s.prime().as_u32())
+    }
+
+    #[getter]
+    pub fn rows(&self, py: Python<'_>) -> PyResult<usize> {
+        self.with_slice_mut(py, |s| s.rows())
+    }
+
+    #[getter]
+    pub fn columns(&self, py: Python<'_>) -> PyResult<usize> {
+        self.with_slice_mut(py, |s| s.columns())
+    }
+
+    /// Return an immutable `FpSlice` over row `i` of the rectangle. The handle
+    /// revalidates against the parent on use.
+    pub fn row(&self, py: Python<'_>, i: usize) -> PyResult<PyFpSlice> {
+        let row = checked_row(i, self.rows_span())? + self.row_start;
+        Ok(PyFpSlice {
+            parent: SliceParent::MatrixRow {
+                matrix: self.parent.clone_ref(py),
+                row,
+            },
+            start: self.col_start,
+            end: self.col_end,
+        })
+    }
+
+    /// Return a mutable `FpSliceMut` over row `i` of the rectangle; mutating
+    /// it writes through to the parent matrix.
+    pub fn row_mut(&self, py: Python<'_>, i: usize) -> PyResult<PyFpSliceMut> {
+        let row = checked_row(i, self.rows_span())? + self.row_start;
+        Ok(PyFpSliceMut {
+            parent: SliceParent::MatrixRow {
+                matrix: self.parent.clone_ref(py),
+                row,
+            },
+            start: self.col_start,
+            end: self.col_end,
+        })
+    }
+
+    /// Restrict the view to rows `row_start..row_end`, returning a new
+    /// `MatrixSliceMut` over the same columns and parent.
+    pub fn row_slice(&self, py: Python<'_>, row_start: usize, row_end: usize) -> PyResult<Self> {
+        checked_range(row_start, row_end, self.rows_span())?;
+        Ok(Self {
+            parent: self.parent.clone_ref(py),
+            row_start: self.row_start + row_start,
+            row_end: self.row_start + row_end,
+            col_start: self.col_start,
+            col_end: self.col_end,
+        })
+    }
+
+    /// Return immutable `FpSlice` handles for every row of the rectangle.
+    ///
+    /// We materialize a list of row handles. TODO: Don't do that.
+    pub fn iter(&self, py: Python<'_>) -> PyResult<Vec<PyFpSlice>> {
+        (0..self.rows_span()).map(|i| self.row(py, i)).collect()
+    }
+
+    /// Return mutable `FpSliceMut` handles for every row of the rectangle.
+    ///
+    /// As with `iter`, this is an eager list of index-based row handles
+    /// rather than a lazy borrowing iterator. Mutating any handle writes
+    /// through to the parent matrix.
+    /// TODO: don't make a Vec of row_muts?
+    pub fn iter_mut(&self, py: Python<'_>) -> PyResult<Vec<PyFpSliceMut>> {
+        (0..self.rows_span()).map(|i| self.row_mut(py, i)).collect()
+    }
+
+    /// Add an identity matrix into the matrix. Requires a square matrix
+    /// otherwise a `ValueError` is raised
+    pub fn add_identity(&self, py: Python<'_>) -> PyResult<()> {
+        if self.rows_span() != self.cols_span() {
+            return Err(PyValueError::new_err(format!(
+                "add_identity requires a square rectangle: {} rows but {} columns",
+                self.rows_span(),
+                self.cols_span()
+            )));
+        }
+        self.with_slice_mut(py, |mut s| s.add_identity())
+    }
+
+    fn to_py(&self, py: Python<'_>) -> PyResult<Vec<Vec<u32>>> {
+        Ok(self.with_slice_mut(py, |s| s.to_vec())?)
+    }
+
+    pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
+        let (prime, rows, columns) =
+            self.with_slice_mut(py, |s| (s.prime().as_u32(), s.rows(), s.columns()))?;
+        Ok(format!("MatrixSliceMut({prime}, {rows}x{columns})"))
+    }
+}
+
+#[pymethods]
+impl PyMatrix {
+    #[new]
+    pub fn new(p: u32, rows: usize, columns: usize) -> PyResult<Self> {
+        Ok(Self(RustMatrix::new(valid_prime(p)?, rows, columns)))
+    }
+
+    #[staticmethod]
+    pub fn from_rows(p: u32, rows: Vec<PyRef<'_, PyFpVector>>, columns: usize) -> PyResult<Self> {
+        let p = valid_prime(p)?;
+        for row in &rows {
+            checked_same_prime(row.0.prime().as_u32(), p.as_u32())?;
+            checked_equal_len(row.0.len(), columns)?;
+        }
+        let input = rows.iter().map(|row| row.0.clone()).collect();
+        Ok(Self(RustMatrix::from_rows(p, input, columns)))
+    }
+
+    #[staticmethod]
+    pub fn from_row(p: u32, row: PyRef<'_, PyFpVector>, columns: usize) -> PyResult<Self> {
+        let p = valid_prime(p)?;
+        checked_same_prime(row.0.prime().as_u32(), p.as_u32())?;
+        checked_equal_len(row.0.len(), columns)?;
+        Ok(Self(RustMatrix::from_row(p, row.0.clone(), columns)))
+    }
+
+    #[staticmethod]
+    pub fn from_py(p: u32, input: Vec<Vec<u32>>) -> PyResult<Self> {
+        let p = valid_prime(p)?;
+        if let Some(first) = input.first() {
+            let columns = first.len();
+            for row in &input {
+                checked_equal_len(row.len(), columns)?;
+            }
+        }
+        Ok(Self(RustMatrix::from_vec(p, &input)))
+    }
+
+    #[staticmethod]
+    pub fn identity(p: u32, dim: usize) -> PyResult<Self> {
+        Ok(Self(RustMatrix::identity(valid_prime(p)?, dim)))
+    }
+
+    #[staticmethod]
+    pub fn from_bytes(p: u32, rows: usize, columns: usize, data: &[u8]) -> PyResult<Self> {
+        RustMatrix::from_bytes(valid_prime(p)?, rows, columns, &mut Cursor::new(data))
+            .map(Self)
+            .map_err(io_err)
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.prime().as_u32()
+    }
+
+    #[getter]
+    pub fn rows(&self) -> usize {
+        self.0.rows()
+    }
+
+    #[getter]
+    pub fn columns(&self) -> usize {
+        self.0.columns()
+    }
+
+    #[getter]
+    pub fn pivots(&self) -> Vec<isize> {
+        self.0.pivots().to_vec()
+    }
+
+    #[getter]
+    pub fn is_zero(&self) -> bool {
+        self.0.is_zero()
+    }
+
+    pub fn to_py(&self) -> Vec<Vec<u32>> {
+        self.0.to_vec()
+    }
+
+    pub fn to_bytes<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
+        serialize_to_pybytes(py, |buffer| self.0.to_bytes(buffer))
+    }
+
+    pub fn row(slf: PyRef<'_, Self>, row: usize) -> PyResult<PyFpSlice> {
+        checked_row(row, slf.0.rows())?;
+        let end = slf.0.columns();
+        let py = slf.py();
+        Ok(PyFpSlice {
+            parent: SliceParent::MatrixRow {
+                matrix: MatrixParent::Matrix(slf.into_pyobject(py)?.unbind()),
+                row,
+            },
+            start: 0,
+            end,
+        })
+    }
+
+    pub fn row_mut(slf: PyRef<'_, Self>, row: usize) -> PyResult<PyFpSliceMut> {
+        checked_row(row, slf.0.rows())?;
+        let end = slf.0.columns();
+        let py = slf.py();
+        Ok(PyFpSliceMut {
+            parent: SliceParent::MatrixRow {
+                matrix: MatrixParent::Matrix(slf.into_pyobject(py)?.unbind()),
+                row,
+            },
+            start: 0,
+            end,
+        })
+    }
+
+    /// Return mutable `FpSliceMut` handles for every row of the matrix,
+    /// one per row `i in 0..rows`.
+    pub fn iter_mut(slf: PyRef<'_, Self>) -> PyResult<Vec<PyFpSliceMut>> {
+        let py = slf.py();
+        let rows = slf.0.rows();
+        let end = slf.0.columns();
+        let parent = slf.into_pyobject(py)?.unbind();
+        Ok((0..rows)
+            .map(|row| PyFpSliceMut {
+                parent: SliceParent::MatrixRow {
+                    matrix: MatrixParent::Matrix(parent.clone_ref(py)),
+                    row,
+                },
+                start: 0,
+                end,
+            })
+            .collect())
+    }
+
+    pub fn set_to_zero(&mut self) {
+        self.0.set_to_zero()
+    }
+
+    pub fn assign(&mut self, other: &Self) -> PyResult<()> {
+        checked_same_prime(self.0.prime().as_u32(), other.0.prime().as_u32())?;
+        checked_equal_len(self.0.rows(), other.0.rows())?;
+        checked_equal_len(self.0.columns(), other.0.columns())?;
+        self.0.assign(&other.0);
+        Ok(())
+    }
+
+    pub fn swap_rows(&mut self, i: usize, j: usize) -> PyResult<()> {
+        checked_row(i, self.0.rows())?;
+        checked_row(j, self.0.rows())?;
+        self.0.swap_rows(i, j);
+        Ok(())
+    }
+
+    pub fn row_op(&mut self, target: usize, source: usize, c: u32) -> PyResult<()> {
+        checked_row(target, self.0.rows())?;
+        checked_row(source, self.0.rows())?;
+        if target == source {
+            return Err(PyValueError::new_err(
+                "target and source rows must be distinct",
+            ));
+        }
+        self.0.safe_row_op(target, source, c);
+        Ok(())
+    }
+
+    pub fn initialize_pivots(&mut self) {
+        self.0.initialize_pivots()
+    }
+
+    pub fn extend_column_dimension(&mut self, columns: usize) {
+        self.0.extend_column_dimension(columns)
+    }
+
+    pub fn extend_column_capacity(&mut self, columns: usize) {
+        self.0.extend_column_capacity(columns)
+    }
+
+    pub fn add_row(slf: PyRef<'_, Self>) -> PyResult<PyFpSliceMut> {
+        let py = slf.py();
+        let parent = slf.into_pyobject(py)?.unbind();
+        let (row, end) = {
+            let mut matrix = parent.try_borrow_mut(py).map_err(borrow_error)?;
+            matrix.0.add_row();
+            (matrix.0.rows() - 1, matrix.0.columns())
+        };
+        Ok(PyFpSliceMut {
+            parent: SliceParent::MatrixRow {
+                matrix: MatrixParent::Matrix(parent),
+                row,
+            },
+            start: 0,
+            end,
+        })
+    }
+
+    /// Return a mutable rectangular view over rows `row_start..row_end` and
+    /// columns `col_start..col_end`.
+    pub fn slice_mut(
+        slf: PyRef<'_, Self>,
+        row_start: usize,
+        row_end: usize,
+        col_start: usize,
+        col_end: usize,
+    ) -> PyResult<PyMatrixSliceMut> {
+        checked_rect(
+            row_start,
+            row_end,
+            col_start,
+            col_end,
+            slf.0.rows(),
+            slf.0.columns(),
+        )?;
+        let py = slf.py();
+        Ok(PyMatrixSliceMut {
+            parent: MatrixParent::Matrix(slf.into_pyobject(py)?.unbind()),
+            row_start,
+            row_end,
+            col_start,
+            col_end,
+        })
+    }
+
+    pub fn trim(&mut self, row_start: usize, row_end: usize, col_start: usize) -> PyResult<()> {
+        checked_range(row_start, row_end, self.0.rows())?;
+        if col_start > self.0.columns() {
+            return Err(PyIndexError::new_err(format!(
+                "column {col_start} out of range for matrix with {} columns",
+                self.0.columns()
+            )));
+        }
+        self.0.trim(row_start, row_end, col_start, false);
+        Ok(())
+    }
+
+    pub fn rotate_down(&mut self, start: usize, end: usize, shift: usize) -> PyResult<()> {
+        checked_range(start, end, self.0.rows())?;
+        if shift > end - start {
+            return Err(PyValueError::new_err(format!(
+                "shift {shift} exceeds range length {}",
+                end - start
+            )));
+        }
+        self.0.rotate_down(start..end, shift);
+        Ok(())
+    }
+
+    pub fn row_reduce(&mut self) -> usize {
+        self.0.row_reduce()
+    }
+
+    pub fn __len__(&self) -> usize {
+        self.0.rows()
+    }
+
+    pub fn __getitem__(slf: PyRef<'_, Self>, row: usize) -> PyResult<PyFpSlice> {
+        Self::row(slf, row)
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!("Matrix({}, {})", self.prime(), self.0)
+    }
+}
+
+impl PySubspace {
+    /// Validate that `vector` matches this subspace's prime and ambient
+    /// dimension, returning an error otherwise.
+    fn check_compatible(&self, vector: &RustFpVector) -> PyResult<()> {
+        checked_same_prime(self.0.prime().as_u32(), vector.prime().as_u32())?;
+        checked_equal_len(vector.len(), self.0.ambient_dimension())?;
+        Ok(())
+    }
+
+    /// Validate that `other` matches this subspace's prime and ambient
+    /// dimension, returning an error otherwise.
+    fn check_compatible_space(&self, other: &RustSubspace) -> PyResult<()> {
+        checked_same_prime(self.0.prime().as_u32(), other.prime().as_u32())?;
+        checked_equal_len(self.0.ambient_dimension(), other.ambient_dimension())?;
+        Ok(())
+    }
+}
+
+#[pymethods]
+impl PySubspace {
+    #[new]
+    pub fn new(p: u32, dim: usize) -> PyResult<Self> {
+        Ok(Self(RustSubspace::new(valid_prime(p)?, dim)))
+    }
+
+    #[staticmethod]
+    pub fn from_matrix(matrix: &PyMatrix) -> Self {
+        Self(RustSubspace::from_matrix(matrix.0.clone()))
+    }
+
+    #[staticmethod]
+    pub fn entire_space(p: u32, dim: usize) -> PyResult<Self> {
+        Ok(Self(RustSubspace::entire_space(valid_prime(p)?, dim)))
+    }
+
+    #[staticmethod]
+    pub fn from_bytes(p: u32, data: &[u8]) -> PyResult<Self> {
+        RustSubspace::from_bytes(valid_prime(p)?, &mut Cursor::new(data))
+            .map(Self)
+            .map_err(io_err)
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.prime().as_u32()
+    }
+
+    #[getter]
+    pub fn dimension(&self) -> usize {
+        self.0.dimension()
+    }
+
+    #[getter]
+    pub fn ambient_dimension(&self) -> usize {
+        self.0.ambient_dimension()
+    }
+
+    /// Test whether `vector` lies in this subspace.
+    pub fn contains(&self, py: Python<'_>, vector: &Bound<'_, PyAny>) -> PyResult<bool> {
+        with_input_slice(py, vector, |slice| {
+            checked_same_prime(self.0.prime().as_u32(), slice.prime().as_u32())?;
+            checked_equal_len(slice.len(), self.0.ambient_dimension())?;
+            Ok(self.0.contains(slice))
+        })
+    }
+
+    pub fn contains_space(&self, other: &Self) -> PyResult<bool> {
+        self.check_compatible_space(&other.0)?;
+        Ok(self.0.contains_space(&other.0))
+    }
+
+    pub fn add_vector(&mut self, vector: &PyFpVector) -> PyResult<usize> {
+        self.check_compatible(&vector.0)?;
+        Ok(self.0.add_vector(vector.0.as_slice()))
+    }
+
+    /// Reduce `vector` in place against this subspace
+    pub fn reduce(&self, vector: &mut PyFpVector) -> PyResult<()> {
+        self.check_compatible(&vector.0)?;
+        self.0.reduce(vector.0.as_slice_mut());
+        Ok(())
+    }
+
+    pub fn sum(&self, other: &Self) -> PyResult<Self> {
+        checked_same_prime(self.0.prime().as_u32(), other.0.prime().as_u32())?;
+        checked_equal_len(self.0.ambient_dimension(), other.0.ambient_dimension())?;
+        Ok(Self(self.0.sum(&other.0)))
+    }
+
+    /// Return the basis of the subspace as a list of owned `FpVector`s.
+    pub fn iter(&self) -> Vec<PyFpVector> {
+        self.0
+            .iter()
+            .map(|row| PyFpVector(row.to_owned()))
+            .collect()
+    }
+
+    /// Return the basis of the subspace as a list of owned `FpVector`s.
+    /// Mirrors upstream `Subspace::basis`.
+    #[getter]
+    pub fn basis(&self) -> Vec<PyFpVector> {
+        self.iter()
+    }
+
+    /// Return a lazy iterator over every vector in the subspace.
+    pub fn iter_all_vectors(&self) -> PySubspaceVectorIterator {
+        let p = u128::from(self.0.prime().as_u32());
+        let dim = self.0.dimension() as u32;
+        let total = p.checked_pow(dim).unwrap_or(u128::MAX);
+        PySubspaceVectorIterator {
+            subspace: self.0.clone(),
+            index: 0,
+            total,
+        }
+    }
+
+    pub fn set_to_zero(&mut self) {
+        self.0.set_to_zero()
+    }
+
+    pub fn set_to_entire(&mut self) {
+        self.0.set_to_entire()
+    }
+
+    pub fn to_bytes<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
+        serialize_to_pybytes(py, |buffer| self.0.to_bytes(buffer))
+    }
+
+    pub fn __len__(&self) -> usize {
+        self.0.dimension()
+    }
+
+    pub fn __contains__(&self, py: Python<'_>, vector: &Bound<'_, PyAny>) -> PyResult<bool> {
+        self.contains(py, vector)
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!(
+            "Subspace({}, dim={}, ambient={})",
+            self.prime(),
+            self.0.dimension(),
+            self.0.ambient_dimension()
+        )
+    }
+}
+
+#[pymethods]
+impl PyQuasiInverse {
+    /// Construct a `QuasiInverse` from an optional `image` (pivot list) and a
+    /// `preimage` matrix.
+    #[new]
+    #[pyo3(signature = (image, preimage))]
+    pub fn new(image: Option<Vec<isize>>, preimage: &PyMatrix) -> Self {
+        Self(RustQuasiInverse::new(image, preimage.0.clone()))
+    }
+
+    /// Deserialize a `QuasiInverse` from bytes produced by [`Self::to_bytes`].
+    #[staticmethod]
+    pub fn from_bytes(p: u32, data: &[u8]) -> PyResult<Self> {
+        RustQuasiInverse::from_bytes(valid_prime(p)?, &mut Cursor::new(data))
+            .map(Self)
+            .map_err(io_err)
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.prime().as_u32()
+    }
+
+    #[getter]
+    pub fn image_dimension(&self) -> usize {
+        self.0.image_dimension()
+    }
+
+    #[getter]
+    pub fn source_dimension(&self) -> usize {
+        self.0.source_dimension()
+    }
+
+    #[getter]
+    pub fn target_dimension(&self) -> usize {
+        self.0.target_dimension()
+    }
+
+    #[getter]
+    pub fn preimage(&self) -> PyMatrix {
+        PyMatrix(self.0.preimage().clone())
+    }
+
+    #[getter]
+    pub fn pivots(&self) -> Option<Vec<isize>> {
+        self.0.pivots().map(<[isize]>::to_vec)
+    }
+
+    /// Apply the quasi-inverse to `input` and add `coeff` times the result
+    /// to `target`.
+    ///
+    /// `input` is a vector in the target space (length `target_dimension`)
+    /// and `target` receives the result in the source space (length
+    /// `source_dimension`). Both accept either an `FpVector` or the
+    /// corresponding slice handle.
+    pub fn apply(
+        &self,
+        py: Python<'_>,
+        target: &Bound<'_, PyAny>,
+        coeff: u32,
+        input: &Bound<'_, PyAny>,
+    ) -> PyResult<()> {
+        // Borrow the input. If the same object is passed as both `input` and
+        // `target`, the nested shared+mutable borrows raises `RuntimeError`.
+        with_input_slice(py, input, |input_slice| {
+            checked_same_prime(self.0.prime().as_u32(), input_slice.prime().as_u32())?;
+            checked_equal_len(input_slice.len(), self.0.target_dimension())?;
+            with_target_slice_mut(py, target, |target_slice| {
+                checked_same_prime(
+                    self.0.prime().as_u32(),
+                    target_slice.as_slice().prime().as_u32(),
+                )?;
+                checked_equal_len(target_slice.as_slice().len(), self.0.source_dimension())?;
+                // Reduce `coeff` mod p before calling upstream to avoid overflow
+                let coeff = coeff % self.0.prime().as_u32();
+                self.0.apply(target_slice, coeff, input_slice);
+                Ok(())
+            })
+        })
+    }
+
+    /// Serialize the quasi-inverse to bytes.
+    pub fn to_bytes<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
+        serialize_to_pybytes(py, |buffer| self.0.to_bytes(buffer))
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!(
+            "QuasiInverse({}, image_dim={}, source_dim={}, target_dim={})",
+            self.prime(),
+            self.0.image_dimension(),
+            self.0.source_dimension(),
+            self.0.target_dimension()
+        )
+    }
+}
+
+impl PySubquotient {
+    /// Validate that `vector` matches this subquotient's prime and ambient
+    /// dimension, returning an error otherwise.
+    fn check_compatible(&self, vector: &RustFpVector) -> PyResult<()> {
+        checked_same_prime(self.0.prime().as_u32(), vector.prime().as_u32())?;
+        checked_equal_len(vector.len(), self.0.ambient_dimension())?;
+        Ok(())
+    }
+}
+
+#[pymethods]
+impl PySubquotient {
+    /// Create a new subquotient of an ambient space of dimension `dim`,
+    /// defaulting to the zero subspace.
+    #[new]
+    pub fn new(p: u32, dim: usize) -> PyResult<Self> {
+        Ok(Self(RustSubquotient::new(valid_prime(p)?, dim)))
+    }
+
+    /// Create a new subquotient of an ambient space of dimension `dim`,
+    /// where the subspace is the full space and the quotient is trivial.
+    #[staticmethod]
+    pub fn new_full(p: u32, dim: usize) -> PyResult<Self> {
+        Ok(Self(RustSubquotient::new_full(valid_prime(p)?, dim)))
+    }
+
+    /// Construct the subquotient `(sub + quotient) / quotient` from a chain
+    /// of subspaces. The two subspaces must share a prime and ambient
+    /// dimension.
+    #[staticmethod]
+    pub fn from_parts(sub: &PySubspace, quotient: &PySubspace) -> PyResult<Self> {
+        checked_same_prime(sub.0.prime().as_u32(), quotient.0.prime().as_u32())?;
+        checked_equal_len(sub.0.ambient_dimension(), quotient.0.ambient_dimension())?;
+        Ok(Self(RustSubquotient::from_parts(
+            sub.0.clone(),
+            quotient.0.clone(),
+        )))
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.prime().as_u32()
+    }
+
+    #[getter]
+    pub fn dimension(&self) -> usize {
+        self.0.dimension()
+    }
+
+    #[getter]
+    pub fn ambient_dimension(&self) -> usize {
+        self.0.ambient_dimension()
+    }
+
+    #[getter]
+    pub fn quotient_dimension(&self) -> usize {
+        self.0.quotient_dimension()
+    }
+
+    #[getter]
+    pub fn subspace_dimension(&self) -> usize {
+        self.0.subspace_dimension()
+    }
+
+    #[getter]
+    pub fn is_empty(&self) -> bool {
+        self.0.is_empty()
+    }
+
+    pub fn zeros(&self) -> PySubspace {
+        PySubspace(self.0.zeros().clone())
+    }
+
+    /// The generators of the subquotient, returned as a list of owned
+    /// `FpVector`s.
+    #[getter]
+    pub fn gens(&self) -> Vec<PyFpVector> {
+        self.0
+            .gens()
+            .map(|row| PyFpVector(row.to_owned()))
+            .collect()
+    }
+
+    /// The generators of the subspace part of the subquotient, returned as
+    /// a list of owned `FpVector`s.
+    #[getter]
+    pub fn subspace_gens(&self) -> Vec<PyFpVector> {
+        self.0
+            .subspace_gens()
+            .map(|row| PyFpVector(row.to_owned()))
+            .collect()
+    }
+
+    /// The pivot columns of the complement to the subspace.
+    #[getter]
+    pub fn complement_pivots(&self) -> Vec<usize> {
+        self.0.complement_pivots().collect()
+    }
+
+    /// The pivot table of the quotient subspace.
+    #[getter]
+    pub fn quotient_pivots(&self) -> Vec<isize> {
+        self.0.quotient_pivots().to_vec()
+    }
+
+    /// Reduce `vector` in place: project it onto a complement of the
+    /// quotient and express it relative to the generators. Returns the list
+    /// of coefficients with respect to the generators.
+    pub fn reduce(&self, vector: &mut PyFpVector) -> PyResult<Vec<u32>> {
+        self.check_compatible(&vector.0)?;
+        Ok(self.0.reduce(vector.0.as_slice_mut()))
+    }
+
+    /// Project `vector` in place onto the complement of the quotient part.
+    pub fn reduce_by_quotient(&self, py: Python<'_>, vector: &Bound<'_, PyAny>) -> PyResult<()> {
+        with_target_slice_mut(py, vector, |slice| {
+            checked_same_prime(self.0.prime().as_u32(), slice.as_slice().prime().as_u32())?;
+            checked_equal_len(slice.as_slice().len(), self.0.ambient_dimension())?;
+            self.0.reduce_by_quotient(slice);
+            Ok(())
+        })
+    }
+
+    /// Add `vector` to the quotient part of the subquotient.
+    pub fn quotient(&mut self, vector: &PyFpVector) -> PyResult<()> {
+        self.check_compatible(&vector.0)?;
+        self.0.quotient(vector.0.as_slice());
+        Ok(())
+    }
+
+    /// Add `vector` as a generator of the subquotient.
+    pub fn add_gen(&mut self, vector: &PyFpVector) -> PyResult<()> {
+        self.check_compatible(&vector.0)?;
+        self.0.add_gen(vector.0.as_slice());
+        Ok(())
+    }
+
+    /// Remove all generators, leaving the quotient part untouched.
+    pub fn clear_gens(&mut self) {
+        self.0.clear_gens()
+    }
+
+    /// Set the subquotient to be the full ambient space quotiented by zero.
+    pub fn set_to_full(&mut self) {
+        self.0.set_to_full()
+    }
+
+    /// Apply `matrix` to each generator of `source`, then reduce the image
+    /// in `target`, returning the coefficient lists.
+    #[staticmethod]
+    pub fn reduce_matrix(
+        matrix: &PyMatrix,
+        source: &Self,
+        target: &Self,
+    ) -> PyResult<Vec<Vec<u32>>> {
+        checked_same_prime(source.0.prime().as_u32(), target.0.prime().as_u32())?;
+        checked_same_prime(source.0.prime().as_u32(), matrix.0.prime().as_u32())?;
+        checked_equal_len(matrix.0.rows(), source.0.ambient_dimension())?;
+        checked_equal_len(matrix.0.columns(), target.0.ambient_dimension())?;
+        Ok(RustSubquotient::reduce_matrix(
+            &matrix.0, &source.0, &target.0,
+        ))
+    }
+
+    pub fn __len__(&self) -> usize {
+        self.0.dimension()
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!(
+            "Subquotient({}, dim={}, ambient={})",
+            self.prime(),
+            self.0.dimension(),
+            self.0.ambient_dimension()
+        )
+    }
+}
+
+impl PyAffineSubspace {
+    /// Validate that `other` matches this affine subspace's prime and
+    /// ambient dimension, returning an error otherwise.
+    fn check_compatible_space(&self, other: &Self) -> PyResult<()> {
+        checked_same_prime(self.prime(), other.prime())?;
+        checked_equal_len(self.ambient_dimension(), other.ambient_dimension())?;
+        Ok(())
+    }
+}
+
+#[pymethods]
+impl PyAffineSubspace {
+    /// Construct an affine subspace `offset + linear_part`.
+    #[new]
+    pub fn new(offset: &PyFpVector, linear_part: &PySubspace) -> Self {
+        Self(RustAffineSubspace::new(
+            offset.0.clone(),
+            linear_part.0.clone(),
+        ))
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.linear_part().prime().as_u32()
+    }
+
+    #[getter]
+    pub fn ambient_dimension(&self) -> usize {
+        self.0.linear_part().ambient_dimension()
+    }
+
+    #[getter]
+    pub fn dimension(&self) -> usize {
+        self.0.linear_part().dimension()
+    }
+
+    /// Return an owned copy of the reduced offset vector.
+    pub fn offset(&self) -> PyFpVector {
+        PyFpVector(self.0.offset().clone())
+    }
+
+    /// Return an owned copy of the linear part `Subspace`
+    pub fn linear_part(&self) -> PySubspace {
+        PySubspace(self.0.linear_part().clone())
+    }
+
+    /// Test whether `vector` lies in this affine subspace.
+    pub fn contains(&self, py: Python<'_>, vector: &Bound<'_, PyAny>) -> PyResult<bool> {
+        with_input_slice(py, vector, |slice| {
+            checked_same_prime(self.prime(), slice.prime().as_u32())?;
+            checked_equal_len(slice.len(), self.ambient_dimension())?;
+            Ok(self.0.contains(slice))
+        })
+    }
+
+    pub fn contains_space(&self, other: &Self) -> PyResult<bool> {
+        self.check_compatible_space(other)?;
+        Ok(self.0.contains_space(&other.0))
+    }
+
+    /// Return the affine subspace spanned by the union of `self` and
+    /// `other`: the sum of the linear parts translated by the sum of the
+    /// offsets.
+    pub fn sum(&self, other: &Self) -> PyResult<Self> {
+        self.check_compatible_space(other)?;
+        Ok(Self(self.0.sum(&other.0)))
+    }
+
+    pub fn __contains__(&self, py: Python<'_>, vector: &Bound<'_, PyAny>) -> PyResult<bool> {
+        self.contains(py, vector)
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!("AffineSubspace({})", self.0)
+    }
+}
+
+#[pymethods]
+impl PySubspaceVectorIterator {
+    pub fn __iter__(slf: PyRefMut<'_, Self>) -> PyRefMut<'_, Self> {
+        slf
+    }
+
+    pub fn __next__(&mut self) -> Option<PyFpVector> {
+        if self.index >= self.total {
+            return None;
+        }
+        let p = u128::from(self.subspace.prime().as_u32());
+        let dim = self.subspace.dimension();
+        // Decode `index` into base-`p` digits, most significant first, to
+        // match the lexicographic order of `combinations` upstream where
+        // the first digit (matching the first basis row) varies slowest.
+        let mut digits = vec![0u32; dim];
+        let mut rem = self.index;
+        for slot in digits.iter_mut().rev() {
+            *slot = (rem % p) as u32;
+            rem /= p;
+        }
+        let mut vector =
+            RustFpVector::new(self.subspace.prime(), self.subspace.ambient_dimension());
+        for (&c, row) in digits.iter().zip(self.subspace.iter()) {
+            vector.as_slice_mut().add(row, c);
+        }
+        self.index += 1;
+        Some(PyFpVector(vector))
+    }
+}
+
+/// `AugmentedMatrix<N>` is a const-generic type, and PyO3 cannot expose a
+/// generic `#[pyclass]`. We bind the two concrete arities used in the codebase
+/// (`N = 2` and `N = 3`) as separate classes `AugmentedMatrix2` and
+/// `AugmentedMatrix3`. To avoid duplicating the shared glue, this
+/// `macro_rules!` macro generates each class from a single definition. However,
+/// `#[pymodule]` cannot see through a `macro_rules!` expansion to auto-collect
+/// the classes, so they are registered explicitly with `add_class` in
+/// `#[pymodule_init]`.
+macro_rules! augmented_matrix_pyclass {
+    ($name:ident, $pyname:literal, $n:literal, $variant:ident, { $($extra:tt)* }) => {
+        /// The inner `AugmentedMatrix<N>` is held in a [`Consumable`] so the
+        /// consuming methods (`into_matrix`, `compute_quasi_inverses`) can
+        /// `take()` it out and run the upstream by-value operation.
+        #[pyclass(name = $pyname)]
+        pub struct $name(pub(crate) Consumable<RustAugmentedMatrix<$n>>);
+
+        #[pymethods]
+        impl $name {
+            /// Construct an `rows x sum(columns)` augmented matrix whose
+            /// column blocks have the given widths. `columns` must contain
+            /// exactly `N` segment widths.
+            #[new]
+            fn new(p: u32, rows: usize, columns: Vec<usize>) -> PyResult<Self> {
+                let len = columns.len();
+                let cols: [usize; $n] = columns.try_into().map_err(|_| {
+                    PyValueError::new_err(format!(
+                        "expected {} segment widths, got {len}",
+                        $n
+                    ))
+                })?;
+                Ok(Self(Consumable::new(
+                    $pyname,
+                    RustAugmentedMatrix::<$n>::new(valid_prime(p)?, rows, cols),
+                )))
+            }
+
+            #[getter]
+            fn prime(&self) -> PyResult<u32> {
+                Ok(self.0.get()?.prime().as_u32())
+            }
+
+            #[getter]
+            fn rows(&self) -> PyResult<usize> {
+                Ok(self.0.get()?.rows())
+            }
+
+            #[getter]
+            fn columns(&self) -> PyResult<usize> {
+                Ok(self.0.get()?.columns())
+            }
+
+            /// Number of column segments (`N`).
+            #[getter]
+            fn segments(&self) -> usize {
+                $n
+            }
+
+            /// The starting column index of each segment.
+            #[getter]
+            fn segment_starts(&self) -> PyResult<Vec<usize>> {
+                Ok(self.0.get()?.start.to_vec())
+            }
+
+            /// The (exclusive) ending column index of each segment.
+            #[getter]
+            fn segment_ends(&self) -> PyResult<Vec<usize>> {
+                Ok(self.0.get()?.end.to_vec())
+            }
+
+            #[getter]
+            fn pivots(&self) -> PyResult<Vec<isize>> {
+                Ok(self.0.get()?.pivots().to_vec())
+            }
+
+            #[getter]
+            fn is_zero(&self) -> PyResult<bool> {
+                Ok(self.0.get()?.is_zero())
+            }
+
+            fn to_py(&self) -> PyResult<Vec<Vec<u32>>> {
+                Ok(self.0.get()?.to_vec())
+            }
+
+            fn row_reduce(&mut self) -> PyResult<usize> {
+                Ok(self.0.get_mut()?.row_reduce())
+            }
+
+            /// Add an identity matrix into the rectangular segment spanning
+            /// segment indices `start..=end`. The segment must be square
+            fn add_identity(&mut self, start: usize, end: usize) -> PyResult<()> {
+                let m = self.0.get_mut()?;
+                let cols = segment_cols(m, start, end)?;
+                if m.rows() != cols {
+                    return Err(PyValueError::new_err(format!(
+                        "add_identity requires a square segment: matrix has {} rows but \
+                         segment [{start}, {end}] has {cols} columns",
+                        m.rows()
+                    )));
+                }
+                m.segment(start, end).add_identity();
+                Ok(())
+            }
+
+            /// Return an owned copy of row `i` restricted to the columns of the
+            /// segment range `start..=end`.
+            /// TODO: Maybe don't copy?
+            fn row_segment(
+                &self,
+                i: usize,
+                start: usize,
+                end: usize,
+            ) -> PyResult<PyFpVector> {
+                let m = self.0.get()?;
+                checked_row(i, m.rows())?;
+                segment_cols(m, start, end)?;
+                Ok(PyFpVector(m.row_segment(i, start, end).to_owned()))
+            }
+
+            /// Return a mutable rectangular view spanning all rows and the
+            /// columns of segment range `start..=end`, as a `MatrixSliceMut`
+            /// over the inner matrix.
+            fn segment(
+                slf: PyRef<'_, Self>,
+                start: usize,
+                end: usize,
+            ) -> PyResult<PyMatrixSliceMut> {
+                let (row_end, col_start, col_end) = {
+                    let m = slf.0.get()?;
+                    segment_cols(m, start, end)?;
+                    (m.rows(), m.start[start], m.end[end])
+                };
+                let py = slf.py();
+                Ok(PyMatrixSliceMut {
+                    parent: MatrixParent::$variant(slf.into_pyobject(py)?.unbind()),
+                    row_start: 0,
+                    row_end,
+                    col_start,
+                    col_end,
+                })
+            }
+
+            /// Return a `FpSliceMut` over row `i` restricted to the
+            /// columns of segment range `start..=end`.
+            fn row_segment_mut(
+                slf: PyRef<'_, Self>,
+                i: usize,
+                start: usize,
+                end: usize,
+            ) -> PyResult<PyFpSliceMut> {
+                let (col_start, col_end) = {
+                    let m = slf.0.get()?;
+                    checked_row(i, m.rows())?;
+                    segment_cols(m, start, end)?;
+                    (m.start[start], m.end[end])
+                };
+                let py = slf.py();
+                Ok(PyFpSliceMut {
+                    parent: SliceParent::MatrixRow {
+                        matrix: MatrixParent::$variant(slf.into_pyobject(py)?.unbind()),
+                        row: i,
+                    },
+                    start: col_start,
+                    end: col_end,
+                })
+            }
+
+            /// Return a `FpSliceMut` over the whole of row `i`
+            /// (all columns, across every segment).
+            fn row_mut(slf: PyRef<'_, Self>, i: usize) -> PyResult<PyFpSliceMut> {
+                let end = {
+                    let m = slf.0.get()?;
+                    checked_row(i, m.rows())?;
+                    m.columns()
+                };
+                let py = slf.py();
+                Ok(PyFpSliceMut {
+                    parent: SliceParent::MatrixRow {
+                        matrix: MatrixParent::$variant(slf.into_pyobject(py)?.unbind()),
+                        row: i,
+                    },
+                    start: 0,
+                    end,
+                })
+            }
+
+            /// Compute the kernel of the augmented matrix (which must be row
+            /// reduced), returning an owned `Subspace`. Available for all
+            /// arities.
+            fn compute_kernel(&self) -> PyResult<PySubspace> {
+                Ok(PySubspace(self.0
+                    .get()?
+                    .compute_kernel()))
+            }
+
+            /// Return the inner `Matrix` as an owned `Matrix`, **consuming**
+            /// this augmented matrix.
+            ///
+            /// After this call the augmented matrix is consumed, so any further
+            /// use raises `RuntimeError`.
+            #[allow(clippy::wrong_self_convention)]
+            fn into_matrix(&mut self) -> PyResult<PyMatrix> {
+                Ok(PyMatrix(self.0.take()?.into_matrix()))
+            }
+
+            fn __repr__(&self) -> String {
+                match self.0.get() {
+                    Ok(m) => format!(
+                        concat!($pyname, "({}, {}x{})"),
+                        m.prime().as_u32(),
+                        m.rows(),
+                        m.columns()
+                    ),
+                    Err(_) => concat!($pyname, "(consumed)").to_string(),
+                }
+            }
+
+            $($extra)*
+        }
+    };
+}
+
+augmented_matrix_pyclass!(PyAugmentedMatrix2, "AugmentedMatrix2", 2, Augmented2, {
+    /// Compute the image of the augmented matrix `[A | I]` (which must be
+    /// row reduced), returning an owned `Subspace`.
+    fn compute_image(&self) -> PyResult<PySubspace> {
+        Ok(PySubspace(self.0.get()?.compute_image()))
+    }
+
+    /// Compute the quasi-inverse of the augmented matrix `[A | I]` (which
+    /// must be row reduced), returning an owned `QuasiInverse`.
+    fn compute_quasi_inverse(&self) -> PyResult<PyQuasiInverse> {
+        Ok(PyQuasiInverse(self.0.get()?.compute_quasi_inverse()))
+    }
+
+    #[staticmethod]
+    pub fn from_py(p: u32, input: Vec<Vec<u32>>) -> PyResult<Self> {
+        let p = valid_prime(p)?;
+        if input.is_empty() {
+            return Err(PyValueError::new_err(
+                "AugmentedMatrix2.from_py() requires at least one row",
+            ));
+        }
+        let columns = input[0].len();
+        for row in &input {
+            checked_equal_len(row.len(), columns)?;
+        }
+        Ok(Self(Consumable::new("AugmentedMatrix2", RustMatrix::augmented_from_vec(p, &input))))
+    }
+});
+
+augmented_matrix_pyclass!(PyAugmentedMatrix3, "AugmentedMatrix3", 3, Augmented3, {
+    /// Compute the two quasi-inverses for a row-reduced augmented matrix of the
+    /// form `[A | B | I]` where `A` is surjective, returning the pair
+    /// `(quasi_inverse_of_A, residual_quasi_inverse)`.
+    ///
+    /// This consumes the augmented matrix. After this call any further use
+    /// raises `RuntimeError`.
+    fn compute_quasi_inverses(&mut self) -> PyResult<(PyQuasiInverse, PyQuasiInverse)> {
+        let (a, b) = self.0.take()?.compute_quasi_inverses();
+        Ok((PyQuasiInverse(a), PyQuasiInverse(b)))
+    }
+});
diff --git a/ext_py/src/fp_mod/vectors.rs b/ext_py/src/fp_mod/vectors.rs
new file mode 100644
index 0000000000..7df73fa782
--- /dev/null
+++ b/ext_py/src/fp_mod/vectors.rs
@@ -0,0 +1,701 @@
+use std::io::Cursor;
+
+use fp::{
+    prime::Prime,
+    vector::{FpSlice as RustFpSlice, FpSliceMut as RustFpSliceMut, FpVector as RustFpVector},
+};
+use pyo3::{exceptions::PyValueError, types::PyBytes};
+
+use super::*;
+
+/// Length of a slice-like operand argument (`FpVector` or `FpSlice`),
+/// computed without retaining a borrow. Used for pre-borrow dimension
+/// checks by the `FpSliceMut` operand-taking methods, which accept either
+/// an `FpVector` or an `FpSlice`.
+pub(crate) fn slice_like_len(operand: &Bound<'_, PyAny>) -> PyResult<usize> {
+    if let Ok(slice) = operand.extract::<PyRef<'_, PyFpSlice>>() {
+        Ok(slice.span())
+    } else if let Ok(vector) = operand.extract::<PyRef<'_, PyFpVector>>() {
+        Ok(vector.0.len())
+    } else {
+        Err(PyValueError::new_err("expected an FpVector or FpSlice"))
+    }
+}
+
+/// Run `f` on the reconstructed immutable slice for `parent[start..end]`,
+/// after revalidating the parent's current dimensions.
+///
+/// Revalidation only guards the parent's current *dimensions* (vector length
+/// or matrix row count and row length). It deliberately does not track
+/// logical-coordinate remapping: an operation like `Matrix::trim` with
+/// `col_start > 0` shifts the data backwards in each row without shrinking it
+/// below the slice's `end`, so a surviving slice silently reads the remapped
+/// columns rather than raising. Preventing that would require tracking the
+/// origin of every coordinate, which is out of scope for the
+/// handle+range design.
+pub(crate) fn with_parent_slice<R>(
+    parent: &SliceParent,
+    start: usize,
+    end: usize,
+    py: Python<'_>,
+    f: impl FnOnce(RustFpSlice<'_>) -> R,
+) -> PyResult<R> {
+    match parent {
+        SliceParent::Vector(v) => {
+            let parent = v.try_borrow(py).map_err(borrow_error)?;
+            checked_range(start, end, parent.0.len())?;
+            Ok(f(parent.0.slice(start, end)))
+        }
+        SliceParent::MatrixRow { matrix, row } => matrix.with_matrix(py, |m| {
+            checked_row(*row, m.rows())?;
+            let full = m.row(*row);
+            checked_range(start, end, full.len())?;
+            Ok(f(full.restrict(start, end)))
+        })?,
+    }
+}
+
+/// Run `f` on the reconstructed mutable slice for `parent[start..end]`,
+/// after revalidating the parent's current dimensions.
+pub(crate) fn with_parent_slice_mut<R>(
+    parent: &SliceParent,
+    start: usize,
+    end: usize,
+    py: Python<'_>,
+    f: impl FnOnce(RustFpSliceMut<'_>) -> R,
+) -> PyResult<R> {
+    match parent {
+        SliceParent::Vector(v) => {
+            let mut parent = v.try_borrow_mut(py).map_err(borrow_error)?;
+            checked_range(start, end, parent.0.len())?;
+            Ok(f(parent.0.slice_mut(start, end)))
+        }
+        SliceParent::MatrixRow { matrix, row } => matrix.with_matrix_mut(py, |m| {
+            checked_row(*row, m.rows())?;
+            // Validate against the actual current row length, matching the
+            // read path (`with_parent_slice`). For a `Matrix` this equals
+            // `columns()`, but deriving it from the row keeps both paths
+            // consistent regardless of that invariant.
+            let row_len = m.row(*row).len();
+            checked_range(start, end, row_len)?;
+            Ok(f(m.row_mut(*row).slice_mut(start, end)))
+        })?,
+    }
+}
+
+pub(crate) fn checked_index(index: usize, len: usize) -> PyResult<usize> {
+    if index < len {
+        Ok(index)
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "index {index} out of range for vector of length {len}"
+        )))
+    }
+}
+
+pub(crate) fn py_index(index: isize, len: usize) -> PyResult<usize> {
+    let index = if index < 0 {
+        len as isize + index
+    } else {
+        index
+    };
+    if index >= 0 && (index as usize) < len {
+        Ok(index as usize)
+    } else {
+        Err(PyIndexError::new_err(format!(
+            "index {index} out of range for vector of length {len}"
+        )))
+    }
+}
+
+#[pyclass(name = "FpVector")]
+pub struct PyFpVector(pub(crate) RustFpVector);
+
+#[pymethods]
+impl PyFpVector {
+    #[new]
+    pub fn new(p: u32, len: usize) -> PyResult<Self> {
+        Ok(Self(RustFpVector::new(valid_prime(p)?, len)))
+    }
+
+    /// Construct a new zero vector of length `len` over `F_p`. Static-method
+    /// alias for the constructor, allowing `FpVector.new(p, len)`.
+    #[staticmethod]
+    #[pyo3(name = "new")]
+    pub fn new_static(p: u32, len: usize) -> PyResult<Self> {
+        Self::new(p, len)
+    }
+
+    #[staticmethod]
+    pub fn new_with_capacity(p: u32, len: usize, capacity: usize) -> PyResult<Self> {
+        Ok(Self(RustFpVector::new_with_capacity(
+            valid_prime(p)?,
+            len,
+            capacity,
+        )))
+    }
+
+    #[staticmethod]
+    pub fn from_py(p: u32, entries: Vec<u32>) -> PyResult<Self> {
+        Ok(Self(RustFpVector::from_slice(valid_prime(p)?, &entries)))
+    }
+
+    #[staticmethod]
+    pub fn from_bytes(p: u32, len: usize, data: &[u8]) -> PyResult<Self> {
+        RustFpVector::from_bytes(valid_prime(p)?, len, &mut Cursor::new(data))
+            .map(Self)
+            .map_err(io_err)
+    }
+
+    #[getter]
+    pub fn prime(&self) -> u32 {
+        self.0.prime().as_u32()
+    }
+
+    #[getter]
+    pub fn len(&self) -> usize {
+        self.0.len()
+    }
+
+    #[getter]
+    pub fn is_empty(&self) -> bool {
+        self.0.is_empty()
+    }
+
+    pub fn entry(&self, index: usize) -> PyResult<u32> {
+        Ok(self.0.entry(checked_index(index, self.0.len())?))
+    }
+
+    #[getter]
+    pub fn density(&self) -> f32 {
+        self.0.density()
+    }
+
+    #[getter]
+    pub fn is_zero(&self) -> bool {
+        self.0.is_zero()
+    }
+
+    #[getter]
+    pub fn first_nonzero(&self) -> Option<(usize, u32)> {
+        self.0.first_nonzero()
+    }
+
+    pub fn iter_nonzero(&self) -> Vec<(usize, u32)> {
+        self.0.as_slice().iter_nonzero().collect()
+    }
+
+    pub fn slice(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSlice> {
+        checked_range(start, end, slf.0.len())?;
+        let py = slf.py();
+        Ok(PyFpSlice {
+            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
+            start,
+            end,
+        })
+    }
+
+    /// Restrict to the sub-range of coordinates `[start, end)`, returning a
+    /// read-only `FpSlice` view. For an `FpVector` (whose coordinates start
+    /// at 0) this mirrors `slice(start, end)`; named to match the analogous
+    /// `FpSlice.restrict` method.
+    pub fn restrict(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSlice> {
+        checked_range(start, end, slf.0.len())?;
+        let py = slf.py();
+        Ok(PyFpSlice {
+            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
+            start,
+            end,
+        })
+    }
+
+    pub fn slice_mut(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSliceMut> {
+        checked_range(start, end, slf.0.len())?;
+        let py = slf.py();
+        Ok(PyFpSliceMut {
+            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
+            start,
+            end,
+        })
+    }
+
+    /// A mutable `FpSliceMut` spanning the whole vector (mirrors upstream
+    /// `FpVector::as_slice_mut`); equivalent to `slice_mut(0, len())`.
+    pub fn as_slice_mut(slf: PyRef<'_, Self>) -> PyResult<PyFpSliceMut> {
+        let end = slf.0.len();
+        let py = slf.py();
+        Ok(PyFpSliceMut {
+            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
+            start: 0,
+            end,
+        })
+    }
+
+    pub fn set_entry(&mut self, index: usize, value: u32) -> PyResult<()> {
+        self.0.set_entry(checked_index(index, self.0.len())?, value);
+        Ok(())
+    }
+
+    pub fn scale(&mut self, c: u32) {
+        self.0.scale(c)
+    }
+
+    pub fn set_to_zero(&mut self) {
+        self.0.set_to_zero()
+    }
+
+    pub fn add_basis_element(&mut self, index: usize, value: u32) -> PyResult<()> {
+        self.0
+            .add_basis_element(checked_index(index, self.0.len())?, value);
+        Ok(())
+    }
+
+    pub fn extend_len(&mut self, len: usize) {
+        self.0.extend_len(len)
+    }
+
+    pub fn set_scratch_vector_size(&mut self, len: usize) {
+        self.0.set_scratch_vector_size(len)
+    }
+
+    pub fn to_bytes<'py>(&self, py: Python<'py>) -> PyResult<Bound<'py, PyBytes>> {
+        serialize_to_pybytes(py, |buffer| self.0.to_bytes(buffer))
+    }
+
+    pub fn update_from_bytes(&mut self, data: &[u8]) -> PyResult<()> {
+        self.0
+            .update_from_bytes(&mut Cursor::new(data))
+            .map_err(io_err)
+    }
+
+    pub fn __len__(&self) -> usize {
+        self.0.len()
+    }
+
+    pub fn __getitem__(&self, index: isize) -> PyResult<u32> {
+        Ok(self.0.entry(py_index(index, self.0.len())?))
+    }
+
+    pub fn __setitem__(&mut self, index: isize, value: u32) -> PyResult<()> {
+        self.0.set_entry(py_index(index, self.0.len())?, value);
+        Ok(())
+    }
+
+    pub fn __iter__(slf: PyRef<'_, Self>) -> PyFpVectorIterator {
+        PyFpVectorIterator {
+            entries: slf.0.iter().collect(),
+            index: 0,
+        }
+    }
+
+    /// Return an owned clone of this vector. Mirrors `FpSlice.to_owned`,
+    /// allowing `.to_owned()` to be used uniformly on both vectors and
+    /// slices.
+    pub fn to_owned(&self) -> PyFpVector {
+        Self(self.0.clone())
+    }
+
+    pub fn __repr__(&self) -> String {
+        format!("FpVector({}, {})", self.prime(), self.0)
+    }
+}
+
+#[pyclass(name = "FpSlice")]
+pub struct PyFpSlice {
+    pub(crate) parent: SliceParent,
+    pub(crate) start: usize,
+    pub(crate) end: usize,
+}
+
+impl PyFpSlice {
+    fn with_slice<R>(&self, py: Python<'_>, f: impl FnOnce(RustFpSlice<'_>) -> R) -> PyResult<R> {
+        with_parent_slice(&self.parent, self.start, self.end, py, f)
+    }
+
+    /// Cached span of the handle, used only for computing index bounds.
+    /// This does NOT revalidate the parent; callers that touch the parent
+    /// go through `with_slice`/`with_slice_mut`, which revalidate.
+    fn span(&self) -> usize {
+        self.end - self.start
+    }
+
+    fn to_owned_checked(&self, py: Python<'_>) -> PyResult<RustFpVector> {
+        self.with_slice(py, |s| s.to_owned())
+    }
+}
+
+#[pymethods]
+impl PyFpSlice {
+    #[getter]
+    pub fn prime(&self, py: Python<'_>) -> PyResult<u32> {
+        self.with_slice(py, |s| s.prime().as_u32())
+    }
+
+    #[getter]
+    pub fn len(&self, py: Python<'_>) -> PyResult<usize> {
+        self.with_slice(py, |s| s.len())
+    }
+
+    #[getter]
+    pub fn is_empty(&self, py: Python<'_>) -> PyResult<bool> {
+        self.with_slice(py, |s| s.is_empty())
+    }
+
+    pub fn entry(&self, py: Python<'_>, index: usize) -> PyResult<u32> {
+        let index = checked_index(index, self.span())?;
+        self.with_slice(py, |s| s.entry(index))
+    }
+
+    pub fn iter(&self, py: Python<'_>) -> PyResult<PyFpVectorIterator> {
+        let entries = self.with_slice(py, |s| s.iter().collect())?;
+        Ok(PyFpVectorIterator { entries, index: 0 })
+    }
+
+    pub fn iter_nonzero(&self, py: Python<'_>) -> PyResult<Vec<(usize, u32)>> {
+        self.with_slice(py, |s| s.iter_nonzero().collect())
+    }
+
+    #[getter]
+    pub fn is_zero(&self, py: Python<'_>) -> PyResult<bool> {
+        self.with_slice(py, |s| s.is_zero())
+    }
+
+    #[getter]
+    pub fn first_nonzero(&self, py: Python<'_>) -> PyResult<Option<(usize, u32)>> {
+        self.with_slice(py, |s| s.first_nonzero())
+    }
+
+    pub fn restrict(&self, py: Python<'_>, start: usize, end: usize) -> PyResult<Self> {
+        checked_range(start, end, self.span())?;
+        Ok(Self {
+            parent: self.parent.clone_ref(py),
+            start: self.start + start,
+            end: self.start + end,
+        })
+    }
+
+    pub fn to_owned(&self, py: Python<'_>) -> PyResult<PyFpVector> {
+        Ok(PyFpVector(self.to_owned_checked(py)?))
+    }
+
+    pub fn __len__(&self, py: Python<'_>) -> PyResult<usize> {
+        self.len(py)
+    }
+
+    pub fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult<u32> {
+        let index = py_index(index, self.span())?;
+        self.with_slice(py, |s| s.entry(index))
+    }
+
+    pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
+        self.with_slice(py, |s| format!("FpSlice({}, {})", s.prime().as_u32(), s))
+    }
+}
+
+#[pyclass(name = "FpSliceMut")]
+pub struct PyFpSliceMut {
+    pub(crate) parent: SliceParent,
+    pub(crate) start: usize,
+    pub(crate) end: usize,
+}
+
+impl PyFpSliceMut {
+    fn with_slice<R>(&self, py: Python<'_>, f: impl FnOnce(RustFpSlice<'_>) -> R) -> PyResult<R> {
+        with_parent_slice(&self.parent, self.start, self.end, py, f)
+    }
+
+    /// Cached span of the handle, used only for computing index bounds.
+    /// This does NOT revalidate the parent; callers that touch the parent
+    /// go through `with_slice`/`with_slice_mut`, which revalidate.
+    fn span(&self) -> usize {
+        self.end - self.start
+    }
+
+    fn with_slice_mut<R>(
+        &self,
+        py: Python<'_>,
+        f: impl FnOnce(RustFpSliceMut<'_>) -> R,
+    ) -> PyResult<R> {
+        with_parent_slice_mut(&self.parent, self.start, self.end, py, f)
+    }
+
+    /// Provide an immutable view of `operand` to `f`, avoiding a deep clone
+    /// in the common case. If `operand` is backed by a *different* Python
+    /// object from `self`, it is borrowed transiently (no clone) so that the
+    /// caller can hold this shared borrow alongside `self`'s mutable borrow.
+    /// Only when they share a backing object (genuine aliasing) does this
+    /// fall back to an owned copy, sidestepping the PyO3 double-borrow.
+    fn with_operand<R>(
+        &self,
+        py: Python<'_>,
+        operand: &PyFpSlice,
+        f: impl FnOnce(RustFpSlice<'_>) -> PyResult<R>,
+    ) -> PyResult<R> {
+        if self.parent.same_object(py, &operand.parent) {
+            let owned = operand.to_owned_checked(py)?;
+            f(owned.as_slice())
+        } else {
+            operand.with_slice(py, f)?
+        }
+    }
+
+    /// Like [`with_operand`], but accepts an operand that is either an
+    /// `FpSlice` or a (full) `FpVector`. An `FpSlice` is routed through
+    /// [`with_operand`] (preserving its clone-on-alias handling). An
+    /// `FpVector` is borrowed transiently and viewed via its full-vector
+    /// slice, falling back to an owned copy only when it shares its backing
+    /// Python object with the target (genuine aliasing), matching the
+    /// `FpSlice` path's behavior.
+    fn with_operand_any<R>(
+        &self,
+        py: Python<'_>,
+        operand: &Bound<'_, PyAny>,
+        f: impl FnOnce(RustFpSlice<'_>) -> PyResult<R>,
+    ) -> PyResult<R> {
+        if let Ok(slice) = operand.extract::<PyRef<'_, PyFpSlice>>() {
+            self.with_operand(py, &slice, f)
+        } else if let Ok(vector) = operand.cast::<PyFpVector>() {
+            if self.parent.same_vector(py, vector) {
+                let owned = vector.try_borrow().map_err(borrow_error)?.0.clone();
+                f(owned.as_slice())
+            } else {
+                let vector = vector.try_borrow().map_err(borrow_error)?;
+                f(vector.0.as_slice())
+            }
+        } else {
+            Err(PyValueError::new_err("expected an FpVector or FpSlice"))
+        }
+    }
+
+    /// Shared "sandwich" for single-operand binary ops: borrow `other`
+    /// (an `FpSlice` or full `FpVector`, via `with_operand_any`), then
+    /// borrow `self` mutably (via `with_slice_mut`), verify the primes
+    /// match, and run `f` with the mutable target and the operand slice.
+    /// Methods with extra pre-checks (length/offset/mask/span) keep those
+    /// in front of this call.
+    fn with_binary_op<R>(
+        &self,
+        py: Python<'_>,
+        other: &Bound<'_, PyAny>,
+        f: impl FnOnce(RustFpSliceMut<'_>, RustFpSlice<'_>) -> R,
+    ) -> PyResult<R> {
+        self.with_operand_any(py, other, |o| {
+            self.with_slice_mut(py, |t| {
+                checked_same_prime(t.prime().as_u32(), o.prime().as_u32())?;
+                Ok(f(t, o))
+            })?
+        })
+    }
+}
+
+#[pymethods]
+impl PyFpSliceMut {
+    #[getter]
+    pub fn prime(&self, py: Python<'_>) -> PyResult<u32> {
+        self.with_slice(py, |s| s.prime().as_u32())
+    }
+
+    #[getter]
+    pub fn len(&self, py: Python<'_>) -> PyResult<usize> {
+        self.with_slice(py, |s| s.len())
+    }
+
+    #[getter]
+    pub fn is_empty(&self, py: Python<'_>) -> PyResult<bool> {
+        self.with_slice(py, |s| s.is_empty())
+    }
+
+    pub fn set_entry(&self, py: Python<'_>, index: usize, value: u32) -> PyResult<()> {
+        let index = checked_index(index, self.span())?;
+        self.with_slice_mut(py, |mut s| s.set_entry(index, value))
+    }
+
+    pub fn set_to_zero(&self, py: Python<'_>) -> PyResult<()> {
+        self.with_slice_mut(py, |mut s| s.set_to_zero())
+    }
+
+    pub fn scale(&self, py: Python<'_>, c: u32) -> PyResult<()> {
+        self.with_slice_mut(py, |mut s| s.scale(c))
+    }
+
+    pub fn add(&self, py: Python<'_>, other: &Bound<'_, PyAny>, c: u32) -> PyResult<()> {
+        checked_equal_len(self.span(), slice_like_len(other)?)?;
+        self.with_binary_op(py, other, |mut target, other_slice| {
+            target.add(other_slice, c);
+        })
+    }
+
+    pub fn assign(&self, py: Python<'_>, other: &Bound<'_, PyAny>) -> PyResult<()> {
+        checked_equal_len(self.span(), slice_like_len(other)?)?;
+        self.with_binary_op(py, other, |mut target, other_slice| {
+            target.assign(other_slice);
+        })
+    }
+
+    pub fn add_tensor(
+        &self,
+        py: Python<'_>,
+        offset: usize,
+        coeff: u32,
+        left: &Bound<'_, PyAny>,
+        right: &Bound<'_, PyAny>,
+    ) -> PyResult<()> {
+        let width = slice_like_len(left)?
+            .checked_mul(slice_like_len(right)?)
+            .and_then(|width| offset.checked_add(width))
+            .ok_or_else(|| PyIndexError::new_err("tensor range overflows usize"))?;
+        checked_range(offset, width, self.span())?;
+        // Borrow each operand transiently, falling back to an owned copy
+        // only for one that shares a backing object with the target. Two
+        // shared borrows coexist fine; only the target's mutable borrow can
+        // collide with an operand that aliases it.
+        self.with_operand_any(py, left, |left_slice| {
+            self.with_operand_any(py, right, |right_slice| {
+                self.with_slice_mut(py, |mut target| {
+                    checked_same_prime(target.prime().as_u32(), left_slice.prime().as_u32())?;
+                    checked_same_prime(target.prime().as_u32(), right_slice.prime().as_u32())?;
+                    target.add_tensor(offset, coeff, left_slice, right_slice);
+                    Ok(())
+                })?
+            })
+        })
+    }
+
+    pub fn add_basis_element(&self, py: Python<'_>, index: usize, value: u32) -> PyResult<()> {
+        let index = checked_index(index, self.span())?;
+        self.with_slice_mut(py, |mut s| s.add_basis_element(index, value))
+    }
+
+    pub fn as_slice(&self, py: Python<'_>) -> PyFpSlice {
+        PyFpSlice {
+            parent: self.parent.clone_ref(py),
+            start: self.start,
+            end: self.end,
+        }
+    }
+
+    pub fn slice_mut(&self, py: Python<'_>, start: usize, end: usize) -> PyResult<Self> {
+        checked_range(start, end, self.span())?;
+        Ok(Self {
+            parent: self.parent.clone_ref(py),
+            start: self.start + start,
+            end: self.start + end,
+        })
+    }
+
+    pub fn to_owned(&self, py: Python<'_>) -> PyResult<PyFpVector> {
+        Ok(PyFpVector(self.with_slice(py, |s| s.to_owned())?))
+    }
+
+    pub fn __len__(&self, py: Python<'_>) -> PyResult<usize> {
+        self.len(py)
+    }
+
+    pub fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult<u32> {
+        let index = py_index(index, self.span())?;
+        self.with_slice(py, |s| s.entry(index))
+    }
+
+    pub fn __setitem__(&self, py: Python<'_>, index: isize, value: u32) -> PyResult<()> {
+        let index = py_index(index, self.span())?;
+        self.with_slice_mut(py, |mut s| s.set_entry(index, value))
+    }
+
+    pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
+        self.with_slice(py, |s| format!("FpSliceMut({}, {})", s.prime().as_u32(), s))
+    }
+}
+
+#[pyclass(name = "FpVectorIterator")]
+pub struct PyFpVectorIterator {
+    entries: Vec<u32>,
+    index: usize,
+}
+
+#[pymethods]
+impl PyFpVectorIterator {
+    pub fn __iter__(slf: PyRefMut<'_, Self>) -> PyRefMut<'_, Self> {
+        slf
+    }
+
+    pub fn __next__(&mut self) -> Option<u32> {
+        let value = self.entries.get(self.index).copied();
+        self.index += usize::from(value.is_some());
+        value
+    }
+}
+
+/// Run `f` on a borrowed immutable slice over a vector-like argument
+/// (`FpVector` or `FpSlice`), holding the shared borrow only for the
+/// duration of the call. This is the read-only sibling of
+/// [`with_target_slice_mut`]: it avoids the deep `FpVector` clone that
+/// [`extract_input_owned`] performs for every immutable input argument.
+///
+/// The transient borrow surfaces as a PyO3 borrow conflict (`RuntimeError`)
+/// if the same object is simultaneously borrowed mutably elsewhere — e.g.
+/// passed as both the input and the mutable target — rather than UB.
+///
+/// # Error taxonomy
+///
+/// We first dispatch on the object's *type* and only then attempt the
+/// borrow, so the two failure modes stay distinct:
+///  * a genuinely wrong type → `ValueError("expected an FpVector or
+///    FpSlice")`, and
+///  * a correct type that is already borrowed mutably elsewhere (aliasing)
+///    → the borrow conflict is propagated verbatim as `RuntimeError`.
+///
+/// Exposed `pub(crate)` so that other binding modules (e.g. `algebra_py`)
+/// reuse it for immutable input element arguments.
+pub(crate) fn with_input_slice<R>(
+    py: Python<'_>,
+    obj: &Bound<'_, PyAny>,
+    f: impl FnOnce(RustFpSlice<'_>) -> PyResult<R>,
+) -> PyResult<R> {
+    if let Ok(vector) = obj.cast::<PyFpVector>() {
+        let vector = vector.try_borrow().map_err(borrow_error)?;
+        f(vector.0.as_slice())
+    } else if let Ok(slice) = obj.cast::<PyFpSlice>() {
+        let slice = slice.try_borrow().map_err(borrow_error)?;
+        slice.with_slice(py, f)?
+    } else {
+        Err(PyValueError::new_err("expected an FpVector or FpSlice"))
+    }
+}
+
+/// Run `f` on the mutable slice backing a vector-like argument
+/// (`FpVector` or `FpSliceMut`), used as an output target.
+///
+/// # Error taxonomy
+///
+/// We dispatch on the object's *type* before attempting the mutable
+/// borrow, so the two failure modes stay distinct:
+///  * a genuinely wrong type → `ValueError("expected an FpVector or
+///    FpSliceMut")`, and
+///  * a correct type that is already borrowed elsewhere (e.g. the same
+///    `FpVector` simultaneously passed as a borrowed input via
+///    [`with_input_slice`] *and* as this mutable target) → the borrow
+///    conflict is propagated verbatim as `RuntimeError`.
+///
+/// Aliasing the mutable target with an input is therefore rejected with a
+/// `RuntimeError` (an intentional API change from the pre-clone-removal
+/// behavior, which silently succeeded by cloning the input first).
+///
+/// Exposed `pub(crate)` so that other binding modules (e.g. `algebra_py`)
+/// can accept a bound `fp_py` result argument for the `multiply_*` family;
+/// the closure receives the reconstructed `FpSliceMut` and may return a
+/// `PyResult` so callers can pre-validate (prime/length) inside the borrow.
+pub(crate) fn with_target_slice_mut<R>(
+    py: Python<'_>,
+    obj: &Bound<'_, PyAny>,
+    f: impl FnOnce(RustFpSliceMut<'_>) -> PyResult<R>,
+) -> PyResult<R> {
+    if let Ok(vector) = obj.cast::<PyFpVector>() {
+        let mut vector = vector.try_borrow_mut().map_err(borrow_error)?;
+        f(vector.0.as_slice_mut())
+    } else if let Ok(slice) = obj.cast::<PyFpSliceMut>() {
+        let slice = slice.try_borrow().map_err(borrow_error)?;
+        slice.with_slice_mut(py, f)?
+    } else {
+        Err(PyValueError::new_err("expected an FpVector or FpSliceMut"))
+    }
+}
diff --git a/ext_py/src/lib.rs b/ext_py/src/lib.rs
new file mode 100644
index 0000000000..995d910d5f
--- /dev/null
+++ b/ext_py/src/lib.rs
@@ -0,0 +1,12 @@
+use pyo3::prelude::*;
+
+mod fp_mod;
+
+pub use fp_mod::fp_py;
+
+#[pymodule]
+#[pyo3(name = "ext")]
+mod ext_py {
+    #[pymodule_export]
+    use super::fp_py;
+}
diff --git a/ext_py/tests/test_affine_subspace.py b/ext_py/tests/test_affine_subspace.py
new file mode 100644
index 0000000000..d1edf60270
--- /dev/null
+++ b/ext_py/tests/test_affine_subspace.py
@@ -0,0 +1,181 @@
+import pytest
+
+from ext import fp
+
+
+def make_subspace(p, rows):
+    s = fp.Subspace(p, len(rows[0]))
+    for row in rows:
+        s.add_vector(fp.FpVector.from_py(p, row))
+    return s
+
+
+def test_construction_and_queries():
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 0, 0])
+    aff = fp.AffineSubspace(offset, linear)
+
+    # Prime is passed/queried as a plain int.
+    assert aff.prime == 2
+    assert aff.ambient_dimension == 3
+    assert aff.dimension == 2
+    assert repr(aff) == "AffineSubspace([1, 0, 0] + {[0, 1, 0], [0, 0, 1]})"
+
+
+def test_offset_and_linear_part_roundtrip():
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 0, 0])
+    aff = fp.AffineSubspace(offset, linear)
+
+    # offset() returns an owned FpVector (the input reduced against the
+    # linear part); here it is unchanged because [1,0,0] is already reduced.
+    stored = aff.offset()
+    assert isinstance(stored, fp.FpVector)
+    assert list(stored) == [1, 0, 0]
+
+    # linear_part() returns an owned Subspace clone.
+    lp = aff.linear_part()
+    assert isinstance(lp, fp.Subspace)
+    assert lp.dimension == 2
+    assert lp.ambient_dimension == 3
+    assert [list(v) for v in lp.iter()] == [[0, 1, 0], [0, 0, 1]]
+
+
+def test_offset_is_reduced():
+    # Offset [1, 1, 0] reduces to [1, 0, 0] against span{[0,1,0],[0,0,1]}.
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 1, 0])
+    aff = fp.AffineSubspace(offset, linear)
+    assert list(aff.offset()) == [1, 0, 0]
+
+
+def test_offset_and_linear_part_are_owned_copies():
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 0, 0])
+    aff = fp.AffineSubspace(offset, linear)
+
+    # Mutating the returned offset must not alias the stored offset.
+    off = aff.offset()
+    off.set_entry(0, 0)
+    off.set_entry(1, 1)
+    assert list(aff.offset()) == [1, 0, 0]
+
+    # Mutating the returned linear part must not alias the stored linear part.
+    lin = aff.linear_part()
+    lin.set_to_zero()
+    assert aff.linear_part().dimension == 2
+    assert aff.dimension == 2
+    assert [list(v) for v in aff.linear_part().iter()] == [[0, 1, 0], [0, 0, 1]]
+
+
+def test_contains_membership():
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 0, 0])
+    aff = fp.AffineSubspace(offset, linear)
+
+    # [1,1,0] = offset + [0,1,0] is in the affine subspace.
+    assert aff.contains(fp.FpVector.from_py(2, [1, 1, 0]))
+    assert fp.FpVector.from_py(2, [1, 1, 0]) in aff
+
+    # [0,1,0] has the wrong first coordinate, so it is not contained.
+    assert not aff.contains(fp.FpVector.from_py(2, [0, 1, 0]))
+    assert fp.FpVector.from_py(2, [0, 1, 0]) not in aff
+
+
+def test_contains_accepts_slice():
+    linear = make_subspace(2, [[0, 1, 0], [0, 0, 1]])
+    offset = fp.FpVector.from_py(2, [1, 0, 0])
+    aff = fp.AffineSubspace(offset, linear)
+
+    v = fp.FpVector.from_py(2, [1, 1, 0])
+    assert aff.contains(v.slice(0, 3))
+
+
+def test_contains_space():
+    # a = origin + span{[0,1,0],[0,0,1]} (a linear subspace through 0).
+    a = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [0, 0, 0]),
+        make_subspace(2, [[0, 1, 0], [0, 0, 1]]),
+    )
+    # b = [0,1,0] + span{[0,1,0]}: both its linear part and offset lie in a.
+    b = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [0, 1, 0]),
+        make_subspace(2, [[0, 1, 0]]),
+    )
+    # c = [1,0,0] + span{[0,1,0]}: its offset is outside a.
+    c = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [1, 0, 0]),
+        make_subspace(2, [[0, 1, 0]]),
+    )
+
+    assert a.contains_space(b)
+    assert not a.contains_space(c)
+
+
+def test_sum_semantics():
+    # a = [1,0,0] + span{[0,1,0]}, b = [0,0,1] + span{[0,0,1]}.
+    a = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [1, 0, 0]),
+        make_subspace(2, [[0, 1, 0]]),
+    )
+    b = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [0, 0, 1]),
+        make_subspace(2, [[0, 0, 1]]),
+    )
+
+    # sum adds the already-stored, already-reduced offsets, not the
+    # constructor arguments. a.offset stays [1,0,0] (already reduced against
+    # span{[0,1,0]}); b's offset [0,0,1] reduces to [0,0,0] against
+    # span{[0,0,1]} when b is constructed. So the sum offset is
+    # a.offset [1,0,0] + b.offset [0,0,0] = [1,0,0], re-reduced against
+    # span{[0,1,0],[0,0,1]} (dim 2) -> [1,0,0]. No [1,0,1] intermediate occurs.
+    s = a.sum(b)
+    assert s.dimension == 2
+    assert list(s.offset()) == [1, 0, 0]
+    # [1,1,1] = offset + [0,1,1] is contained; [0,0,0] is not.
+    assert s.contains(fp.FpVector.from_py(2, [1, 1, 1]))
+    assert not s.contains(fp.FpVector.from_py(2, [0, 0, 0]))
+
+
+@pytest.mark.xfail(reason="panics")
+def test_invalid_construction_raises():
+    linear = make_subspace(2, [[0, 1, 0]])
+
+    # Offset length must match the linear part's ambient dimension.
+    with pytest.raises(ValueError):
+        fp.AffineSubspace(fp.FpVector.from_py(2, [1, 0]), linear)
+
+    # Offset prime must match the linear part's prime.
+    with pytest.raises(ValueError):
+        fp.AffineSubspace(fp.FpVector.from_py(3, [1, 0, 0]), linear)
+
+
+def test_invalid_inputs_raise():
+    aff = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [1, 0, 0]),
+        make_subspace(2, [[0, 1, 0]]),
+    )
+
+    # contains: wrong prime / dimension.
+    with pytest.raises(ValueError):
+        aff.contains(fp.FpVector.from_py(3, [1, 0, 0]))
+    with pytest.raises(ValueError):
+        aff.contains(fp.FpVector.from_py(2, [1, 0]))
+
+    # contains_space / sum: incompatible prime or ambient dimension.
+    other_prime = fp.AffineSubspace(
+        fp.FpVector.from_py(3, [1, 0, 0]),
+        make_subspace(3, [[0, 1, 0]]),
+    )
+    other_dim = fp.AffineSubspace(
+        fp.FpVector.from_py(2, [1, 0, 0, 0]),
+        make_subspace(2, [[0, 1, 0, 0]]),
+    )
+    with pytest.raises(ValueError):
+        aff.contains_space(other_prime)
+    with pytest.raises(ValueError):
+        aff.contains_space(other_dim)
+    with pytest.raises(ValueError):
+        aff.sum(other_prime)
+    with pytest.raises(ValueError):
+        aff.sum(other_dim)
diff --git a/ext_py/tests/test_augmented_matrix.py b/ext_py/tests/test_augmented_matrix.py
new file mode 100644
index 0000000000..4f1d4349aa
--- /dev/null
+++ b/ext_py/tests/test_augmented_matrix.py
@@ -0,0 +1,249 @@
+import pytest
+
+from ext import fp
+
+
+def test_import_both_classes():
+    assert hasattr(fp, "AugmentedMatrix2")
+    assert hasattr(fp, "AugmentedMatrix3")
+
+
+def test_construction_and_queries():
+    m = fp.AugmentedMatrix2(3, 2, [2, 2])
+    # Prime is passed/queried as a plain int.
+    assert m.prime == 3
+    assert m.rows == 2
+    assert m.segments == 2
+    assert m.is_zero
+    # The first segment starts at column 0; the matrix has as many columns as
+    # the final segment's end.
+    starts = m.segment_starts
+    ends = m.segment_ends
+    assert starts[0] == 0
+    assert ends[-1] == m.columns
+    assert len(starts) == 2
+    assert len(ends) == 2
+
+    m3 = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    assert m3.segments == 3
+    assert len(m3.segment_starts) == 3
+
+
+def test_invalid_segment_widths_raise():
+    # Wrong number of segment widths.
+    with pytest.raises(ValueError):
+        fp.AugmentedMatrix2(3, 2, [2])
+    with pytest.raises(ValueError):
+        fp.AugmentedMatrix2(3, 2, [2, 2, 2])
+    with pytest.raises(ValueError):
+        fp.AugmentedMatrix3(3, 2, [2, 2])
+    # Invalid prime.
+    with pytest.raises(ValueError):
+        fp.AugmentedMatrix2(4, 2, [2, 2])
+
+
+def test_add_identity_and_invalid_segments():
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    m.add_identity(1, 1)
+    start1 = m.segment_starts[1]
+    rows = m.to_py()
+    assert rows[0][start1] == 1
+    assert rows[1][start1 + 1] == 1
+
+    # Out-of-range segment indices raise IndexError.
+    with pytest.raises(IndexError):
+        m.add_identity(0, 2)
+    with pytest.raises(IndexError):
+        m.row_segment(0, 0, 2)
+
+    # Non-square segment (2 rows, width-3 segment) raises ValueError.
+    wide = fp.AugmentedMatrix2(2, 2, [3, 3])
+    with pytest.raises(ValueError):
+        wide.add_identity(0, 0)
+
+
+def test_row_segment_returns_owned_fpvector():
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    m.add_identity(1, 1)
+    seg = m.row_segment(0, 1, 1)
+    assert isinstance(seg, fp.FpVector)
+    assert len(seg) == 2
+    assert seg[0] == 1
+    assert seg[1] == 0
+
+
+def test_row_mut_writes_through_to_parent():
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    rm = m.row_mut(0)
+    assert isinstance(rm, fp.FpSliceMut)
+    assert len(rm) == m.columns
+    rm.set_entry(0, 1)
+    rm.slice_mut(2, 4).set_entry(1, 1)
+    assert m.to_py()[0][0] == 1
+    assert m.to_py()[0][3] == 1
+    with pytest.raises(IndexError):
+        m.row_mut(2)
+
+
+def test_into_matrix_returns_matrix_and_consumes():
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    m.add_identity(1, 1)
+    cols = m.columns
+    inner = m.into_matrix()
+    assert isinstance(inner, fp.Matrix)
+    assert inner.rows == 2
+    assert inner.columns == cols
+    # `into_matrix` now consumes the augmented matrix: any further use raises
+    # RuntimeError, and calling it again also raises.
+    with pytest.raises(RuntimeError):
+        m.rows
+    with pytest.raises(RuntimeError):
+        m.into_matrix()
+
+
+def test_compute_quasi_inverses_consumes():
+    m = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    m.add_identity(0, 0)
+    m.add_identity(1, 1)
+    m.add_identity(2, 2)
+    m.row_reduce()
+    a, b = m.compute_quasi_inverses()
+    assert isinstance(a, fp.QuasiInverse)
+    # The augmented matrix is consumed afterwards.
+    with pytest.raises(RuntimeError):
+        m.compute_quasi_inverses()
+    with pytest.raises(RuntimeError):
+        m.to_py()
+
+
+def test_augmented_matrix2_compute_image_and_quasi_inverse():
+    # [I | I] over F2: A is the 2x2 identity, the second block is the identity.
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    m.add_identity(0, 0)
+    m.add_identity(1, 1)
+    m.row_reduce()
+
+    image = m.compute_image()
+    assert isinstance(image, fp.Subspace)
+    assert image.prime == 2
+    assert image.dimension == 2
+
+    qi = m.compute_quasi_inverse()
+    assert isinstance(qi, fp.QuasiInverse)
+    assert qi.prime == 2
+    assert qi.source_dimension == 2
+    # A is the identity, so its quasi-inverse preimage is the identity too.
+    assert qi.preimage.to_py() == [[1, 0], [0, 1]]
+
+    kernel = m.compute_kernel()
+    assert isinstance(kernel, fp.Subspace)
+    # A is full rank, so the kernel is trivial.
+    assert kernel.dimension == 0
+
+
+def test_augmented_matrix3_compute_kernel_and_quasi_inverses():
+    # [A | B | I] with all square identity blocks over F3.
+    m = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    m.add_identity(0, 0)
+    m.add_identity(1, 1)
+    m.add_identity(2, 2)
+    m.row_reduce()
+
+    kernel = m.compute_kernel()
+    assert isinstance(kernel, fp.Subspace)
+    assert kernel.prime == 3
+
+    a, b = m.compute_quasi_inverses()
+    assert isinstance(a, fp.QuasiInverse)
+    assert isinstance(b, fp.QuasiInverse)
+    assert a.prime == 3
+    assert b.prime == 3
+    # A = I is full-rank 2->2, so its quasi-inverse maps F3^2 -> F3^2.
+    assert a.source_dimension == 2
+    assert a.target_dimension == 2
+    # The residual quasi-inverse inverts B (= I) and is itself a full-rank
+    # 2->2 map.
+    assert b.source_dimension == 2
+    assert b.target_dimension == 2
+
+
+@pytest.mark.xfail(reason="panics")
+def test_compute_methods_require_row_reduce():
+    # Reproduction from the review: calling compute_* on a freshly constructed
+    # (not row-reduced) augmented matrix must raise ValueError, not panic.
+    m2 = fp.AugmentedMatrix2(2, 2, [2, 2])
+    with pytest.raises(ValueError):
+        m2.compute_kernel()
+    with pytest.raises(ValueError):
+        m2.compute_image()
+    with pytest.raises(ValueError):
+        m2.compute_quasi_inverse()
+
+    m3 = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    with pytest.raises(ValueError):
+        m3.compute_kernel()
+    with pytest.raises(ValueError):
+        m3.compute_quasi_inverses()
+
+    # The failed (not-row-reduced) compute_* call must NOT consume the matrix;
+    # it is still usable afterwards (parity with the Rust test).
+    assert m2.rows == 2
+    m2.row_reduce()
+    assert isinstance(m2.compute_kernel(), fp.Subspace)
+    assert m3.rows == 2
+    m3.row_reduce()
+    a, _ = m3.compute_quasi_inverses()
+    assert isinstance(a, fp.QuasiInverse)
+
+
+def test_outstanding_handles_after_into_matrix_raise_not_panic():
+    # Cross-feature hazard: a `segment(...)` and a `row_segment_mut(...)` handle
+    # outstanding when the matrix is consumed via into_matrix() must raise a
+    # clean RuntimeError (consumed) when used, never panic.
+    m = fp.AugmentedMatrix2(3, 2, [2, 2])
+    m.add_identity(1, 1)
+    seg = m.segment(0, 1)
+    rseg = m.row_segment_mut(0, 0, 0)
+    m.into_matrix()
+    with pytest.raises(RuntimeError):
+        seg.rows
+    with pytest.raises(RuntimeError):
+        rseg.len
+    with pytest.raises(RuntimeError):
+        rseg.to_owned()
+
+
+def test_outstanding_handles_after_compute_quasi_inverses_raise_not_panic():
+    m = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    m.add_identity(0, 0)
+    m.add_identity(1, 1)
+    m.add_identity(2, 2)
+    seg = m.segment(0, 2)
+    rseg = m.row_segment_mut(1, 0, 0)
+    m.row_reduce()
+    m.compute_quasi_inverses()
+    with pytest.raises(RuntimeError):
+        seg.rows
+    with pytest.raises(RuntimeError):
+        rseg.len
+
+
+def test_into_matrix_repr_and_value_correctness():
+    # Build a small known matrix via add_identity in the identity segment.
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    m.add_identity(1, 1)
+    expected = m.to_py()
+    assert "consumed" not in repr(m)
+
+    inner = m.into_matrix()
+    # The returned matrix's entries equal the pre-consume contents.
+    assert inner.to_py() == expected
+    # __repr__ reflects the consumed state.
+    assert "consumed" in repr(m)
+
+
+def test_repr():
+    m = fp.AugmentedMatrix2(3, 2, [2, 2])
+    assert repr(m).startswith("AugmentedMatrix2(3, 2x")
+    m3 = fp.AugmentedMatrix3(2, 1, [1, 1, 1])
+    assert repr(m3).startswith("AugmentedMatrix3(2, 1x")
diff --git a/ext_py/tests/test_fp_vector_slices.py b/ext_py/tests/test_fp_vector_slices.py
new file mode 100644
index 0000000000..28ffd0ec1a
--- /dev/null
+++ b/ext_py/tests/test_fp_vector_slices.py
@@ -0,0 +1,329 @@
+import pytest
+
+from ext import fp
+
+
+def test_fp_vector_slice_queries_and_to_owned():
+    v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
+
+    s = v.slice(1, 4)
+
+    assert s.prime == 5
+    assert len(s) == 3
+    assert not s.is_empty
+    assert s.entry(1) == 2
+    assert s[1] == 2
+    assert s[-1] == 0
+    assert not s.is_zero
+    assert s.first_nonzero == (0, 1)
+    assert list(s.iter()) == [1, 2, 0]
+    assert s.iter_nonzero() == [(0, 1), (1, 2)]
+    assert repr(s) == "FpSlice(5, [1, 2, 0])"
+
+    restricted = s.restrict(1, 3)
+    assert len(restricted) == 2
+    assert [restricted[i] for i in range(len(restricted))] == [2, 0]
+
+    owned = s.to_owned()
+    assert isinstance(owned, fp.FpVector)
+    assert repr(owned) == "FpVector(5, [1, 2, 0])"
+
+
+def test_fp_vector_restrict_returns_subslice():
+    v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
+
+    r = v.restrict(1, 4)
+    assert isinstance(r, fp.FpSlice)
+    assert len(r) == 3
+    assert r.iter_nonzero() == [(0, 1), (1, 2)]
+    # Mirrors slice(start, end) for an FpVector.
+    assert r.iter_nonzero() == v.slice(1, 4).iter_nonzero()
+    assert repr(r.to_owned()) == "FpVector(5, [1, 2, 0])"
+
+    with pytest.raises(IndexError):
+        v.restrict(0, 6)
+    with pytest.raises(IndexError):
+        v.restrict(3, 2)
+
+
+def test_fp_vector_iter_nonzero():
+    v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
+    assert v.iter_nonzero() == [(1, 1), (2, 2), (4, 4)]
+    # Matches the full-vector slice.
+    assert v.iter_nonzero() == v.slice(0, len(v)).iter_nonzero()
+
+
+def test_fp_vector_to_owned_is_independent_clone():
+    v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
+    owned = v.to_owned()
+    assert isinstance(owned, fp.FpVector)
+    assert list(owned) == list(v)
+    # Mutating the source does not affect the clone.
+    v.set_entry(0, 3)
+    assert owned[0] == 0
+
+
+def test_fp_vector_slice_range_and_index_errors():
+    v = fp.FpVector(3, 4)
+    s = v.slice(1, 3)
+
+    with pytest.raises(IndexError):
+        v.slice(3, 2)
+    with pytest.raises(IndexError):
+        v.slice(0, 5)
+    with pytest.raises(IndexError):
+        v.slice_mut(0, 5)
+    with pytest.raises(IndexError):
+        s.entry(2)
+    with pytest.raises(IndexError):
+        s[2]
+    with pytest.raises(IndexError):
+        s[-3]
+    with pytest.raises(IndexError):
+        s.restrict(1, 3)
+
+
+def test_fp_slice_handles_parent_lifetime_after_original_reference_deleted():
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4])
+    s = v.slice(1, 3)
+    sm = v.slice_mut(2, 4)
+
+    del v
+
+    assert list(s.iter()) == [2, 3]
+    sm.add_basis_element(0, 4)
+    assert s[1] == 2
+    assert list(sm.as_slice().iter()) == [2, 4]
+
+
+def test_stale_slice_handles_raise_python_exception_after_parent_shrink():
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4])
+    s = v.slice(1, 4)
+    sm = v.slice_mut(1, 4)
+
+    v.set_scratch_vector_size(2)
+
+    with pytest.raises(IndexError):
+        s.prime
+    with pytest.raises(IndexError):
+        s[0]
+    with pytest.raises(IndexError):
+        repr(s)
+    with pytest.raises(IndexError):
+        sm[0]
+    with pytest.raises(IndexError):
+        sm.set_entry(0, 1)
+    with pytest.raises(IndexError):
+        sm.set_to_zero()
+
+
+def test_fp_slice_mut_updates_parent_and_as_slice():
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0])
+    s = v.slice_mut(1, 4)
+
+    assert s.prime == 5
+    assert len(s) == 3
+    assert s[0] == 2
+    assert s[-1] == 4
+    assert repr(s) == "FpSliceMut(5, [2, 3, 4])"
+
+    s.set_entry(0, 7)
+    assert v[1] == 2
+
+    s[1] = 9
+    assert v[2] == 4
+
+    s.add_basis_element(2, 3)
+    assert v[3] == 2
+
+    s.scale(2)
+    assert [v[i] for i in range(len(v))] == [1, 4, 3, 4, 0]
+
+    as_slice = s.as_slice()
+    assert len(as_slice) == 3
+    assert [as_slice[i] for i in range(len(as_slice))] == [4, 3, 4]
+
+    sub = s.slice_mut(1, 3)
+    sub.set_to_zero()
+    assert [v[i] for i in range(len(v))] == [1, 4, 0, 0, 0]
+
+
+def test_fp_slice_mut_add_assign_and_masked_methods():
+    v = fp.FpVector.from_py(5, [1, 0, 0, 0, 0, 0])
+    target = v.slice_mut(1, 4)
+    source_v = fp.FpVector.from_py(5, [1, 2, 3, 4])
+    source = source_v.slice(0, 3)
+
+    target.add(source, 2)
+    assert [v[i] for i in range(len(v))] == [1, 2, 4, 1, 0, 0]
+
+    target.assign(fp.FpVector.from_py(5, [4, 0, 1]).slice(0, 3))
+    assert [v[i] for i in range(len(v))] == [1, 4, 0, 1, 0, 0]
+
+
+def test_fp_slice_mut_add_tensor():
+    v = fp.FpVector(5, 6)
+
+    v.slice_mut(0, 6).add_tensor(
+        1,
+        2,
+        fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
+        fp.FpVector.from_py(5, [3, 4]).slice(0, 2),
+    )
+
+    assert [v[i] for i in range(len(v))] == [0, 1, 3, 2, 1, 0]
+
+def test_fp_slice_mut_assign_distinct_and_aliased():
+    vt = fp.FpVector.from_py(5, [1, 2, 3])
+    vo = fp.FpVector.from_py(5, [4, 0, 1])
+    vt.slice_mut(0, 3).assign(vo.slice(0, 3))
+    assert [vt[i] for i in range(3)] == [4, 0, 1]
+
+    # Aliased clone fallback over the same vector (non-overlapping halves).
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0, 1])
+    v.slice_mut(0, 3).assign(v.slice(3, 6))
+    assert [v[i] for i in range(6)] == [4, 0, 1, 4, 0, 1]
+
+
+def test_fp_slice_mut_add_tensor_distinct_and_aliased():
+    # Distinct operands, distinct target (no-clone path).
+    vt = fp.FpVector(5, 6)
+    left = fp.FpVector.from_py(5, [1, 2])
+    right = fp.FpVector.from_py(5, [3, 4])
+    vt.slice_mut(0, 6).add_tensor(1, 2, left.slice(0, 2), right.slice(0, 2))
+    assert [vt[i] for i in range(6)] == [0, 1, 3, 2, 1, 0]
+
+    # Dual-operand left == right over the SAME vector (two shared borrows of
+    # one operand coexist; neither aliases the distinct target).
+    vt2 = fp.FpVector(5, 6)
+    w = fp.FpVector.from_py(5, [1, 2])
+    vt2.slice_mut(0, 6).add_tensor(1, 2, w.slice(0, 2), w.slice(0, 2))
+    assert [vt2[i] for i in range(6)] == [0, 2, 4, 4, 3, 0]
+
+    # An operand aliasing the target: the clone fallback must read the operand's
+    # ORIGINAL contents even though the target overwrites them mid-op.
+    v = fp.FpVector.from_py(5, [1, 2, 0, 0, 0, 0])
+    r = fp.FpVector.from_py(5, [3, 4])
+    v.slice_mut(0, 6).add_tensor(1, 2, v.slice(0, 2), r.slice(0, 2))
+    assert [v[i] for i in range(6)] == [1, 3, 3, 2, 1, 0]
+
+
+def test_fp_slice_mut_matrix_row_aliasing_clone_fallback():
+    # Two slices over different rows of the SAME matrix are conservatively
+    # treated as aliased (the whole Matrix pyclass is borrowed as a unit), so
+    # this takes the clone fallback. The result must still be correct.
+    m = fp.Matrix.from_py(5, [[1, 1, 1], [2, 2, 2]])
+    m.row_mut(0).add(m.row(1), 1)
+    assert m.to_py() == [[3, 3, 3], [2, 2, 2]]
+
+    # assign across rows of one matrix likewise clones but stays correct.
+    m2 = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1]])
+    m2.row_mut(0).assign(m2.row(1))
+    assert m2.to_py() == [[4, 0, 1], [4, 0, 1]]
+
+
+def test_fp_slice_mut_accepts_fp_vector_operand():
+    # The operand-taking FpSliceMut methods accept either an FpSlice or a
+    # (full) FpVector, coercing the FpVector via its full-vector slice.
+    v = fp.FpVector.from_py(5, [1, 0, 0, 0, 0, 0])
+    target = v.slice_mut(1, 4)
+    operand = fp.FpVector.from_py(5, [1, 2, 3])
+
+    target.add(operand, 2)
+    assert [v[i] for i in range(len(v))] == [1, 2, 4, 1, 0, 0]
+
+    target.assign(operand)
+    assert [v[i] for i in range(len(v))] == [1, 1, 2, 3, 0, 0]
+
+    # An FpVector operand whose backing object IS the target's parent takes the
+    # clone fallback and must still be correct.
+    w = fp.FpVector.from_py(5, [1, 2, 3])
+    w.slice_mut(0, 3).add(w, 1)
+    assert [w[i] for i in range(len(w))] == [2, 4, 1]
+
+    # add_tensor accepts FpVector operands too.
+    out = fp.FpVector(5, 6)
+    left = fp.FpVector.from_py(5, [1, 2])
+    right = fp.FpVector.from_py(5, [3, 4])
+    out.slice_mut(0, 6).add_tensor(1, 2, left, right)
+    assert [out[i] for i in range(6)] == [0, 1, 3, 2, 1, 0]
+
+    # A matrix row_mut accepting an FpVector operand (the secondary_massey.py
+    # pattern: out.slice_mut(...).add(g, 1) where g is an FpVector).
+    m = fp.Matrix.from_py(5, [[1, 1, 1], [2, 2, 2]])
+    m.row_mut(0).add(fp.FpVector.from_py(5, [3, 0, 1]), 1)
+    assert m.to_py() == [[4, 1, 2], [2, 2, 2]]
+
+
+def test_fp_slice_mut_new_method_errors_are_python_exceptions():
+    target_v = fp.FpVector(5, 3)
+    target = target_v.slice_mut(0, 3)
+    same_prime_short = fp.FpVector(5, 2).slice(0, 2)
+    other_prime = fp.FpVector(7, 3).slice(0, 3)
+
+    with pytest.raises(ValueError):
+        target.add(same_prime_short, 1)
+    with pytest.raises(ValueError):
+        target.add(other_prime, 1)
+    with pytest.raises(IndexError):
+        target.add_tensor(
+            1,
+            1,
+            fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
+            fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
+        )
+
+
+def test_fp_slice_mut_to_owned():
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0])
+    s = v.slice_mut(1, 4)
+
+    owned = s.to_owned()
+    assert isinstance(owned, fp.FpVector)
+    assert owned.prime == 5
+    assert repr(owned) == "FpVector(5, [2, 3, 4])"
+
+    # to_owned is a copy: mutating the source does not change the owned vector.
+    s.set_entry(0, 0)
+    assert repr(owned) == "FpVector(5, [2, 3, 4])"
+
+
+def test_fp_slice_len_revalidates_after_parent_shrink():
+    v = fp.FpVector.from_py(5, [1, 2, 3, 4])
+    s = v.slice(1, 4)
+    sm = v.slice_mut(1, 4)
+
+    assert len(s) == 3
+    assert len(sm) == 3
+    assert not s.is_empty
+
+    v.set_scratch_vector_size(2)
+
+    with pytest.raises(IndexError):
+        len(s)
+    with pytest.raises(IndexError):
+        s.is_empty
+    with pytest.raises(IndexError):
+        len(sm)
+    with pytest.raises(IndexError):
+        sm.is_empty
+
+
+def test_fp_slice_mut_index_and_range_errors():
+    v = fp.FpVector(2, 3)
+    s = v.slice_mut(1, 3)
+
+    with pytest.raises(IndexError):
+        s.set_entry(2, 1)
+    with pytest.raises(IndexError):
+        s.add_basis_element(2, 1)
+    with pytest.raises(IndexError):
+        s[2]
+    with pytest.raises(IndexError):
+        s[-3]
+    with pytest.raises(IndexError):
+        s[2] = 1
+    with pytest.raises(IndexError):
+        s.slice_mut(2, 1)
+    with pytest.raises(IndexError):
+        s.slice_mut(0, 3)
diff --git a/ext_py/tests/test_matrix.py b/ext_py/tests/test_matrix.py
new file mode 100644
index 0000000000..46a83a9eed
--- /dev/null
+++ b/ext_py/tests/test_matrix.py
@@ -0,0 +1,386 @@
+import pytest
+
+from ext import fp
+
+
+def test_matrix_construction_and_queries():
+    m = fp.Matrix(7, 2, 3)
+    assert m.prime == 7
+    assert m.rows == 2
+    assert m.columns == 3
+    assert m.is_zero
+    assert len(m) == 2
+    assert m.to_py() == [[0, 0, 0], [0, 0, 0]]
+    assert repr(m).startswith("Matrix(7, ")
+
+
+@pytest.mark.xfail(reason="Needs fix")
+def test_matrix_to_py_zero_dimensions():
+    # Regression: a matrix with zero columns and nonzero rows shouldn't panic
+    assert fp.Matrix(2, 1, 0).to_py() == [[]]
+    assert fp.Matrix(2, 3, 0).to_py() == [[], [], []]
+    # Zero rows (with or without columns) is the empty list of rows.
+    assert fp.Matrix(2, 0, 3).to_py() == []
+    assert fp.Matrix(2, 0, 0).to_py() == []
+    # Sibling row-materializing access must also not panic on zero columns.
+    assert fp.Matrix(2, 1, 0).rows == 1
+    assert fp.Matrix(2, 1, 0).columns == 0
+    assert list(fp.Matrix(2, 1, 0).row(0)) == []
+
+
+def test_matrix_from_py_and_identity():
+    m = fp.Matrix.from_py(7, [[1, 3, 6], [0, 3, 4]])
+    assert m.to_py() == [[1, 3, 6], [0, 3, 4]]
+    assert not m.is_zero
+
+    ident = fp.Matrix.identity(5, 3)
+    assert ident.to_py() == [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
+
+
+def test_matrix_from_rows_and_from_row():
+    r0 = fp.FpVector.from_py(5, [1, 2, 3])
+    r1 = fp.FpVector.from_py(5, [4, 0, 1])
+    m = fp.Matrix.from_rows(5, [r0, r1], 3)
+    assert m.to_py() == [[1, 2, 3], [4, 0, 1]]
+
+    single = fp.Matrix.from_row(5, r0, 3)
+    assert single.to_py() == [[1, 2, 3]]
+
+
+def test_matrix_augmented_from_py():
+    py_m = [[1, 3, 6], [0, 3, 4]]
+    m = fp.AugmentedMatrix2.from_py(7, py_m)
+    assert m.rows == 2
+    assert m.columns == 12
+    assert m.segment_starts == [0, 10]
+    assert m.segment_ends == [3, 12]
+    assert m.segment(0, 0).to_py() == py_m
+    assert m.segment(1, 1).to_py() == [[1, 0], [0, 1]]
+
+
+def test_prime_is_int():
+    m = fp.Matrix(5, 1, 1)
+    assert isinstance(m.prime, int)
+
+
+def test_invalid_prime_and_dims():
+    with pytest.raises(ValueError):
+        fp.Matrix(1, 2, 2)
+    with pytest.raises(ValueError):
+        fp.Matrix.from_py(4, [[1, 2]])
+    with pytest.raises(ValueError):
+        fp.Matrix.from_py(7, [[1, 2], [3]])
+
+
+def test_row_access_and_getitem():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1]])
+    row = m.row(1)
+    assert row.prime == 5
+    assert len(row) == 3
+    assert row.entry(0) == 4
+    assert row[2] == 1
+    assert row[-1] == 1
+    assert not row.is_zero
+    assert row.first_nonzero == (0, 4)
+    assert list(row.iter()) == [4, 0, 1]
+    assert row.iter_nonzero() == [(0, 4), (2, 1)]
+    assert m[0].to_owned().prime == 5
+    assert list(m[0].iter()) == [1, 2, 3]
+
+    with pytest.raises(IndexError):
+        m.row(2)
+    with pytest.raises(IndexError):
+        row.entry(3)
+
+
+def test_row_mut_reflects_in_parent():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1]])
+    rm = m.row_mut(0)
+    rm.set_entry(0, 9)
+    assert m.to_py()[0] == [4, 2, 3]
+    rm[1] = 3
+    assert m.row(0)[1] == 3
+    rm.scale(2)
+    assert m.to_py()[0] == [3, 1, 1]
+    rm.set_to_zero()
+    assert m.to_py()[0] == [0, 0, 0]
+    rm.add_basis_element(2, 1)
+    assert m.to_py()[0] == [0, 0, 1]
+
+    with pytest.raises(IndexError):
+        rm.set_entry(3, 1)
+
+
+def test_row_mut_add_slice():
+    m = fp.Matrix.from_py(5, [[1, 2, 3]])
+    other = fp.FpVector.from_py(5, [1, 1, 1])
+    m.row_mut(0).add(other.slice(0, 3), 2)
+    assert m.to_py()[0] == [3, 4, 0]
+
+
+def test_iter_mut_writes_through_and_counts_rows():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1]])
+    rows = list(m.iter_mut())
+    assert len(rows) == m.rows == 2
+    assert type(rows[0]) is type(m.row_mut(0))
+    for i, out in enumerate(rows):
+        out.set_entry(i, 3)
+    assert m.to_py() == [[3, 2, 3], [4, 3, 1]]
+    # A zero-row matrix yields an empty iterator.
+    assert list(fp.Matrix(5, 0, 3).iter_mut()) == []
+
+
+def test_mutators():
+    m = fp.Matrix.from_py(5, [[1, 2], [3, 4]])
+    m.swap_rows(0, 1)
+    assert m.to_py() == [[3, 4], [1, 2]]
+    with pytest.raises(IndexError):
+        m.swap_rows(0, 2)
+
+    m.row_op(0, 1, 1)
+    assert m.to_py() == [[4, 1], [1, 2]]
+    with pytest.raises(ValueError):
+        m.row_op(0, 0, 1)
+
+    m.set_to_zero()
+    assert m.is_zero
+
+
+def test_add_row_extends_matrix():
+    m = fp.Matrix(5, 1, 2)
+    new = m.add_row()
+    assert m.rows == 2
+    new.set_entry(0, 3)
+    assert m.to_py()[1] == [3, 0]
+
+
+def test_assign_requires_matching_shape():
+    m = fp.Matrix.from_py(5, [[1, 2], [3, 4]])
+    other = fp.Matrix.from_py(5, [[1, 1], [1, 1]])
+    m.assign(other)
+    assert m.to_py() == [[1, 1], [1, 1]]
+
+    mismatch = fp.Matrix(5, 3, 2)
+    with pytest.raises(ValueError):
+        m.assign(mismatch)
+    diff_prime = fp.Matrix(7, 2, 2)
+    with pytest.raises(ValueError):
+        m.assign(diff_prime)
+
+
+def test_row_reduce_rank():
+    m = fp.Matrix.from_py(2, [[1, 1, 0], [0, 1, 1], [1, 0, 1]])
+    assert m.row_reduce() == 2
+
+    pivots = m.pivots
+    assert isinstance(pivots, list)
+
+
+def test_extend_columns_and_pivots():
+    m = fp.Matrix(5, 2, 2)
+    m.extend_column_dimension(4)
+    assert m.columns == 4
+    m.initialize_pivots()
+    assert m.pivots == [-1, -1, -1, -1]
+
+
+def test_trim_and_rotate():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1], [2, 2, 2]])
+    m.trim(0, 2, 1)
+    assert m.to_py() == [[2, 3], [0, 1]]
+
+    n = fp.Matrix.from_py(5, [[1, 0], [2, 0], [3, 0]])
+    n.rotate_down(0, 3, 1)
+    assert n.to_py() == [[3, 0], [1, 0], [2, 0]]
+
+
+def test_bytes_roundtrip():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 2]])
+    data = m.to_bytes()
+    n = fp.Matrix.from_bytes(5, 2, 3, data)
+    assert n.to_py() == m.to_py()
+
+
+def test_stale_row_handle_after_trim_raises():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1], [2, 2, 2]])
+    row = m.row(2)
+    assert row.entry(0) == 2
+    m.trim(0, 1, 0)
+    assert m.rows == 1
+    with pytest.raises(IndexError):
+        row.entry(0)
+    with pytest.raises(IndexError):
+        row[0]
+
+
+def test_stale_row_mut_handle_after_trim_raises():
+    m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1], [2, 2, 2]])
+    rm = m.row_mut(2)
+    rm.set_entry(0, 1)
+    m.trim(0, 1, 0)
+    with pytest.raises(IndexError):
+        rm.set_entry(0, 1)
+
+
+def test_row_returns_same_type_as_vector_slice():
+    m = fp.Matrix.from_py(5, [[1, 2, 3]])
+    v = fp.FpVector.from_py(5, [1, 2, 3])
+    assert type(m.row(0)) is type(v.slice(0, 3))
+    assert type(m.row_mut(0)) is type(v.slice_mut(0, 3))
+    assert type(m[0]) is type(v.slice(0, 3))
+
+
+def test_row_slice_restrict_and_to_owned():
+    m = fp.Matrix.from_py(5, [[1, 2, 3, 4]])
+    row = m.row(0)
+    sub = row.restrict(1, 3)
+    assert len(sub) == 2
+    assert [sub[i] for i in range(len(sub))] == [2, 3]
+    owned = row.to_owned()
+    assert isinstance(owned, fp.FpVector)
+    assert repr(owned) == "FpVector(5, [1, 2, 3, 4])"
+    assert repr(row).startswith("FpSlice(5, ")
+
+
+def test_row_mut_to_owned_and_slice_mut():
+    m = fp.Matrix.from_py(5, [[1, 2, 3, 4]])
+    rm = m.row_mut(0)
+    owned = rm.to_owned()
+    assert owned.prime == 5
+    sub = rm.slice_mut(0, 2)
+    sub.scale(2)
+    assert m.to_py()[0] == [2, 4, 3, 4]
+    assert repr(rm).startswith("FpSliceMut(5, ")
+
+
+def test_row_len_revalidates_after_column_shrink():
+    # `trim` with col_start > 0 reduces the number of columns, so a previously
+    # created row slice whose `end` exceeds the new column count is stale.
+    m = fp.Matrix.from_py(5, [[1, 2, 3, 4], [5, 6, 7, 8]])
+    row = m.row(0)
+    rm = m.row_mut(1)
+    assert len(row) == 4
+    assert len(rm) == 4
+
+    # Drop two leading columns: columns goes 4 -> 2.
+    m.trim(0, 2, 2)
+    assert m.columns == 2
+
+    with pytest.raises(IndexError):
+        len(row)
+    with pytest.raises(IndexError):
+        row.entry(0)
+    with pytest.raises(IndexError):
+        row[0]
+    with pytest.raises(IndexError):
+        len(rm)
+    with pytest.raises(IndexError):
+        rm.set_entry(0, 1)
+
+
+def test_content_shift_staleness_after_trim_is_documented_behavior():
+    # A sub-slice that survives a `col_start > 0` trim (its `end` stays within
+    # the new column count) does NOT raise: revalidation guards parent
+    # dimensions only, not logical-coordinate remapping. The slice now reads the
+    # shifted columns. This pins the documented semantics.
+    m = fp.Matrix.from_py(5, [[1, 2, 3, 4]])
+    # restrict to absolute indices 1..3, i.e. logical columns [2, 3].
+    sub = m.row(0).restrict(1, 3)
+    assert [sub[i] for i in range(len(sub))] == [2, 3]
+
+    # Trim one leading column: columns goes 4 -> 3, data shifts left by one.
+    # Row was [1, 2, 3, 4] -> [2, 3, 4]. The sub-slice's range 1..3 still fits.
+    m.trim(0, 1, 1)
+    assert m.columns == 3
+    assert m.to_py()[0] == [2, 3, 4]
+
+    # The handle does not raise; it now reads the remapped indices 1..3 -> [3, 4]
+    # instead of the original [2, 3]. Length is unchanged but data has shifted.
+    assert len(sub) == 2
+    assert [sub[i] for i in range(len(sub))] == [3, 4]
+
+
+def test_interacting_mutable_slices_over_same_parent():
+    # Borrow conflicts cannot arise from holding two slice handles: each call
+    # reconstructs and re-borrows the parent for the duration of that single
+    # call, and no parent borrow is held across re-entry into Python. So two
+    # FpSliceMut handles over the same matrix interleave safely rather than
+    # panicking. We assert the safe interleaving here; a genuine borrow conflict
+    # is not reachable from Python with this design.
+    m = fp.Matrix.from_py(5, [[1, 2], [3, 4]])
+    a = m.row_mut(0)
+    b = m.row_mut(1)
+
+    a.set_entry(0, 0)
+    b.set_entry(1, 0)
+    a.scale(2)
+
+    assert m.to_py() == [[0, 4], [3, 0]]
+
+
+# --- compute_kernel / compute_image ----------------------------------------
+
+# Example from the Rust `compute_kernel`/`compute_image` docs at p = 3.
+_KI_ROWS = [
+    [1, 2, 1, 1, 0],
+    [1, 0, 2, 1, 1],
+    [2, 2, 0, 2, 1],
+]
+
+
+def test_compute_kernel_dimensions_and_membership():
+    m = fp.AugmentedMatrix2.from_py(3, _KI_ROWS)
+    m.row_reduce()
+    ker = m.compute_kernel()
+    assert ker.prime == 3
+    # The kernel here is the left null space (combinations of the 3 rows that
+    # vanish), so it is ambient-dimension 3 and one-dimensional.
+    assert ker.ambient_dimension == len(_KI_ROWS)
+    assert ker.dimension == 1
+    basis = ker.basis
+    assert len(basis) == 1
+    # Kernel basis row from the upstream doc example.
+    assert list(basis[0]) == [1, 1, 2]
+    # Verify each kernel row r satisfies sum_i r[i] * row_i == 0.
+    orig = _KI_ROWS
+    ncols = len(orig[0])
+    for row in basis:
+        product = fp.FpVector.new(3, ncols)
+        for i in range(len(orig)):
+            c = row.entry(i)
+            if c:
+                for j in range(ncols):
+                    product.set_entry(j, (product.entry(j) + c * orig[i][j]) % 3)
+        assert product.is_zero
+
+
+def test_compute_image_dimensions_and_membership():
+    m = fp.AugmentedMatrix2.from_py(3, _KI_ROWS)
+    m.row_reduce()
+    img = m.compute_image()
+    assert img.prime == 3
+    assert img.ambient_dimension == len(_KI_ROWS[0])
+    # Image of a rank-3 map into a 5-dim target has dimension 3... but the
+    # doc example restricts to the target block; the image basis matches the
+    # upstream doc rows.
+    rows = [list(r) for r in img.basis]
+    assert [1, 0, 2, 1, 1] in rows
+    assert [0, 1, 1, 0, 1] in rows
+    # The recorded image rows are members of the image subspace.
+    for r in img.basis:
+        assert img.contains(r)
+
+
+
+@pytest.mark.xfail(reason="Raises PanicException")
+def test_compute_kernel_requires_row_reduce():
+    m = fp.AugmentedMatrix2.from_py(3, [[1, 0, 1], [0, 1, 1]])
+    with pytest.raises(ValueError):
+        m.compute_kernel()
+
+
+@pytest.mark.xfail(reason="Raises PanicException")
+def test_compute_image_requires_row_reduce():
+    m = fp.AugmentedMatrix2.from_py(3, [[1, 0, 1], [0, 1, 1]])
+    with pytest.raises(ValueError):
+        m.compute_image()
diff --git a/ext_py/tests/test_matrix_slice_mut.py b/ext_py/tests/test_matrix_slice_mut.py
new file mode 100644
index 0000000000..55185da1b3
--- /dev/null
+++ b/ext_py/tests/test_matrix_slice_mut.py
@@ -0,0 +1,196 @@
+import pytest
+
+from ext import fp
+
+
+def make_matrix(rows):
+    """Build an F5 matrix from a list-of-lists via row_mut.set_entry."""
+    p = 5
+    m = fp.Matrix(p, len(rows), len(rows[0]))
+    for i, row in enumerate(rows):
+        rm = m.row_mut(i)
+        for j, v in enumerate(row):
+            rm.set_entry(j, v)
+    return m
+
+
+def test_slice_mut_construction_and_queries():
+    m = make_matrix([[1, 2, 3, 4], [0, 1, 2, 3], [4, 3, 2, 1]])
+    rect = m.slice_mut(0, 2, 1, 3)
+    # prime returned as a plain int.
+    assert rect.prime == 5
+    assert rect.rows == 2
+    assert rect.columns == 2
+    assert repr(rect) == "MatrixSliceMut(5, 2x2)"
+
+
+def test_row_and_row_slice_read():
+    m = make_matrix([[1, 2, 3, 4], [0, 1, 2, 3], [4, 3, 2, 1]])
+    rect = m.slice_mut(0, 3, 1, 3)
+    # row(i) is the column-shifted view into the rectangle.
+    row0 = rect.row(0)
+    assert isinstance(row0, fp.FpSlice)
+    assert len(row0) == 2
+    assert row0[0] == 2
+    assert row0[1] == 3
+
+    # row_slice restricts the row range, keeping the columns.
+    sub = rect.row_slice(1, 3)
+    assert sub.rows == 2
+    assert sub.columns == 2
+    assert sub.row(0)[0] == 1  # original row 1, column 1
+
+
+def test_row_mut_writes_through_to_parent():
+    m = make_matrix([[1, 2, 3, 4], [0, 1, 2, 3]])
+    rect = m.slice_mut(0, 2, 1, 4)
+    rm = rect.row_mut(0)
+    assert isinstance(rm, fp.FpSliceMut)
+    rm.set_entry(0, 0)  # rectangle column 0 == matrix column 1
+    assert m.to_py()[0] == [1, 0, 3, 4]
+
+
+def test_iter_mut_reflects_in_parent():
+    m = make_matrix([[1, 1, 1, 1], [2, 2, 2, 2]])
+    rect = m.slice_mut(0, 2, 0, 4)
+    # iter yields read-only row handles.
+    for r in rect.iter():
+        assert isinstance(r, fp.FpSlice)
+    # iter_mut yields mutable handles that write through.
+    for r in rect.iter_mut():
+        r.set_entry(0, 0)
+    assert m.to_py() == [[0, 1, 1, 1], [0, 2, 2, 2]]
+
+
+def test_add_identity():
+    m = fp.Matrix(3, 2, 4)
+    rect = m.slice_mut(0, 2, 2, 4)
+    rect.add_identity()
+    assert m.to_py() == [[0, 0, 1, 0], [0, 0, 0, 1]]
+
+    # Non-square rectangle raises ValueError.
+    wide = m.slice_mut(0, 2, 0, 4)
+    with pytest.raises(ValueError):
+        wide.add_identity()
+
+
+def test_invalid_rectangle_raises():
+    m = make_matrix([[1, 2, 3], [4, 0, 1]])
+    with pytest.raises(IndexError):
+        m.slice_mut(0, 5, 0, 1)  # too many rows
+    with pytest.raises(IndexError):
+        m.slice_mut(0, 1, 0, 9)  # too many columns
+    with pytest.raises(IndexError):
+        m.slice_mut(2, 1, 0, 1)  # inverted row range
+
+
+def test_out_of_range_row_index_raises():
+    m = make_matrix([[1, 2, 3], [4, 0, 1]])
+    rect = m.slice_mut(0, 2, 0, 3)
+    with pytest.raises(IndexError):
+        rect.row(2)
+    with pytest.raises(IndexError):
+        rect.row_mut(2)
+
+
+def test_stale_handle_after_parent_shrinks_raises():
+    m = make_matrix([[1, 1, 0], [0, 1, 1]])
+    rect = m.slice_mut(0, 2, 0, 3)
+    # A square sub-rectangle, valid before the parent shrinks.
+    square = m.slice_mut(0, 2, 0, 2)
+    # Trim the parent to a single row; the 2-row rectangles are now stale.
+    m.trim(0, 1, 0)
+    with pytest.raises(IndexError):
+        rect.rows
+    # The square rectangle passes its shape check but fails revalidation.
+    with pytest.raises(IndexError):
+        square.add_identity()
+    # A row handle taken before the shrink also raises on use.
+    rm = m.slice_mut(0, 1, 0, 3).row_mut(0)
+    # (still valid: row 0 survives)
+    rm.set_entry(0, 1)
+
+
+def test_augmented_segment_mutates_and_reads_back():
+    m = fp.AugmentedMatrix2(2, 2, [2, 2])
+    seg = m.segment(1, 1)
+    assert isinstance(seg, fp.MatrixSliceMut)
+    assert seg.rows == 2
+    assert seg.columns == 2
+    seg.add_identity()
+    start1 = m.segment_starts[1]
+    rows = m.to_py()
+    assert rows[0][start1] == 1
+    assert rows[1][start1 + 1] == 1
+
+
+def test_augmented_row_segment_mut_writes_through():
+    m = fp.AugmentedMatrix2(3, 2, [2, 2])
+    row = m.row_segment_mut(0, 0, 0)
+    assert isinstance(row, fp.FpSliceMut)
+    row.set_entry(0, 2)
+    assert m.to_py()[0][0] == 2
+
+
+def test_augmented3_segment_and_row_segment_mut_write_through():
+    # Exercises the MatrixParent::Augmented3 arm: take a segment / row segment
+    # from an AugmentedMatrix3, mutate through it, and observe the change via
+    # the augmented matrix (to_py / row_segment).
+    m = fp.AugmentedMatrix3(3, 2, [2, 2, 2])
+    seg = m.segment(1, 1)
+    assert isinstance(seg, fp.MatrixSliceMut)
+    assert seg.rows == 2
+    assert seg.columns == 2
+    seg.add_identity()
+    start1 = m.segment_starts[1]
+    rows = m.to_py()
+    assert rows[0][start1] == 1
+    assert rows[1][start1 + 1] == 1
+
+    row = m.row_segment_mut(0, 0, 0)
+    assert isinstance(row, fp.FpSliceMut)
+    row.set_entry(0, 2)
+    assert m.to_py()[0][0] == 2
+    # Observable via the read-only row_segment accessor too.
+    assert list(m.row_segment(0, 0, 0))[0] == 2
+
+
+def test_stale_slice_mut_repr_raises():
+    # __repr__ revalidates via with_slice_mut, so a stale handle raises rather
+    # than returning a string.
+    m = make_matrix([[1, 1, 0], [0, 1, 1]])
+    rect = m.slice_mut(0, 2, 0, 3)
+    m.trim(0, 1, 0)
+    with pytest.raises(IndexError):
+        repr(rect)
+
+
+def test_augmented_segment_builds_nontrivial_compute_values():
+    # Mirrors the F3 doctest of Matrix::compute_image / compute_quasi_inverse,
+    # but built entirely through the Python segment-mut API: place A in segment
+    # 0 entry by entry, the identity in segment 1, then row reduce and verify
+    # the committed expected image and preimage. This closes the earlier gap
+    # where the Python layer could not set interior augmented-matrix entries.
+    a = [
+        [1, 2, 1, 1, 0],
+        [1, 0, 2, 1, 1],
+        [2, 2, 0, 2, 1],
+    ]
+    m = fp.AugmentedMatrix2(3, 3, [5, 3])
+    # Fill A (segment 0) via row_segment_mut.
+    for i, arow in enumerate(a):
+        rm = m.row_segment_mut(i, 0, 0)
+        for j, v in enumerate(arow):
+            rm.set_entry(j, v)
+    # Identity into segment 1 via segment(...).
+    m.segment(1, 1).add_identity()
+    m.row_reduce()
+
+    image = m.compute_image()
+    assert image.dimension == 2
+    image_rows = [list(v) for v in image.iter()]
+    assert image_rows == [[1, 0, 2, 1, 1], [0, 1, 1, 0, 1]]
+
+    qi = m.compute_quasi_inverse()
+    assert qi.source_dimension == 3
+    assert qi.preimage.to_py() == [[0, 1, 0], [0, 2, 2]]
diff --git a/ext_py/tests/test_quasi_inverse.py b/ext_py/tests/test_quasi_inverse.py
new file mode 100644
index 0000000000..622b429f7a
--- /dev/null
+++ b/ext_py/tests/test_quasi_inverse.py
@@ -0,0 +1,251 @@
+import pytest
+
+from ext import fp
+
+
+def make_qi():
+    # Mirrors the Rust `test_stream_qi` example at p = 2.
+    preimage = fp.Matrix.from_py(
+        2,
+        [
+            [1, 0, 1, 1],
+            [1, 1, 0, 0],
+            [0, 1, 0, 1],
+            [1, 1, 1, 0],
+        ],
+    )
+    return fp.QuasiInverse([0, -1, 1, -1, 2, 3], preimage)
+
+
+def test_import_and_construction():
+    qi = make_qi()
+    assert qi.prime == 2
+    assert qi.image_dimension == 4
+    assert qi.source_dimension == 4
+    assert qi.target_dimension == 6
+    assert "QuasiInverse(2" in repr(qi)
+
+
+def test_pivots_and_preimage():
+    qi = make_qi()
+    assert qi.pivots == [0, -1, 1, -1, 2, 3]
+    assert qi.preimage.to_py() == [
+        [1, 0, 1, 1],
+        [1, 1, 0, 0],
+        [0, 1, 0, 1],
+        [1, 1, 1, 0],
+    ]
+
+
+def test_apply_known_example():
+    qi = make_qi()
+    v = fp.FpVector.from_py(2, [1, 1, 0, 0, 1, 0])
+    out = fp.FpVector(2, 4)
+    qi.apply(out, 1, v)
+    assert list(out) == [1, 1, 1, 0]
+
+
+def test_apply_accepts_slices():
+    qi = make_qi()
+    v = fp.FpVector.from_py(2, [1, 1, 0, 0, 1, 0])
+    out = fp.FpVector(2, 4)
+    qi.apply(out.slice_mut(0, 4), 1, v.slice(0, 6))
+    assert list(out) == [1, 1, 1, 0]
+
+
+def test_apply_dimension_and_prime_mismatch():
+    qi = make_qi()
+    out = fp.FpVector(2, 4)
+    # Wrong input length (target_dimension is 6).
+    bad_input = fp.FpVector(2, 5)
+    with pytest.raises(ValueError):
+        qi.apply(out, 1, bad_input)
+    # Wrong output length (source_dimension is 4).
+    good_input = fp.FpVector(2, 6)
+    with pytest.raises(ValueError):
+        qi.apply(fp.FpVector(2, 3), 1, good_input)
+    # Prime mismatch on input.
+    with pytest.raises(ValueError):
+        qi.apply(out, 1, fp.FpVector(3, 6))
+
+
+def test_apply_wrong_type():
+    qi = make_qi()
+    with pytest.raises(ValueError):
+        qi.apply(123, 1, fp.FpVector(2, 6))
+
+
+def make_identity_qi():
+    # Square identity preimage at p=3, so source==target==3 and a single
+    # length-3 vector is a valid argument both as input and as output target.
+    preimage = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    return fp.QuasiInverse([0, 1, 2], preimage)
+
+
+def test_apply_input_target_aliasing_raises_runtimeerror():
+    # Passing the SAME bare FpVector as both the mutable target and the input
+    # is an aliasing conflict: it must raise RuntimeError (borrow conflict),
+    # NOT the generic ValueError reserved for wrong-type arguments.
+    qi = make_identity_qi()
+    v = fp.FpVector.from_py(3, [1, 2, 1])
+    with pytest.raises(RuntimeError):
+        qi.apply(v, 1, v)
+    with pytest.raises(Exception) as excinfo:
+        qi.apply(v, 1, v)
+    assert not isinstance(excinfo.value, ValueError)
+
+
+def test_apply_wrong_type_is_valueerror_not_runtimeerror():
+    # A genuine wrong-type argument still raises ValueError (not RuntimeError).
+    qi = make_identity_qi()
+    good = fp.FpVector(3, 3)
+    with pytest.raises(ValueError):
+        qi.apply(123, 1, good)
+    with pytest.raises(ValueError):
+        qi.apply(fp.Matrix.from_py(3, [[1, 0, 0]]), 1, good)
+
+
+def test_apply_distinct_objects_regression():
+    # The normal distinct-objects call still produces the known value.
+    qi = make_identity_qi()
+    v = fp.FpVector.from_py(3, [1, 2, 1])
+    out = fp.FpVector(3, 3)
+    qi.apply(out, 2, v)
+    assert list(out) == [2, 1, 2]
+
+
+def test_bytes_roundtrip():
+    qi = make_qi()
+    data = qi.to_bytes()
+    assert isinstance(data, bytes)
+    restored = fp.QuasiInverse.from_bytes(2, data)
+    assert restored.source_dimension == qi.source_dimension
+    assert restored.target_dimension == qi.target_dimension
+    assert restored.image_dimension == qi.image_dimension
+
+    v = fp.FpVector.from_py(2, [1, 1, 0, 0, 1, 0])
+    out = fp.FpVector(2, 4)
+    restored.apply(out, 1, v)
+    assert list(out) == [1, 1, 1, 0]
+
+
+def test_from_bytes_malformed():
+    with pytest.raises(RuntimeError):
+        fp.QuasiInverse.from_bytes(2, b"\x00\x01\x02")
+
+
+def test_compute_quasi_inverse_from_matrix():
+    # Example from the Rust `compute_quasi_inverse` doc at p = 3.
+    rows = [
+        [1, 2, 1, 1, 0],
+        [1, 0, 2, 1, 1],
+        [2, 2, 0, 2, 1],
+    ]
+    m = fp.AugmentedMatrix2.from_py(3, rows)
+    m.row_reduce()
+    qi = m.compute_quasi_inverse()
+    assert qi.prime == 3
+    assert qi.source_dimension == 3
+    assert qi.preimage.to_py() == [[0, 1, 0], [0, 2, 2]]
+
+
+@pytest.mark.xfail(reason="panics")
+def test_compute_quasi_inverse_requires_row_reduce():
+    # Without row_reduce the matrix has uninitialized (empty) pivots, which
+    # upstream would slice out of bounds and panic. The guard turns that into
+    # a clean ValueError.
+    m = fp.AugmentedMatrix2.from_py(3, [[1, 0, 1], [0, 1, 1]])
+    with pytest.raises(ValueError):
+        m.compute_quasi_inverse()
+
+
+@pytest.mark.xfail(reason="panics")
+def test_inconsistent_image_raises():
+    # preimage has 4 rows; supply an image with 5 non-negative pivots.
+    preimage = fp.Matrix.from_py(
+        2,
+        [
+            [1, 0, 1, 1],
+            [1, 1, 0, 0],
+            [0, 1, 0, 1],
+            [1, 1, 1, 0],
+        ],
+    )
+    with pytest.raises(ValueError):
+        fp.QuasiInverse([0, 1, 2, 3, 0], preimage)
+
+
+@pytest.mark.xfail(reason="panics")
+def test_pivot_out_of_range_raises():
+    # preimage has 4 rows; a non-negative pivot of 4 is an invalid row index.
+    preimage = fp.Matrix.from_py(
+        2,
+        [
+            [1, 0, 1, 1],
+            [1, 1, 0, 0],
+            [0, 1, 0, 1],
+            [1, 1, 1, 0],
+        ],
+    )
+    with pytest.raises(ValueError):
+        fp.QuasiInverse([0, 1, 2, 4], preimage)
+
+
+def test_apply_with_coeff_odd_prime():
+    # At p = 3, identity-style preimage so the result is `coeff * input`.
+    preimage = fp.Matrix.from_py(
+        3,
+        [
+            [1, 0, 0],
+            [0, 1, 0],
+            [0, 0, 1],
+        ],
+    )
+    qi = fp.QuasiInverse([0, 1, 2], preimage)
+    v = fp.FpVector.from_py(3, [1, 2, 1])
+
+    out = fp.FpVector(3, 3)
+    qi.apply(out, 2, v)
+    # 2 * [1, 2, 1] mod 3 = [2, 1, 2]
+    assert list(out) == [2, 1, 2]
+
+    # A large coeff must reduce mod p and not overflow/panic.
+    out2 = fp.FpVector(3, 3)
+    qi.apply(out2, 0xFFFF_FFFF, v)
+    # 0xFFFFFFFF mod 3 == 0, so result is the zero vector.
+    assert list(out2) == [0, 0, 0]
+
+    out3 = fp.FpVector(3, 3)
+    qi.apply(out3, 2 + 3 * 10, v)  # coeff = 32 ≡ 2 mod 3
+    assert list(out3) == [2, 1, 2]
+
+
+def test_none_image_construction_and_roundtrip():
+    preimage = fp.Matrix.from_py(
+        2,
+        [
+            [1, 0, 0],
+            [0, 1, 0],
+            [0, 0, 1],
+        ],
+    )
+    qi = fp.QuasiInverse(None, preimage)
+    assert qi.pivots is None
+    # With a None (identity) image, target_dimension == image_dimension.
+    assert qi.image_dimension == 3
+    assert qi.source_dimension == 3
+    assert qi.target_dimension == 3
+
+    v = fp.FpVector.from_py(2, [1, 0, 1])
+    out = fp.FpVector(2, 3)
+    qi.apply(out, 1, v)
+    assert list(out) == [1, 0, 1]
+
+    # A None image is serialized as an explicit identity pivot list and so
+    # round-trips to Some([0, 1, 2, ...]) rather than None.
+    restored = fp.QuasiInverse.from_bytes(2, qi.to_bytes())
+    assert restored.pivots == [0, 1, 2]
+    assert restored.target_dimension == 3
+    out2 = fp.FpVector(2, 3)
+    restored.apply(out2, 1, v)
+    assert list(out2) == [1, 0, 1]
diff --git a/ext_py/tests/test_subquotient.py b/ext_py/tests/test_subquotient.py
new file mode 100644
index 0000000000..7bceb13a78
--- /dev/null
+++ b/ext_py/tests/test_subquotient.py
@@ -0,0 +1,228 @@
+import pytest
+
+from ext import fp
+
+
+def test_construction_and_queries():
+    sq = fp.Subquotient(3, 5)
+    assert sq.prime == 3
+    assert isinstance(sq.prime, int)
+    assert sq.ambient_dimension == 5
+    assert sq.dimension == 0
+    assert len(sq) == 0
+    assert sq.is_empty
+    assert repr(sq) == "Subquotient(3, dim=0, ambient=5)"
+
+
+def test_new_full():
+    sq = fp.Subquotient.new_full(2, 4)
+    assert sq.dimension == 4
+    assert sq.ambient_dimension == 4
+    assert sq.quotient_dimension == 4
+    assert len(sq.gens) == 4
+
+
+def test_invalid_prime_raises():
+    with pytest.raises(ValueError):
+        fp.Subquotient(4, 3)
+    with pytest.raises(ValueError):
+        fp.Subquotient.new_full(4, 3)
+
+
+def test_add_gen_quotient_reduce_and_gens():
+    # Mirrors the upstream `test_add_gen` example at p = 3, dim = 5.
+    sq = fp.Subquotient(3, 5)
+    sq.quotient(fp.FpVector.from_py(3, [1, 1, 0, 0, 1]))
+    sq.quotient(fp.FpVector.from_py(3, [0, 2, 0, 0, 1]))
+    sq.add_gen(fp.FpVector.from_py(3, [1, 1, 0, 0, 0]))
+    sq.add_gen(fp.FpVector.from_py(3, [0, 1, 0, 0, 0]))
+
+    assert sq.dimension == 1
+    gens = sq.gens
+    assert len(gens) == 1
+    assert list(gens[0]) == [0, 0, 0, 0, 1]
+
+    zeros = sq.zeros()
+    assert isinstance(zeros, fp.Subspace)
+    assert zeros.dimension == 2
+
+    # reduce returns the coefficients w.r.t. the generators and mutates the
+    # vector in place.
+    elt = fp.FpVector.from_py(3, [2, 0, 0, 0, 0])
+    assert sq.reduce(elt) == [2]
+
+    # complement + quotient + gens cover the ambient space.
+    assert (
+        sq.zeros().dimension + len(sq.gens) + len(sq.complement_pivots)
+        == sq.ambient_dimension
+    )
+
+
+def test_subspace_gens_quotient_pivots_and_dimension():
+    # Non-trivial subquotient (non-empty quotient and one generator), mirroring
+    # the upstream `test_add_gen` example at p = 3, dim = 5. After the calls the
+    # upstream `Display` is:
+    #   Generators: [0, 0, 0, 0, 1]
+    #   Zeros:      [1, 0, 0, 0, 2]
+    #               [0, 1, 0, 0, 2]
+    sq = fp.Subquotient(3, 5)
+    sq.quotient(fp.FpVector.from_py(3, [1, 1, 0, 0, 1]))
+    sq.quotient(fp.FpVector.from_py(3, [0, 2, 0, 0, 1]))
+    sq.add_gen(fp.FpVector.from_py(3, [1, 1, 0, 0, 0]))
+    sq.add_gen(fp.FpVector.from_py(3, [0, 1, 0, 0, 0]))
+
+    # dimension is the subspace-part generator count; quotient (zeros) dim is 2.
+    assert sq.dimension == 1
+    assert sq.zeros().dimension == 2
+
+    # subspace_dimension == self.dimension + quotient.dimension per upstream
+    # `subquotient.rs::subspace_dimension`.
+    assert sq.subspace_dimension == sq.dimension + sq.zeros().dimension
+    assert sq.subspace_dimension == 3
+
+    # subspace_gens chains gens() with the quotient's basis vectors (upstream
+    # `subspace_gens` = `gens().chain(quotient.iter())`).
+    subspace_gens = [list(v) for v in sq.subspace_gens]
+    assert subspace_gens == [
+        [0, 0, 0, 0, 1],
+        [1, 0, 0, 0, 2],
+        [0, 1, 0, 0, 2],
+    ]
+
+    # quotient_pivots is the quotient subspace's pivot table: pivots[col] = row
+    # index of the pivot in that column, else -1. Quotient pivots are in cols 0,1.
+    assert sq.quotient_pivots == [0, 1, -1, -1, -1]
+
+
+def test_clear_gens_keeps_quotient():
+    sq = fp.Subquotient(3, 5)
+    sq.quotient(fp.FpVector.from_py(3, [1, 1, 0, 0, 1]))
+    sq.add_gen(fp.FpVector.from_py(3, [0, 1, 0, 0, 0]))
+    assert sq.dimension >= 1
+    sq.clear_gens()
+    assert sq.dimension == 0
+    assert sq.zeros().dimension == 1
+
+
+def test_set_to_full():
+    sq = fp.Subquotient(2, 3)
+    sq.set_to_full()
+    # `set_to_full` makes the gens the entire space and clears the quotient,
+    # but (matching upstream) does not update the cached `dimension` counter.
+    assert sq.zeros().dimension == 0
+    assert len(sq.gens) == 3
+
+    # Stale-`dimension` quirk: `set_to_full` makes gens the entire space and
+    # clears the quotient, but upstream does NOT update the cached `dimension`
+    # counter. So on a fresh Subquotient these are inconsistent today:
+    #   dimension()/len(sq) report 0 (stale), while gens() actually has 3 rows.
+    # Pin the surprising current behavior so a future upstream fix (syncing the
+    # cached dimension) trips this test and prompts a revisit.
+    assert sq.dimension == 0
+    assert len(sq) == 0
+    assert len(sq.gens) == 3
+
+
+def test_from_parts():
+    sub = fp.Subspace(2, 3)
+    sub.add_vector(fp.FpVector.from_py(2, [1, 0, 0]))
+    sub.add_vector(fp.FpVector.from_py(2, [0, 1, 0]))
+    quot = fp.Subspace(2, 3)
+    quot.add_vector(fp.FpVector.from_py(2, [1, 0, 0]))
+
+    sq = fp.Subquotient.from_parts(sub, quot)
+    assert sq.dimension == 1
+    assert sq.ambient_dimension == 3
+
+
+def test_from_parts_mismatch_raises():
+    sub = fp.Subspace(2, 3)
+    bad = fp.Subspace(2, 4)
+    with pytest.raises(ValueError):
+        fp.Subquotient.from_parts(sub, bad)
+    other_prime = fp.Subspace(3, 3)
+    with pytest.raises(ValueError):
+        fp.Subquotient.from_parts(sub, other_prime)
+
+
+def test_invalid_vector_inputs_raise():
+    sq = fp.Subquotient(3, 3)
+    with pytest.raises(ValueError):
+        sq.quotient(fp.FpVector.from_py(5, [1, 0, 0]))
+    with pytest.raises(ValueError):
+        sq.add_gen(fp.FpVector.from_py(3, [1, 0]))
+    with pytest.raises(ValueError):
+        sq.reduce(fp.FpVector.from_py(3, [1, 0]))
+
+
+def test_reduce_by_quotient():
+    sq = fp.Subquotient(3, 3)
+    sq.quotient(fp.FpVector.from_py(3, [1, 0, 0]))
+    v = fp.FpVector.from_py(3, [1, 1, 0])
+    sq.reduce_by_quotient(v)
+    assert list(v) == [0, 1, 0]
+
+
+def test_reduce_by_quotient_slice_mut():
+    sq = fp.Subquotient(3, 3)
+    sq.quotient(fp.FpVector.from_py(3, [1, 0, 0]))
+    # Full-row slice: reduction must write through to the matrix.
+    m = fp.Matrix.from_py(3, [[1, 1, 0]])
+    sq.reduce_by_quotient(m.row_mut(0))
+    assert list(m.row(0)) == [0, 1, 0]
+    # Sub-range slice: only the targeted columns are reduced in place.
+    m2 = fp.Matrix.from_py(3, [[2, 1, 1, 0]])
+    row = m2.row_mut(0)
+    sq.reduce_by_quotient(row.slice_mut(1, 4))
+    assert list(m2.row(0)) == [2, 0, 1, 0]
+
+
+def test_reduce_matrix():
+    source = fp.Subquotient.new_full(3, 2)
+    target = fp.Subquotient.new_full(3, 2)
+    # identity matrix maps source ambient (rows) to target ambient (cols).
+    m = fp.Matrix.from_py(3, [[1, 0], [0, 1]])
+    result = fp.Subquotient.reduce_matrix(m, source, target)
+    assert len(result) == source.dimension
+
+
+def test_reduce_matrix_values_with_nontrivial_quotient():
+    # source = full space of dim 2 at p = 3, so gens() = [1, 0] and [0, 1].
+    source = fp.Subquotient.new_full(3, 2)
+
+    # target has a non-trivial quotient: quotient kills coordinate 1, generator
+    # is [1, 0]. So target.reduce projects out column 1 and reads coeff at col 0.
+    target = fp.Subquotient(3, 2)
+    target.quotient(fp.FpVector.from_py(3, [0, 1]))
+    target.add_gen(fp.FpVector.from_py(3, [1, 0]))
+
+    # Non-identity matrix. `Matrix.apply` computes (input row-vector) * matrix:
+    # result = sum_i input[i] * row(i). So gen [1,0] -> row 0 = [2, 1];
+    # gen [0,1] -> row 1 = [0, 1].
+    m = fp.Matrix.from_py(3, [[2, 1], [0, 1]])
+
+    # Reducing images in target: drop col 1 then read coeff at col 0.
+    #   [2, 1] -> quotient -> [2, 0] -> coeff [2]
+    #   [0, 1] -> quotient -> [0, 0] -> coeff [0]
+    result = fp.Subquotient.reduce_matrix(m, source, target)
+    assert result == [[2], [0]]
+
+
+def test_reduce_matrix_dimension_mismatches_raise():
+    source = fp.Subquotient.new_full(3, 2)
+    target = fp.Subquotient.new_full(3, 2)
+
+    # rows != source.ambient_dimension (3 rows vs ambient 2).
+    bad_rows = fp.Matrix.from_py(3, [[1, 0], [0, 1], [0, 0]])
+    with pytest.raises(ValueError):
+        fp.Subquotient.reduce_matrix(bad_rows, source, target)
+
+    # columns != target.ambient_dimension (3 cols vs ambient 2).
+    bad_cols = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 0]])
+    with pytest.raises(ValueError):
+        fp.Subquotient.reduce_matrix(bad_cols, source, target)
+
+    # prime mismatch between matrix and subquotients.
+    bad_prime = fp.Matrix.from_py(5, [[1, 0], [0, 1]])
+    with pytest.raises(ValueError):
+        fp.Subquotient.reduce_matrix(bad_prime, source, target)
diff --git a/ext_py/tests/test_subspace.py b/ext_py/tests/test_subspace.py
new file mode 100644
index 0000000000..b564efc070
--- /dev/null
+++ b/ext_py/tests/test_subspace.py
@@ -0,0 +1,173 @@
+import pytest
+
+from ext import fp
+
+
+def test_subspace_construction_and_queries():
+    s = fp.Subspace(3, 3)
+    assert s.prime == 3
+    assert s.ambient_dimension == 3
+    assert s.dimension == 0
+    assert len(s) == 0
+    assert repr(s) == "Subspace(3, dim=0, ambient=3)"
+
+
+def test_subspace_entire_space():
+    s = fp.Subspace.entire_space(2, 3)
+    assert s.dimension == 3
+    assert s.ambient_dimension == 3
+
+
+def test_subspace_from_matrix():
+    m = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 2]])
+    s = fp.Subspace.from_matrix(m)
+    assert s.dimension == 2
+    assert s.ambient_dimension == 3
+
+
+def test_add_vector_and_contains():
+    s = fp.Subspace(3, 3)
+    v = fp.FpVector.from_py(3, [1, 0, 0])
+    assert s.add_vector(v) == 1
+    assert s.contains(v)
+    assert v in s
+
+    w = fp.FpVector.from_py(3, [0, 1, 0])
+    assert not s.contains(w)
+    assert w not in s
+
+    # A scalar multiple is still contained.
+    twice = fp.FpVector.from_py(3, [2, 0, 0])
+    assert s.contains(twice)
+
+
+def test_contains_space_and_sum():
+    a = fp.Subspace(3, 3)
+    a.add_vector(fp.FpVector.from_py(3, [1, 0, 0]))
+    b = fp.Subspace(3, 3)
+    b.add_vector(fp.FpVector.from_py(3, [0, 1, 0]))
+
+    assert not a.contains_space(b)
+
+    # Incompatible prime/ambient dimension raise ValueError rather than panic.
+    with pytest.raises(ValueError):
+        a.contains_space(fp.Subspace(5, 3))
+    with pytest.raises(ValueError):
+        a.contains_space(fp.Subspace(3, 4))
+
+    # The sum of two complementary lines is their 2-dimensional span.
+    s = a.sum(b)
+    assert s.dimension == 2
+    assert s.contains_space(a)
+    assert s.contains_space(b)
+    assert s.contains(fp.FpVector.from_py(3, [1, 0, 0]))
+    assert s.contains(fp.FpVector.from_py(3, [0, 1, 0]))
+    assert s.dimension <= s.ambient_dimension
+
+    # Overlapping subspaces: the sum's dimension is the union's rank.
+    c = fp.Subspace(3, 3)
+    c.add_vector(fp.FpVector.from_py(3, [1, 0, 0]))
+    overlap = a.sum(c)
+    assert overlap.dimension == 1
+    assert overlap.contains_space(a)
+
+
+def test_reduce_in_place():
+    s = fp.Subspace(3, 3)
+    s.add_vector(fp.FpVector.from_py(3, [1, 0, 0]))
+
+    v = fp.FpVector.from_py(3, [2, 1, 0])
+    s.reduce(v)
+    assert v[0] == 0
+    assert v[1] == 1
+
+
+def test_set_to_zero_and_entire():
+    s = fp.Subspace.entire_space(2, 3)
+    assert s.dimension == 3
+    s.set_to_zero()
+    assert s.dimension == 0
+    s.set_to_entire()
+    assert s.dimension == 3
+
+
+def test_iter_returns_basis_vectors():
+    m = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 2]])
+    s = fp.Subspace.from_matrix(m)
+    basis = s.iter()
+    assert [list(v) for v in basis] == [[1, 0, 0], [0, 1, 2]]
+
+
+def test_basis_matches_iter():
+    m = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 2]])
+    s = fp.Subspace.from_matrix(m)
+    assert [list(v) for v in s.basis] == [[1, 0, 0], [0, 1, 2]]
+    assert [list(v) for v in s.basis] == [list(v) for v in s.iter()]
+
+
+def test_contains_accepts_slice():
+    m = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 2]])
+    s = fp.Subspace.from_matrix(m)
+    v = fp.FpVector.from_py(3, [1, 0, 0])
+    # An owned FpVector and an FpSlice over it behave identically.
+    assert s.contains(v)
+    assert s.contains(v.slice(0, 3))
+    w = fp.FpVector.from_py(3, [0, 0, 1])
+    assert not s.contains(w.slice(0, 3))
+    # Slice of the wrong ambient dimension is rejected.
+    with pytest.raises(ValueError):
+        s.contains(w.slice(0, 2))
+
+
+def test_iter_all_vectors():
+    m = fp.Matrix.from_py(3, [[1, 0, 0], [0, 1, 2]])
+    s = fp.Subspace.from_matrix(m)
+
+    # iter_all_vectors returns a lazy iterator, not a list.
+    it = s.iter_all_vectors()
+    assert not isinstance(it, list)
+    assert iter(it) is it
+
+    # Iterating via a for-loop yields every vector in the subspace.
+    collected = []
+    for v in s.iter_all_vectors():
+        collected.append(list(v))
+    assert len(collected) == 9
+
+    # list(...) over the iterator yields the same set/count.
+    vectors = sorted(list(v) for v in s.iter_all_vectors())
+    assert len(vectors) == 9
+    assert sorted(collected) == vectors
+    assert [0, 0, 0] in vectors
+    assert [1, 1, 2] in vectors
+
+
+def test_bytes_roundtrip():
+    s = fp.Subspace(3, 3)
+    s.add_vector(fp.FpVector.from_py(3, [1, 0, 0]))
+    data = s.to_bytes()
+    assert isinstance(data, bytes)
+    restored = fp.Subspace.from_bytes(3, data)
+    assert restored.dimension == 1
+    assert restored.contains(fp.FpVector.from_py(3, [1, 0, 0]))
+
+
+def test_invalid_inputs_raise():
+    s = fp.Subspace(3, 3)
+
+    with pytest.raises(ValueError):
+        s.contains(fp.FpVector.from_py(5, [1, 0, 0]))
+
+    with pytest.raises(ValueError):
+        s.add_vector(fp.FpVector.from_py(3, [1, 0]))
+
+    with pytest.raises(ValueError):
+        fp.Subspace(1, 3)
+
+
+def test_bytes_error_taxonomy_is_consistent():
+    # Subspace mirrors FpVector/Matrix: both to_bytes and from_bytes map
+    # serialization failures to RuntimeError, so malformed bytes raise
+    # RuntimeError (not ValueError).
+    with pytest.raises(RuntimeError):
+        fp.Subspace.from_bytes(3, b"\x00\x01\x02")

From 32eb39925fa00b412224a6ecba983ed7ef795cc1 Mon Sep 17 00:00:00 2001
From: Hood Chatham <roberthoodchatham@gmail.com>
Date: Thu, 9 Jul 2026 18:11:00 -0700
Subject: [PATCH 2/4] Add mut and const, remove as_slice and as_slice_mut

---
 ext_py/src/fp_mod/vectors.rs          | 28 ++++++++++++++++-----------
 ext_py/tests/test_fp_vector_slices.py |  7 +++----
 2 files changed, 20 insertions(+), 15 deletions(-)

diff --git a/ext_py/src/fp_mod/vectors.rs b/ext_py/src/fp_mod/vectors.rs
index 7df73fa782..fcf7b91c04 100644
--- a/ext_py/src/fp_mod/vectors.rs
+++ b/ext_py/src/fp_mod/vectors.rs
@@ -219,9 +219,23 @@ impl PyFpVector {
         })
     }
 
-    /// A mutable `FpSliceMut` spanning the whole vector (mirrors upstream
-    /// `FpVector::as_slice_mut`); equivalent to `slice_mut(0, len())`.
-    pub fn as_slice_mut(slf: PyRef<'_, Self>) -> PyResult<PyFpSliceMut> {
+    /// A read-only `FpSlice` spanning the whole vector; equivalent to
+    /// `slice(0, len())`.
+    #[getter]
+    pub fn r#const(slf: PyRef<'_, Self>) -> PyResult<PyFpSlice> {
+        let end = slf.0.len();
+        let py = slf.py();
+        Ok(PyFpSlice {
+            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
+            start: 0,
+            end,
+        })
+    }
+
+    /// A mutable `FpSliceMut` spanning the whole vector; equivalent to
+    /// `slice_mut(0, len())`.
+    #[getter]
+    pub fn r#mut(slf: PyRef<'_, Self>) -> PyResult<PyFpSliceMut> {
         let end = slf.0.len();
         let py = slf.py();
         Ok(PyFpSliceMut {
@@ -566,14 +580,6 @@ impl PyFpSliceMut {
         self.with_slice_mut(py, |mut s| s.add_basis_element(index, value))
     }
 
-    pub fn as_slice(&self, py: Python<'_>) -> PyFpSlice {
-        PyFpSlice {
-            parent: self.parent.clone_ref(py),
-            start: self.start,
-            end: self.end,
-        }
-    }
-
     pub fn slice_mut(&self, py: Python<'_>, start: usize, end: usize) -> PyResult<Self> {
         checked_range(start, end, self.span())?;
         Ok(Self {
diff --git a/ext_py/tests/test_fp_vector_slices.py b/ext_py/tests/test_fp_vector_slices.py
index 28ffd0ec1a..713194d93a 100644
--- a/ext_py/tests/test_fp_vector_slices.py
+++ b/ext_py/tests/test_fp_vector_slices.py
@@ -93,7 +93,7 @@ def test_fp_slice_handles_parent_lifetime_after_original_reference_deleted():
     assert list(s.iter()) == [2, 3]
     sm.add_basis_element(0, 4)
     assert s[1] == 2
-    assert list(sm.as_slice().iter()) == [2, 4]
+    assert list(sm.iter()) == [2, 4]
 
 
 def test_stale_slice_handles_raise_python_exception_after_parent_shrink():
@@ -139,9 +139,8 @@ def test_fp_slice_mut_updates_parent_and_as_slice():
     s.scale(2)
     assert [v[i] for i in range(len(v))] == [1, 4, 3, 4, 0]
 
-    as_slice = s.as_slice()
-    assert len(as_slice) == 3
-    assert [as_slice[i] for i in range(len(as_slice))] == [4, 3, 4]
+    assert len(s) == 3
+    assert [s[i] for i in range(len(s))] == [4, 3, 4]
 
     sub = s.slice_mut(1, 3)
     sub.set_to_zero()

From 98a90b64e34c71309430b1a38b0c6fea391b2720 Mon Sep 17 00:00:00 2001
From: Hood Chatham <roberthoodchatham@gmail.com>
Date: Thu, 9 Jul 2026 18:14:48 -0700
Subject: [PATCH 3/4] Make FpVectorSlice and FpVectorSliceMut iterable

---
 ext_py/src/fp_mod/vectors.rs          | 15 ++++++++++-----
 ext_py/tests/test_fp_vector_slices.py |  6 +++---
 ext_py/tests/test_matrix.py           |  4 ++--
 3 files changed, 15 insertions(+), 10 deletions(-)

diff --git a/ext_py/src/fp_mod/vectors.rs b/ext_py/src/fp_mod/vectors.rs
index fcf7b91c04..3654d23978 100644
--- a/ext_py/src/fp_mod/vectors.rs
+++ b/ext_py/src/fp_mod/vectors.rs
@@ -360,11 +360,6 @@ impl PyFpSlice {
         self.with_slice(py, |s| s.entry(index))
     }
 
-    pub fn iter(&self, py: Python<'_>) -> PyResult<PyFpVectorIterator> {
-        let entries = self.with_slice(py, |s| s.iter().collect())?;
-        Ok(PyFpVectorIterator { entries, index: 0 })
-    }
-
     pub fn iter_nonzero(&self, py: Python<'_>) -> PyResult<Vec<(usize, u32)>> {
         self.with_slice(py, |s| s.iter_nonzero().collect())
     }
@@ -401,6 +396,11 @@ impl PyFpSlice {
         self.with_slice(py, |s| s.entry(index))
     }
 
+    pub fn __iter__(&self, py: Python<'_>) -> PyResult<PyFpVectorIterator> {
+        let entries = self.with_slice(py, |s| s.iter().collect())?;
+        Ok(PyFpVectorIterator { entries, index: 0 })
+    }
+
     pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
         self.with_slice(py, |s| format!("FpSlice({}, {})", s.prime().as_u32(), s))
     }
@@ -607,6 +607,11 @@ impl PyFpSliceMut {
         self.with_slice_mut(py, |mut s| s.set_entry(index, value))
     }
 
+    pub fn __iter__(&self, py: Python<'_>) -> PyResult<PyFpVectorIterator> {
+        let entries = self.with_slice(py, |s| s.iter().collect())?;
+        Ok(PyFpVectorIterator { entries, index: 0 })
+    }
+
     pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
         self.with_slice(py, |s| format!("FpSliceMut({}, {})", s.prime().as_u32(), s))
     }
diff --git a/ext_py/tests/test_fp_vector_slices.py b/ext_py/tests/test_fp_vector_slices.py
index 713194d93a..31a6782880 100644
--- a/ext_py/tests/test_fp_vector_slices.py
+++ b/ext_py/tests/test_fp_vector_slices.py
@@ -16,7 +16,7 @@ def test_fp_vector_slice_queries_and_to_owned():
     assert s[-1] == 0
     assert not s.is_zero
     assert s.first_nonzero == (0, 1)
-    assert list(s.iter()) == [1, 2, 0]
+    assert list(s) == [1, 2, 0]
     assert s.iter_nonzero() == [(0, 1), (1, 2)]
     assert repr(s) == "FpSlice(5, [1, 2, 0])"
 
@@ -90,10 +90,10 @@ def test_fp_slice_handles_parent_lifetime_after_original_reference_deleted():
 
     del v
 
-    assert list(s.iter()) == [2, 3]
+    assert list(s) == [2, 3]
     sm.add_basis_element(0, 4)
     assert s[1] == 2
-    assert list(sm.iter()) == [2, 4]
+    assert list(sm) == [2, 4]
 
 
 def test_stale_slice_handles_raise_python_exception_after_parent_shrink():
diff --git a/ext_py/tests/test_matrix.py b/ext_py/tests/test_matrix.py
index 46a83a9eed..77cdf25995 100644
--- a/ext_py/tests/test_matrix.py
+++ b/ext_py/tests/test_matrix.py
@@ -82,10 +82,10 @@ def test_row_access_and_getitem():
     assert row[-1] == 1
     assert not row.is_zero
     assert row.first_nonzero == (0, 4)
-    assert list(row.iter()) == [4, 0, 1]
+    assert list(row) == [4, 0, 1]
     assert row.iter_nonzero() == [(0, 4), (2, 1)]
     assert m[0].to_owned().prime == 5
-    assert list(m[0].iter()) == [1, 2, 3]
+    assert list(m[0]) == [1, 2, 3]
 
     with pytest.raises(IndexError):
         m.row(2)

From 9facd7abcf728dd8e0c076c651343bd911e79687 Mon Sep 17 00:00:00 2001
From: Hood Chatham <roberthoodchatham@gmail.com>
Date: Thu, 9 Jul 2026 21:55:48 -0700
Subject: [PATCH 4/4] Use .mut and .const and getitem slices instead of
 .slice() and .slice_mut()

---
 ext_py/src/fp_mod.rs                  |   6 +-
 ext_py/src/fp_mod/matrices.rs         |   7 +-
 ext_py/src/fp_mod/vectors.rs          | 109 ++++++++++++-------------
 ext_py/tests/test_affine_subspace.py  |   2 +-
 ext_py/tests/test_augmented_matrix.py |   4 +-
 ext_py/tests/test_fp_vector_slices.py | 113 +++++++++-----------------
 ext_py/tests/test_matrix.py           |  20 ++---
 ext_py/tests/test_quasi_inverse.py    |   2 +-
 ext_py/tests/test_subquotient.py      |   2 +-
 ext_py/tests/test_subspace.py         |   6 +-
 10 files changed, 117 insertions(+), 154 deletions(-)

diff --git a/ext_py/src/fp_mod.rs b/ext_py/src/fp_mod.rs
index fb62c387f5..c0db24c983 100644
--- a/ext_py/src/fp_mod.rs
+++ b/ext_py/src/fp_mod.rs
@@ -294,10 +294,12 @@ pub mod fp_py {
         PyAffineSubspace, PyMatrix, PyMatrixSliceMut, PyQuasiInverse, PySubquotient, PySubspace,
         PySubspaceVectorIterator,
     };
-    use super::matrices::{PyAugmentedMatrix2, PyAugmentedMatrix3};
     #[pymodule_export]
     pub use super::vectors::{PyFpSlice, PyFpSliceMut, PyFpVector, PyFpVectorIterator};
-    use super::*;
+    use super::{
+        matrices::{PyAugmentedMatrix2, PyAugmentedMatrix3},
+        *,
+    };
 
     #[pyclass(name = "Fp", frozen, from_py_object)]
     #[derive(Clone, Copy)]
diff --git a/ext_py/src/fp_mod/matrices.rs b/ext_py/src/fp_mod/matrices.rs
index 2b0c01d7a4..6511a4d5e0 100644
--- a/ext_py/src/fp_mod/matrices.rs
+++ b/ext_py/src/fp_mod/matrices.rs
@@ -234,7 +234,7 @@ impl PyMatrixSliceMut {
     }
 
     fn to_py(&self, py: Python<'_>) -> PyResult<Vec<Vec<u32>>> {
-        Ok(self.with_slice_mut(py, |s| s.to_vec())?)
+        self.with_slice_mut(py, |s| s.to_vec())
     }
 
     pub fn __repr__(&self, py: Python<'_>) -> PyResult<String> {
@@ -1278,7 +1278,10 @@ augmented_matrix_pyclass!(PyAugmentedMatrix2, "AugmentedMatrix2", 2, Augmented2,
         for row in &input {
             checked_equal_len(row.len(), columns)?;
         }
-        Ok(Self(Consumable::new("AugmentedMatrix2", RustMatrix::augmented_from_vec(p, &input))))
+        Ok(Self(Consumable::new(
+            "AugmentedMatrix2",
+            RustMatrix::augmented_from_vec(p, &input),
+        )))
     }
 });
 
diff --git a/ext_py/src/fp_mod/vectors.rs b/ext_py/src/fp_mod/vectors.rs
index 3654d23978..560a0f73f3 100644
--- a/ext_py/src/fp_mod/vectors.rs
+++ b/ext_py/src/fp_mod/vectors.rs
@@ -4,7 +4,11 @@ use fp::{
     prime::Prime,
     vector::{FpSlice as RustFpSlice, FpSliceMut as RustFpSliceMut, FpVector as RustFpVector},
 };
-use pyo3::{exceptions::PyValueError, types::PyBytes};
+use pyo3::{
+    exceptions::PyValueError,
+    types::{PyAny, PyAnyMethods, PyBytes, PySlice},
+    IntoPyObject, Py,
+};
 
 use super::*;
 
@@ -108,6 +112,23 @@ pub(crate) fn py_index(index: isize, len: usize) -> PyResult<usize> {
     }
 }
 
+/// Resolve a Python subscript that may be either an integer index or a
+/// step-1 `slice`. Returns `Some((start, end))` for a slice (after clamping
+/// to `len`), or `None` for a plain integer (left for the caller to
+/// normalize via [`py_index`]). Slices with a step other than 1 raise
+/// `ValueError`.
+fn slice_step_one(obj: &Bound<'_, PyAny>, len: usize) -> PyResult<Option<(usize, usize)>> {
+    if let Ok(slice) = obj.cast::<PySlice>() {
+        let indices = slice.indices(len as isize)?;
+        if indices.step != 1 {
+            return Err(PyValueError::new_err("slice step must be 1"));
+        }
+        Ok(Some((indices.start as usize, indices.stop as usize)))
+    } else {
+        Ok(None)
+    }
+}
+
 #[pyclass(name = "FpVector")]
 pub struct PyFpVector(pub(crate) RustFpVector);
 
@@ -185,40 +206,6 @@ impl PyFpVector {
         self.0.as_slice().iter_nonzero().collect()
     }
 
-    pub fn slice(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSlice> {
-        checked_range(start, end, slf.0.len())?;
-        let py = slf.py();
-        Ok(PyFpSlice {
-            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
-            start,
-            end,
-        })
-    }
-
-    /// Restrict to the sub-range of coordinates `[start, end)`, returning a
-    /// read-only `FpSlice` view. For an `FpVector` (whose coordinates start
-    /// at 0) this mirrors `slice(start, end)`; named to match the analogous
-    /// `FpSlice.restrict` method.
-    pub fn restrict(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSlice> {
-        checked_range(start, end, slf.0.len())?;
-        let py = slf.py();
-        Ok(PyFpSlice {
-            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
-            start,
-            end,
-        })
-    }
-
-    pub fn slice_mut(slf: PyRef<'_, Self>, start: usize, end: usize) -> PyResult<PyFpSliceMut> {
-        checked_range(start, end, slf.0.len())?;
-        let py = slf.py();
-        Ok(PyFpSliceMut {
-            parent: SliceParent::Vector(slf.into_pyobject(py)?.unbind()),
-            start,
-            end,
-        })
-    }
-
     /// A read-only `FpSlice` spanning the whole vector; equivalent to
     /// `slice(0, len())`.
     #[getter]
@@ -355,11 +342,6 @@ impl PyFpSlice {
         self.with_slice(py, |s| s.is_empty())
     }
 
-    pub fn entry(&self, py: Python<'_>, index: usize) -> PyResult<u32> {
-        let index = checked_index(index, self.span())?;
-        self.with_slice(py, |s| s.entry(index))
-    }
-
     pub fn iter_nonzero(&self, py: Python<'_>) -> PyResult<Vec<(usize, u32)>> {
         self.with_slice(py, |s| s.iter_nonzero().collect())
     }
@@ -391,9 +373,18 @@ impl PyFpSlice {
         self.len(py)
     }
 
-    pub fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult<u32> {
-        let index = py_index(index, self.span())?;
-        self.with_slice(py, |s| s.entry(index))
+    pub fn __getitem__(&self, py: Python<'_>, index: &Bound<'_, PyAny>) -> PyResult<Py<PyAny>> {
+        if let Some((start, end)) = slice_step_one(index, self.span())? {
+            let sub = self.restrict(py, start, end)?;
+            Ok(Py::new(py, sub)?.into_any())
+        } else {
+            let i = py_index(index.extract::<isize>()?, self.span())?;
+            Ok(self
+                .with_slice(py, |s| s.entry(i))?
+                .into_pyobject(py)?
+                .into_any()
+                .unbind())
+        }
     }
 
     pub fn __iter__(&self, py: Python<'_>) -> PyResult<PyFpVectorIterator> {
@@ -580,15 +571,6 @@ impl PyFpSliceMut {
         self.with_slice_mut(py, |mut s| s.add_basis_element(index, value))
     }
 
-    pub fn slice_mut(&self, py: Python<'_>, start: usize, end: usize) -> PyResult<Self> {
-        checked_range(start, end, self.span())?;
-        Ok(Self {
-            parent: self.parent.clone_ref(py),
-            start: self.start + start,
-            end: self.start + end,
-        })
-    }
-
     pub fn to_owned(&self, py: Python<'_>) -> PyResult<PyFpVector> {
         Ok(PyFpVector(self.with_slice(py, |s| s.to_owned())?))
     }
@@ -597,9 +579,26 @@ impl PyFpSliceMut {
         self.len(py)
     }
 
-    pub fn __getitem__(&self, py: Python<'_>, index: isize) -> PyResult<u32> {
-        let index = py_index(index, self.span())?;
-        self.with_slice(py, |s| s.entry(index))
+    pub fn __getitem__(&self, py: Python<'_>, index: &Bound<'_, PyAny>) -> PyResult<Py<PyAny>> {
+        if let Some((start, end)) = slice_step_one(index, self.span())? {
+            checked_range(start, end, self.span())?;
+            Ok(Py::new(
+                py,
+                Self {
+                    parent: self.parent.clone_ref(py),
+                    start: self.start + start,
+                    end: self.start + end,
+                },
+            )?
+            .into_any())
+        } else {
+            let i = py_index(index.extract::<isize>()?, self.span())?;
+            Ok(self
+                .with_slice(py, |s| s.entry(i))?
+                .into_pyobject(py)?
+                .into_any()
+                .unbind())
+        }
     }
 
     pub fn __setitem__(&self, py: Python<'_>, index: isize, value: u32) -> PyResult<()> {
diff --git a/ext_py/tests/test_affine_subspace.py b/ext_py/tests/test_affine_subspace.py
index d1edf60270..b26fbb7713 100644
--- a/ext_py/tests/test_affine_subspace.py
+++ b/ext_py/tests/test_affine_subspace.py
@@ -88,7 +88,7 @@ def test_contains_accepts_slice():
     aff = fp.AffineSubspace(offset, linear)
 
     v = fp.FpVector.from_py(2, [1, 1, 0])
-    assert aff.contains(v.slice(0, 3))
+    assert aff.contains(v.const[0:3])
 
 
 def test_contains_space():
diff --git a/ext_py/tests/test_augmented_matrix.py b/ext_py/tests/test_augmented_matrix.py
index 4f1d4349aa..e78e5cfa4b 100644
--- a/ext_py/tests/test_augmented_matrix.py
+++ b/ext_py/tests/test_augmented_matrix.py
@@ -77,8 +77,8 @@ def test_row_mut_writes_through_to_parent():
     rm = m.row_mut(0)
     assert isinstance(rm, fp.FpSliceMut)
     assert len(rm) == m.columns
-    rm.set_entry(0, 1)
-    rm.slice_mut(2, 4).set_entry(1, 1)
+    rm[0] = 1
+    rm[2:4][1] = 1
     assert m.to_py()[0][0] == 1
     assert m.to_py()[0][3] == 1
     with pytest.raises(IndexError):
diff --git a/ext_py/tests/test_fp_vector_slices.py b/ext_py/tests/test_fp_vector_slices.py
index 31a6782880..39085ce7fb 100644
--- a/ext_py/tests/test_fp_vector_slices.py
+++ b/ext_py/tests/test_fp_vector_slices.py
@@ -6,12 +6,11 @@
 def test_fp_vector_slice_queries_and_to_owned():
     v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
 
-    s = v.slice(1, 4)
+    s = v.const[1:4]
 
     assert s.prime == 5
     assert len(s) == 3
     assert not s.is_empty
-    assert s.entry(1) == 2
     assert s[1] == 2
     assert s[-1] == 0
     assert not s.is_zero
@@ -20,37 +19,16 @@ def test_fp_vector_slice_queries_and_to_owned():
     assert s.iter_nonzero() == [(0, 1), (1, 2)]
     assert repr(s) == "FpSlice(5, [1, 2, 0])"
 
-    restricted = s.restrict(1, 3)
-    assert len(restricted) == 2
-    assert [restricted[i] for i in range(len(restricted))] == [2, 0]
-
     owned = s.to_owned()
     assert isinstance(owned, fp.FpVector)
     assert repr(owned) == "FpVector(5, [1, 2, 0])"
 
 
-def test_fp_vector_restrict_returns_subslice():
-    v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
-
-    r = v.restrict(1, 4)
-    assert isinstance(r, fp.FpSlice)
-    assert len(r) == 3
-    assert r.iter_nonzero() == [(0, 1), (1, 2)]
-    # Mirrors slice(start, end) for an FpVector.
-    assert r.iter_nonzero() == v.slice(1, 4).iter_nonzero()
-    assert repr(r.to_owned()) == "FpVector(5, [1, 2, 0])"
-
-    with pytest.raises(IndexError):
-        v.restrict(0, 6)
-    with pytest.raises(IndexError):
-        v.restrict(3, 2)
-
-
 def test_fp_vector_iter_nonzero():
     v = fp.FpVector.from_py(5, [0, 1, 7, 0, 4])
     assert v.iter_nonzero() == [(1, 1), (2, 2), (4, 4)]
     # Matches the full-vector slice.
-    assert v.iter_nonzero() == v.slice(0, len(v)).iter_nonzero()
+    assert v.iter_nonzero() == v.const[0:len(v)].iter_nonzero()
 
 
 def test_fp_vector_to_owned_is_independent_clone():
@@ -59,34 +37,26 @@ def test_fp_vector_to_owned_is_independent_clone():
     assert isinstance(owned, fp.FpVector)
     assert list(owned) == list(v)
     # Mutating the source does not affect the clone.
-    v.set_entry(0, 3)
+    v[0] = 3
     assert owned[0] == 0
 
 
 def test_fp_vector_slice_range_and_index_errors():
     v = fp.FpVector(3, 4)
-    s = v.slice(1, 3)
+    s = v.const[1:3]
 
     with pytest.raises(IndexError):
-        v.slice(3, 2)
-    with pytest.raises(IndexError):
-        v.slice(0, 5)
-    with pytest.raises(IndexError):
-        v.slice_mut(0, 5)
-    with pytest.raises(IndexError):
-        s.entry(2)
+        v.const[3:2]
     with pytest.raises(IndexError):
         s[2]
     with pytest.raises(IndexError):
         s[-3]
-    with pytest.raises(IndexError):
-        s.restrict(1, 3)
 
 
 def test_fp_slice_handles_parent_lifetime_after_original_reference_deleted():
     v = fp.FpVector.from_py(5, [1, 2, 3, 4])
-    s = v.slice(1, 3)
-    sm = v.slice_mut(2, 4)
+    s = v.const[1:3]
+    sm = v.mut[2:4]
 
     del v
 
@@ -98,8 +68,8 @@ def test_fp_slice_handles_parent_lifetime_after_original_reference_deleted():
 
 def test_stale_slice_handles_raise_python_exception_after_parent_shrink():
     v = fp.FpVector.from_py(5, [1, 2, 3, 4])
-    s = v.slice(1, 4)
-    sm = v.slice_mut(1, 4)
+    s = v.const[1:4]
+    sm = v.mut[1:4]
 
     v.set_scratch_vector_size(2)
 
@@ -112,14 +82,14 @@ def test_stale_slice_handles_raise_python_exception_after_parent_shrink():
     with pytest.raises(IndexError):
         sm[0]
     with pytest.raises(IndexError):
-        sm.set_entry(0, 1)
+        sm[0] = 1
     with pytest.raises(IndexError):
         sm.set_to_zero()
 
 
 def test_fp_slice_mut_updates_parent_and_as_slice():
     v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0])
-    s = v.slice_mut(1, 4)
+    s = v.mut[1:4]
 
     assert s.prime == 5
     assert len(s) == 3
@@ -127,7 +97,7 @@ def test_fp_slice_mut_updates_parent_and_as_slice():
     assert s[-1] == 4
     assert repr(s) == "FpSliceMut(5, [2, 3, 4])"
 
-    s.set_entry(0, 7)
+    s[0] = 7
     assert v[1] == 2
 
     s[1] = 9
@@ -142,32 +112,32 @@ def test_fp_slice_mut_updates_parent_and_as_slice():
     assert len(s) == 3
     assert [s[i] for i in range(len(s))] == [4, 3, 4]
 
-    sub = s.slice_mut(1, 3)
+    sub = s[1:3]
     sub.set_to_zero()
     assert [v[i] for i in range(len(v))] == [1, 4, 0, 0, 0]
 
 
 def test_fp_slice_mut_add_assign_and_masked_methods():
     v = fp.FpVector.from_py(5, [1, 0, 0, 0, 0, 0])
-    target = v.slice_mut(1, 4)
+    target = v.mut[1:4]
     source_v = fp.FpVector.from_py(5, [1, 2, 3, 4])
-    source = source_v.slice(0, 3)
+    source = source_v.const[0:3]
 
     target.add(source, 2)
     assert [v[i] for i in range(len(v))] == [1, 2, 4, 1, 0, 0]
 
-    target.assign(fp.FpVector.from_py(5, [4, 0, 1]).slice(0, 3))
+    target.assign(fp.FpVector.from_py(5, [4, 0, 1]).const[0:3])
     assert [v[i] for i in range(len(v))] == [1, 4, 0, 1, 0, 0]
 
 
 def test_fp_slice_mut_add_tensor():
     v = fp.FpVector(5, 6)
 
-    v.slice_mut(0, 6).add_tensor(
+    v.mut.add_tensor(
         1,
         2,
-        fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
-        fp.FpVector.from_py(5, [3, 4]).slice(0, 2),
+        fp.FpVector.from_py(5, [1, 2]).const,
+        fp.FpVector.from_py(5, [3, 4]).const,
     )
 
     assert [v[i] for i in range(len(v))] == [0, 1, 3, 2, 1, 0]
@@ -175,12 +145,12 @@ def test_fp_slice_mut_add_tensor():
 def test_fp_slice_mut_assign_distinct_and_aliased():
     vt = fp.FpVector.from_py(5, [1, 2, 3])
     vo = fp.FpVector.from_py(5, [4, 0, 1])
-    vt.slice_mut(0, 3).assign(vo.slice(0, 3))
+    vt.mut.assign(vo.const)
     assert [vt[i] for i in range(3)] == [4, 0, 1]
 
     # Aliased clone fallback over the same vector (non-overlapping halves).
     v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0, 1])
-    v.slice_mut(0, 3).assign(v.slice(3, 6))
+    v.mut[0:3].assign(v.const[3:6])
     assert [v[i] for i in range(6)] == [4, 0, 1, 4, 0, 1]
 
 
@@ -189,21 +159,21 @@ def test_fp_slice_mut_add_tensor_distinct_and_aliased():
     vt = fp.FpVector(5, 6)
     left = fp.FpVector.from_py(5, [1, 2])
     right = fp.FpVector.from_py(5, [3, 4])
-    vt.slice_mut(0, 6).add_tensor(1, 2, left.slice(0, 2), right.slice(0, 2))
+    vt.mut.add_tensor(1, 2, left.const, right.const)
     assert [vt[i] for i in range(6)] == [0, 1, 3, 2, 1, 0]
 
     # Dual-operand left == right over the SAME vector (two shared borrows of
     # one operand coexist; neither aliases the distinct target).
     vt2 = fp.FpVector(5, 6)
     w = fp.FpVector.from_py(5, [1, 2])
-    vt2.slice_mut(0, 6).add_tensor(1, 2, w.slice(0, 2), w.slice(0, 2))
+    vt2.mut.add_tensor(1, 2, w.const, w.const)
     assert [vt2[i] for i in range(6)] == [0, 2, 4, 4, 3, 0]
 
     # An operand aliasing the target: the clone fallback must read the operand's
     # ORIGINAL contents even though the target overwrites them mid-op.
     v = fp.FpVector.from_py(5, [1, 2, 0, 0, 0, 0])
     r = fp.FpVector.from_py(5, [3, 4])
-    v.slice_mut(0, 6).add_tensor(1, 2, v.slice(0, 2), r.slice(0, 2))
+    v.mut.add_tensor(1, 2, v.const[0:2], r.const[0:2])
     assert [v[i] for i in range(6)] == [1, 3, 3, 2, 1, 0]
 
 
@@ -225,7 +195,7 @@ def test_fp_slice_mut_accepts_fp_vector_operand():
     # The operand-taking FpSliceMut methods accept either an FpSlice or a
     # (full) FpVector, coercing the FpVector via its full-vector slice.
     v = fp.FpVector.from_py(5, [1, 0, 0, 0, 0, 0])
-    target = v.slice_mut(1, 4)
+    target = v.mut[1:4]
     operand = fp.FpVector.from_py(5, [1, 2, 3])
 
     target.add(operand, 2)
@@ -237,18 +207,17 @@ def test_fp_slice_mut_accepts_fp_vector_operand():
     # An FpVector operand whose backing object IS the target's parent takes the
     # clone fallback and must still be correct.
     w = fp.FpVector.from_py(5, [1, 2, 3])
-    w.slice_mut(0, 3).add(w, 1)
+    w.mut[0:3].add(w, 1)
     assert [w[i] for i in range(len(w))] == [2, 4, 1]
 
     # add_tensor accepts FpVector operands too.
     out = fp.FpVector(5, 6)
     left = fp.FpVector.from_py(5, [1, 2])
     right = fp.FpVector.from_py(5, [3, 4])
-    out.slice_mut(0, 6).add_tensor(1, 2, left, right)
+    out.mut.add_tensor(1, 2, left, right)
     assert [out[i] for i in range(6)] == [0, 1, 3, 2, 1, 0]
 
-    # A matrix row_mut accepting an FpVector operand (the secondary_massey.py
-    # pattern: out.slice_mut(...).add(g, 1) where g is an FpVector).
+    # A matrix row_mut accepting an FpVector operand
     m = fp.Matrix.from_py(5, [[1, 1, 1], [2, 2, 2]])
     m.row_mut(0).add(fp.FpVector.from_py(5, [3, 0, 1]), 1)
     assert m.to_py() == [[4, 1, 2], [2, 2, 2]]
@@ -256,9 +225,9 @@ def test_fp_slice_mut_accepts_fp_vector_operand():
 
 def test_fp_slice_mut_new_method_errors_are_python_exceptions():
     target_v = fp.FpVector(5, 3)
-    target = target_v.slice_mut(0, 3)
-    same_prime_short = fp.FpVector(5, 2).slice(0, 2)
-    other_prime = fp.FpVector(7, 3).slice(0, 3)
+    target = target_v.mut
+    same_prime_short = fp.FpVector(5, 2).const
+    other_prime = fp.FpVector(7, 3).const
 
     with pytest.raises(ValueError):
         target.add(same_prime_short, 1)
@@ -268,14 +237,14 @@ def test_fp_slice_mut_new_method_errors_are_python_exceptions():
         target.add_tensor(
             1,
             1,
-            fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
-            fp.FpVector.from_py(5, [1, 2]).slice(0, 2),
+            fp.FpVector.from_py(5, [1, 2]).const,
+            fp.FpVector.from_py(5, [1, 2]).const,
         )
 
 
 def test_fp_slice_mut_to_owned():
     v = fp.FpVector.from_py(5, [1, 2, 3, 4, 0])
-    s = v.slice_mut(1, 4)
+    s = v.mut[1:4]
 
     owned = s.to_owned()
     assert isinstance(owned, fp.FpVector)
@@ -283,14 +252,14 @@ def test_fp_slice_mut_to_owned():
     assert repr(owned) == "FpVector(5, [2, 3, 4])"
 
     # to_owned is a copy: mutating the source does not change the owned vector.
-    s.set_entry(0, 0)
+    s[0] = 0
     assert repr(owned) == "FpVector(5, [2, 3, 4])"
 
 
 def test_fp_slice_len_revalidates_after_parent_shrink():
     v = fp.FpVector.from_py(5, [1, 2, 3, 4])
-    s = v.slice(1, 4)
-    sm = v.slice_mut(1, 4)
+    s = v.const[1:4]
+    sm = v.mut[1:4]
 
     assert len(s) == 3
     assert len(sm) == 3
@@ -310,10 +279,10 @@ def test_fp_slice_len_revalidates_after_parent_shrink():
 
 def test_fp_slice_mut_index_and_range_errors():
     v = fp.FpVector(2, 3)
-    s = v.slice_mut(1, 3)
+    s = v.mut[1:3]
 
     with pytest.raises(IndexError):
-        s.set_entry(2, 1)
+        s[2] = 1
     with pytest.raises(IndexError):
         s.add_basis_element(2, 1)
     with pytest.raises(IndexError):
@@ -322,7 +291,3 @@ def test_fp_slice_mut_index_and_range_errors():
         s[-3]
     with pytest.raises(IndexError):
         s[2] = 1
-    with pytest.raises(IndexError):
-        s.slice_mut(2, 1)
-    with pytest.raises(IndexError):
-        s.slice_mut(0, 3)
diff --git a/ext_py/tests/test_matrix.py b/ext_py/tests/test_matrix.py
index 77cdf25995..816a3cd650 100644
--- a/ext_py/tests/test_matrix.py
+++ b/ext_py/tests/test_matrix.py
@@ -77,7 +77,7 @@ def test_row_access_and_getitem():
     row = m.row(1)
     assert row.prime == 5
     assert len(row) == 3
-    assert row.entry(0) == 4
+    assert row[0] == 4
     assert row[2] == 1
     assert row[-1] == 1
     assert not row.is_zero
@@ -89,8 +89,6 @@ def test_row_access_and_getitem():
 
     with pytest.raises(IndexError):
         m.row(2)
-    with pytest.raises(IndexError):
-        row.entry(3)
 
 
 def test_row_mut_reflects_in_parent():
@@ -114,7 +112,7 @@ def test_row_mut_reflects_in_parent():
 def test_row_mut_add_slice():
     m = fp.Matrix.from_py(5, [[1, 2, 3]])
     other = fp.FpVector.from_py(5, [1, 1, 1])
-    m.row_mut(0).add(other.slice(0, 3), 2)
+    m.row_mut(0).add(other.const, 2)
     assert m.to_py()[0] == [3, 4, 0]
 
 
@@ -204,11 +202,9 @@ def test_bytes_roundtrip():
 def test_stale_row_handle_after_trim_raises():
     m = fp.Matrix.from_py(5, [[1, 2, 3], [4, 0, 1], [2, 2, 2]])
     row = m.row(2)
-    assert row.entry(0) == 2
+    assert row[0] == 2
     m.trim(0, 1, 0)
     assert m.rows == 1
-    with pytest.raises(IndexError):
-        row.entry(0)
     with pytest.raises(IndexError):
         row[0]
 
@@ -225,9 +221,9 @@ def test_stale_row_mut_handle_after_trim_raises():
 def test_row_returns_same_type_as_vector_slice():
     m = fp.Matrix.from_py(5, [[1, 2, 3]])
     v = fp.FpVector.from_py(5, [1, 2, 3])
-    assert type(m.row(0)) is type(v.slice(0, 3))
-    assert type(m.row_mut(0)) is type(v.slice_mut(0, 3))
-    assert type(m[0]) is type(v.slice(0, 3))
+    assert type(m.row(0)) is type(v.const)
+    assert type(m.row_mut(0)) is type(v.mut)
+    assert type(m[0]) is type(v.const)
 
 
 def test_row_slice_restrict_and_to_owned():
@@ -247,7 +243,7 @@ def test_row_mut_to_owned_and_slice_mut():
     rm = m.row_mut(0)
     owned = rm.to_owned()
     assert owned.prime == 5
-    sub = rm.slice_mut(0, 2)
+    sub = rm[0:2]
     sub.scale(2)
     assert m.to_py()[0] == [2, 4, 3, 4]
     assert repr(rm).startswith("FpSliceMut(5, ")
@@ -268,8 +264,6 @@ def test_row_len_revalidates_after_column_shrink():
 
     with pytest.raises(IndexError):
         len(row)
-    with pytest.raises(IndexError):
-        row.entry(0)
     with pytest.raises(IndexError):
         row[0]
     with pytest.raises(IndexError):
diff --git a/ext_py/tests/test_quasi_inverse.py b/ext_py/tests/test_quasi_inverse.py
index 622b429f7a..a5a4760c8b 100644
--- a/ext_py/tests/test_quasi_inverse.py
+++ b/ext_py/tests/test_quasi_inverse.py
@@ -49,7 +49,7 @@ def test_apply_accepts_slices():
     qi = make_qi()
     v = fp.FpVector.from_py(2, [1, 1, 0, 0, 1, 0])
     out = fp.FpVector(2, 4)
-    qi.apply(out.slice_mut(0, 4), 1, v.slice(0, 6))
+    qi.apply(out.mut, 1, v.const)
     assert list(out) == [1, 1, 1, 0]
 
 
diff --git a/ext_py/tests/test_subquotient.py b/ext_py/tests/test_subquotient.py
index 7bceb13a78..d0aa5332a4 100644
--- a/ext_py/tests/test_subquotient.py
+++ b/ext_py/tests/test_subquotient.py
@@ -173,7 +173,7 @@ def test_reduce_by_quotient_slice_mut():
     # Sub-range slice: only the targeted columns are reduced in place.
     m2 = fp.Matrix.from_py(3, [[2, 1, 1, 0]])
     row = m2.row_mut(0)
-    sq.reduce_by_quotient(row.slice_mut(1, 4))
+    sq.reduce_by_quotient(row[1:4])
     assert list(m2.row(0)) == [2, 0, 1, 0]
 
 
diff --git a/ext_py/tests/test_subspace.py b/ext_py/tests/test_subspace.py
index b564efc070..62ca59d8fa 100644
--- a/ext_py/tests/test_subspace.py
+++ b/ext_py/tests/test_subspace.py
@@ -111,12 +111,12 @@ def test_contains_accepts_slice():
     v = fp.FpVector.from_py(3, [1, 0, 0])
     # An owned FpVector and an FpSlice over it behave identically.
     assert s.contains(v)
-    assert s.contains(v.slice(0, 3))
+    assert s.contains(v.const)
     w = fp.FpVector.from_py(3, [0, 0, 1])
-    assert not s.contains(w.slice(0, 3))
+    assert not s.contains(w.const)
     # Slice of the wrong ambient dimension is rejected.
     with pytest.raises(ValueError):
-        s.contains(w.slice(0, 2))
+        s.contains(w.const[0:2])
 
 
 def test_iter_all_vectors():