diff --git a/crates/pecos-qec/src/distance.rs b/crates/pecos-qec/src/distance.rs
index 511deb4d8..a954d69a2 100644
--- a/crates/pecos-qec/src/distance.rs
+++ b/crates/pecos-qec/src/distance.rs
@@ -344,19 +344,33 @@ pub fn find_min_weight_logicals(
 pub fn find_min_weight_logicals_with_info(
     code: &StabilizerCodeSpec,
     config: &DistanceSearchConfig,
+) -> Vec<LogicalOperatorInfo> {
+    find_shortest_logicals(code, config, 0)
+}
+
+/// Find all logical operators from the minimum weight through `delta` weights above it.
+///
+/// The search always starts at weight 1. Once the minimum logical weight is found,
+/// collection continues through `minimum_weight + delta`, subject to
+/// `config.max_weight`.
+#[must_use]
+pub fn find_shortest_logicals(
+    code: &StabilizerCodeSpec,
+    config: &DistanceSearchConfig,
+    delta: usize,
 ) -> Vec<LogicalOperatorInfo> {
     let max_weight = config.max_weight.unwrap_or(code.num_qubits());
     let mut results = Vec::new();
-    let mut found_distance = None;
+    let mut found_distance: Option<usize> = None;
 
     // Build indices once for O(weight) lookups instead of O(num_stabilizers * weight)
     let stab_index = code.build_stabilizer_index();
     let log_index = code.build_logical_index();
 
     for weight in 1..=max_weight {
-        // If we've found logical operators and this weight is larger, stop
+        // If we've searched through the requested range above the minimum, stop.
         if let Some(d) = found_distance
-            && weight > d
+            && weight > d.saturating_add(delta)
         {
             break;
         }
@@ -437,6 +451,39 @@ mod tests {
         )
     }
 
+    fn five_qubit_code() -> StabilizerCodeSpec {
+        // The [[5,1,3]] perfect code
+        // Stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
+        let stab1 = pauli_string(&[(Pauli::X, 0), (Pauli::Z, 1), (Pauli::Z, 2), (Pauli::X, 3)]);
+        let stab2 = pauli_string(&[(Pauli::X, 1), (Pauli::Z, 2), (Pauli::Z, 3), (Pauli::X, 4)]);
+        let stab3 = pauli_string(&[(Pauli::X, 0), (Pauli::X, 2), (Pauli::Z, 3), (Pauli::Z, 4)]);
+        let stab4 = pauli_string(&[(Pauli::Z, 0), (Pauli::X, 1), (Pauli::X, 3), (Pauli::Z, 4)]);
+
+        // Logical operators for [[5,1,3]]: Z = ZZZZZ, X = XXXXX
+        let logical_z = pauli_string(&[
+            (Pauli::Z, 0),
+            (Pauli::Z, 1),
+            (Pauli::Z, 2),
+            (Pauli::Z, 3),
+            (Pauli::Z, 4),
+        ]);
+        let logical_x = pauli_string(&[
+            (Pauli::X, 0),
+            (Pauli::X, 1),
+            (Pauli::X, 2),
+            (Pauli::X, 3),
+            (Pauli::X, 4),
+        ]);
+
+        StabilizerCodeSpec::new(
+            5,
+            vec![stab1, stab2, stab3, stab4],
+            vec![logical_z],
+            vec![logical_x],
+        )
+        .unwrap()
+    }
+
     #[test]
     fn test_weighted_pauli_iterator_weight_1() {
         let iter = WeightedPauliIterator::new(3, 1, false);
@@ -533,36 +580,7 @@ mod tests {
 
     #[test]
     fn test_five_qubit_code_distance() {
-        // The [[5,1,3]] perfect code
-        // Stabilizers: XZZXI, IXZZX, XIXZZ, ZXIXZ
-        let stab1 = pauli_string(&[(Pauli::X, 0), (Pauli::Z, 1), (Pauli::Z, 2), (Pauli::X, 3)]);
-        let stab2 = pauli_string(&[(Pauli::X, 1), (Pauli::Z, 2), (Pauli::Z, 3), (Pauli::X, 4)]);
-        let stab3 = pauli_string(&[(Pauli::X, 0), (Pauli::X, 2), (Pauli::Z, 3), (Pauli::Z, 4)]);
-        let stab4 = pauli_string(&[(Pauli::Z, 0), (Pauli::X, 1), (Pauli::X, 3), (Pauli::Z, 4)]);
-
-        // Logical operators for [[5,1,3]]: Z = ZZZZZ, X = XXXXX
-        let logical_z = pauli_string(&[
-            (Pauli::Z, 0),
-            (Pauli::Z, 1),
-            (Pauli::Z, 2),
-            (Pauli::Z, 3),
-            (Pauli::Z, 4),
-        ]);
-        let logical_x = pauli_string(&[
-            (Pauli::X, 0),
-            (Pauli::X, 1),
-            (Pauli::X, 2),
-            (Pauli::X, 3),
-            (Pauli::X, 4),
-        ]);
-
-        let code = StabilizerCodeSpec::new(
-            5,
-            vec![stab1, stab2, stab3, stab4],
-            vec![logical_z],
-            vec![logical_x],
-        )
-        .unwrap();
+        let code = five_qubit_code();
 
         // Verify the code is valid
         assert!(code.verify().is_ok());
@@ -576,6 +594,36 @@ mod tests {
         assert_eq!(result.distance, 3);
     }
 
+    #[test]
+    fn test_five_qubit_shortest_logicals_respect_logical_weight_spectrum() {
+        let code = five_qubit_code();
+        let config = DistanceSearchConfig::default();
+        let minimum = find_min_weight_logicals_with_info(&code, &config);
+        let delta_one = find_shortest_logicals(&code, &config, 1);
+        let delta_two = find_shortest_logicals(&code, &config, 2);
+
+        assert_eq!(minimum.len(), 30);
+        assert_eq!(delta_one.len(), 30);
+        assert!(
+            delta_one
+                .iter()
+                .map(|info| &info.operator)
+                .eq(minimum.iter().map(|info| &info.operator))
+        );
+
+        assert_eq!(delta_two.len(), 48);
+        assert_eq!(delta_two.iter().filter(|info| info.weight == 3).count(), 30);
+        assert_eq!(delta_two.iter().filter(|info| info.weight == 5).count(), 18);
+        assert!(delta_two.iter().all(|info| matches!(info.weight, 3 | 5)));
+        assert!(
+            delta_two
+                .iter()
+                .take(minimum.len())
+                .map(|info| &info.operator)
+                .eq(minimum.iter().map(|info| &info.operator))
+        );
+    }
+
     #[test]
     fn test_logical_equivalence_tracking() {
         // 3-qubit bit flip code
diff --git a/crates/pecos-qec/src/lib.rs b/crates/pecos-qec/src/lib.rs
index 50100c2e2..77823cfcf 100644
--- a/crates/pecos-qec/src/lib.rs
+++ b/crates/pecos-qec/src/lib.rs
@@ -69,16 +69,19 @@ pub mod fault_tolerance;
 pub mod geometry;
 pub mod logical_discovery;
 pub mod mem_stab;
+pub mod parity_check_matrix;
 pub mod stabilizer_code;
 pub mod stabilizer_code_spec;
 pub mod surface;
 
 pub use dem_stab::{DemStabError, DemStabShotBatch, DemStabSim, DemStabSimBuilder};
 pub use mem_stab::{MemStabError, MemStabSim, MemStabSimBuilder};
+pub use parity_check_matrix::{ParityCheckMatrix, ParityCheckMatrixError};
 
 pub use distance::{
     DistanceResult, DistanceSearchConfig, LogicalOperatorInfo, WeightedPauliIterator,
     calculate_distance, find_min_weight_logicals, find_min_weight_logicals_with_info,
+    find_shortest_logicals,
 };
 pub use fault_tolerance::dem_builder::{
     DecomposedFault, DemBuilder, DemBuilderError, DemOutput, DetectorDef, DetectorErrorModel,
diff --git a/crates/pecos-qec/src/parity_check_matrix.rs b/crates/pecos-qec/src/parity_check_matrix.rs
new file mode 100644
index 000000000..2b02a379c
--- /dev/null
+++ b/crates/pecos-qec/src/parity_check_matrix.rs
@@ -0,0 +1,188 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+//! Role-neutral binary parity-check matrices for QEC codes.
+
+use pecos_core::{Pauli, PauliString, QuarterPhase, QubitId};
+use pecos_quantum::F2Matrix;
+use thiserror::Error;
+
+/// Errors that can occur when constructing a [`ParityCheckMatrix`].
+#[derive(Debug, Clone, PartialEq, Eq, Error)]
+pub enum ParityCheckMatrixError {
+    /// No rows were supplied, so the matrix width cannot be inferred.
+    #[error(
+        "cannot infer parity-check matrix width from empty input; use ParityCheckMatrix::zeros"
+    )]
+    EmptyRows,
+    /// A row has a different width from the first row.
+    #[error("parity-check matrix row {row} has {actual} columns, expected {expected}")]
+    RaggedRows {
+        /// Index of the mismatched row.
+        row: usize,
+        /// Width inferred from the first row.
+        expected: usize,
+        /// Actual width of the mismatched row.
+        actual: usize,
+    },
+    /// A dense entry was not binary.
+    #[error("parity-check matrix entry at row {row}, column {column} is {value}, expected 0 or 1")]
+    InvalidEntry {
+        /// Row containing the invalid entry.
+        row: usize,
+        /// Column containing the invalid entry.
+        column: usize,
+        /// Invalid value.
+        value: u8,
+    },
+}
+
+/// A role-neutral binary matrix whose rows are checks and columns are qubits.
+///
+/// Whether rows become X-type or Z-type stabilizers is chosen only when
+/// converting the matrix; that role is not stored in this type.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub struct ParityCheckMatrix {
+    matrix: F2Matrix,
+}
+
+impl ParityCheckMatrix {
+    /// Constructs a validated parity-check matrix from dense binary rows.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error for empty input, ragged rows, or non-binary entries.
+    pub fn from_dense(rows: Vec<Vec<u8>>) -> Result<Self, ParityCheckMatrixError> {
+        let Some(first) = rows.first() else {
+            return Err(ParityCheckMatrixError::EmptyRows);
+        };
+        let num_qubits = first.len();
+        for (row_index, row) in rows.iter().enumerate() {
+            if row.len() != num_qubits {
+                return Err(ParityCheckMatrixError::RaggedRows {
+                    row: row_index,
+                    expected: num_qubits,
+                    actual: row.len(),
+                });
+            }
+            for (column, &value) in row.iter().enumerate() {
+                if value > 1 {
+                    return Err(ParityCheckMatrixError::InvalidEntry {
+                        row: row_index,
+                        column,
+                        value,
+                    });
+                }
+            }
+        }
+        Ok(Self {
+            matrix: F2Matrix::from_rows(rows),
+        })
+    }
+
+    /// Constructs an all-zero matrix with an explicit number of qubits.
+    #[must_use]
+    pub fn zeros(num_checks: usize, num_qubits: usize) -> Self {
+        Self {
+            matrix: F2Matrix::zeros(num_checks, num_qubits),
+        }
+    }
+
+    /// Returns the number of checks (rows).
+    #[must_use]
+    pub fn num_checks(&self) -> usize {
+        self.matrix.num_rows()
+    }
+
+    /// Returns the number of qubits (columns).
+    #[must_use]
+    pub fn num_qubits(&self) -> usize {
+        self.matrix.num_cols()
+    }
+
+    /// Returns the rank over GF(2).
+    #[must_use]
+    pub fn rank(&self) -> usize {
+        self.matrix.row_reduce().1.len()
+    }
+
+    /// Returns dense copies of all rows.
+    #[must_use]
+    pub fn rows(&self) -> Vec<Vec<u8>> {
+        self.matrix.rows()
+    }
+
+    /// Returns a dense copy of one row, or `None` if the index is out of range.
+    #[must_use]
+    pub fn row(&self, index: usize) -> Option<Vec<u8>> {
+        (index < self.num_checks()).then(|| self.matrix.row(index))
+    }
+
+    /// Converts rows to stabilizers made of Pauli X operators, with phase `+1`.
+    ///
+    /// “X stabilizers” means stabilizers made of X, not stabilizers that detect
+    /// X errors.
+    #[must_use]
+    pub fn to_x_stabilizers(&self) -> Vec<PauliString> {
+        self.to_stabilizers(Pauli::X)
+    }
+
+    /// Converts rows to stabilizers made of Pauli Z operators, with phase `+1`.
+    ///
+    /// “Z stabilizers” means stabilizers made of Z, not stabilizers that detect
+    /// Z errors.
+    #[must_use]
+    pub fn to_z_stabilizers(&self) -> Vec<PauliString> {
+        self.to_stabilizers(Pauli::Z)
+    }
+
+    pub(crate) fn matrix(&self) -> &F2Matrix {
+        &self.matrix
+    }
+
+    fn to_stabilizers(&self, pauli: Pauli) -> Vec<PauliString> {
+        (0..self.num_checks())
+            .map(|row| {
+                let paulis = (0..self.num_qubits())
+                    .filter(|&qubit| self.matrix.get(row, qubit) == 1)
+                    .map(|qubit| (pauli, QubitId::new(qubit)))
+                    .collect();
+                PauliString::with_phase_and_paulis(QuarterPhase::PlusOne, paulis)
+            })
+            .collect()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn rank_detects_dependent_rows() {
+        let matrix =
+            ParityCheckMatrix::from_dense(vec![vec![1, 1, 0], vec![0, 1, 1], vec![1, 0, 1]])
+                .unwrap();
+
+        assert_eq!(matrix.num_checks(), 3);
+        assert_eq!(matrix.rank(), 2);
+    }
+
+    #[test]
+    fn zero_row_matrix_preserves_width() {
+        let matrix = ParityCheckMatrix::zeros(0, 9);
+        assert_eq!(matrix.num_checks(), 0);
+        assert_eq!(matrix.num_qubits(), 9);
+        assert!(matrix.rows().is_empty());
+    }
+}
diff --git a/crates/pecos-qec/src/stabilizer_code_spec.rs b/crates/pecos-qec/src/stabilizer_code_spec.rs
index e1d753106..44ce78f04 100644
--- a/crates/pecos-qec/src/stabilizer_code_spec.rs
+++ b/crates/pecos-qec/src/stabilizer_code_spec.rs
@@ -17,7 +17,9 @@
 // Allow similar names for logical_xs/logical_zs - these are intentional and meaningful
 #![allow(clippy::similar_names)]
 
+use crate::parity_check_matrix::ParityCheckMatrix;
 use pecos_core::{PauliOperator, PauliString};
+use pecos_quantum::{PauliSequence, SymplecticMatrix};
 use std::collections::BTreeSet;
 use thiserror::Error;
 
@@ -44,6 +46,22 @@ pub enum StabilizerCodeSpecError {
     #[error("Logical X{0} and Z{1} anticommute (should commute for different logical qubits)")]
     CrossLogicalAnticommute(usize, usize),
 
+    /// Stabilizer generators are linearly dependent over GF(2).
+    #[error("Stabilizer generators are dependent: rank {rank}, count {count}")]
+    DependentStabilizers { rank: usize, count: usize },
+
+    /// A typed matrix width does not match the builder width.
+    #[error("{matrix} matrix has {actual} qubits, expected {expected}")]
+    MatrixWidthMismatch {
+        matrix: &'static str,
+        expected: usize,
+        actual: usize,
+    },
+
+    /// A CSS X row and Z row are not orthogonal over GF(2).
+    #[error("CSS X row {x_row} and Z row {z_row} are not orthogonal")]
+    CssRowsNotOrthogonal { x_row: usize, z_row: usize },
+
     /// Invalid code parameters.
     #[error("Invalid code: {0}")]
     InvalidCode(String),
@@ -75,6 +93,42 @@ pub struct StabilizerCodeSpec {
     distance: Option<usize>,
 }
 
+impl std::fmt::Display for StabilizerCodeSpec {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        writeln!(f, "[[{}, {}]]", self.num_qubits, self.num_logical_qubits())?;
+
+        writeln!(f, "Stabilizer generators:")?;
+        for stabilizer in &self.stabilizers {
+            writeln!(f, "{}", stabilizer.to_dense_str(Some(self.num_qubits)))?;
+        }
+
+        writeln!(f, "Destabilizer generators:")?;
+        for destabilizer in &self.destabilizers {
+            writeln!(f, "{}", destabilizer.to_dense_str(Some(self.num_qubits)))?;
+        }
+
+        writeln!(f, "Logical operators:")?;
+        for (index, (logical_z, logical_x)) in
+            self.logical_zs.iter().zip(&self.logical_xs).enumerate()
+        {
+            writeln!(
+                f,
+                "Z{}: {}",
+                index + 1,
+                logical_z.to_dense_str(Some(self.num_qubits))
+            )?;
+            writeln!(
+                f,
+                "X{}: {}",
+                index + 1,
+                logical_x.to_dense_str(Some(self.num_qubits))
+            )?;
+        }
+
+        Ok(())
+    }
+}
+
 /// Column-based index for efficient commutation checking.
 ///
 /// For each qubit, tracks which operators have X or Z on that qubit.
@@ -167,7 +221,8 @@ impl StabilizerCodeSpec {
     /// - `logical_xs`: Logical X operators (one per logical qubit)
     ///
     /// # Errors
-    /// Returns an error if the logical X and Z vectors have different lengths.
+    /// Returns an error if the logical X and Z vectors have different lengths,
+    /// or if the stabilizer generators are linearly dependent.
     pub fn new(
         num_qubits: usize,
         stabilizers: Vec<PauliString>,
@@ -180,6 +235,15 @@ impl StabilizerCodeSpec {
             ));
         }
 
+        let stabilizer_count = stabilizers.len();
+        let stabilizer_rank = PauliSequence::new(stabilizers.clone()).rank();
+        if stabilizer_rank != stabilizer_count {
+            return Err(StabilizerCodeSpecError::DependentStabilizers {
+                rank: stabilizer_rank,
+                count: stabilizer_count,
+            });
+        }
+
         Ok(Self {
             num_qubits,
             stabilizers,
@@ -200,7 +264,8 @@ impl StabilizerCodeSpec {
     /// - `logical_xs`: Logical X operators (one per logical qubit)
     ///
     /// # Errors
-    /// Returns an error if the logical X and Z vectors have different lengths.
+    /// Returns an error if the logical X and Z vectors have different lengths,
+    /// or if the stabilizer generators are linearly dependent.
     pub fn with_destabilizers(
         num_qubits: usize,
         stabilizers: Vec<PauliString>,
@@ -208,35 +273,20 @@ impl StabilizerCodeSpec {
         logical_zs: Vec<PauliString>,
         logical_xs: Vec<PauliString>,
     ) -> Result<Self> {
-        if logical_zs.len() != logical_xs.len() {
-            return Err(StabilizerCodeSpecError::InvalidCode(
-                "Number of logical X and Z operators must match".to_string(),
-            ));
-        }
-
-        Ok(Self {
-            num_qubits,
-            stabilizers,
-            destabilizers,
-            logical_zs,
-            logical_xs,
-            distance: None,
-        })
+        let mut code = Self::new(num_qubits, stabilizers, logical_zs, logical_xs)?;
+        code.destabilizers = destabilizers;
+        Ok(code)
     }
 
     /// Creates a stabilizer code from just the stabilizers.
     ///
     /// The logical operators can be added later.
-    #[must_use]
-    pub fn from_stabilizers(num_qubits: usize, stabilizers: Vec<PauliString>) -> Self {
-        Self {
-            num_qubits,
-            stabilizers,
-            destabilizers: Vec::new(),
-            logical_zs: Vec::new(),
-            logical_xs: Vec::new(),
-            distance: None,
-        }
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if the stabilizer generators are linearly dependent.
+    pub fn from_stabilizers(num_qubits: usize, stabilizers: Vec<PauliString>) -> Result<Self> {
+        Self::new(num_qubits, stabilizers, Vec::new(), Vec::new())
     }
 
     /// Creates a builder for constructing a stabilizer code.
@@ -716,7 +766,7 @@ impl StabilizerCodeSpec {
     /// let mut code = StabilizerCodeSpec::from_stabilizers(3, vec![
     ///     Zs([0, 1]),  // ZZI
     ///     Zs([1, 2]),  // IZZ
-    /// ]);
+    /// ]).unwrap();
     ///
     /// // Discover logical operators
     /// code.discover_logicals().unwrap();
@@ -775,16 +825,12 @@ impl StabilizerCodeSpec {
     /// The resulting code has stabilizer generators but no logical operators
     /// or destabilizers. Use [`discover_logicals`](Self::discover_logicals)
     /// to compute them.
-    #[must_use]
-    pub fn from_stabilizer_group(group: &pecos_quantum::PauliStabilizerGroup) -> Self {
-        Self {
-            num_qubits: group.num_qubits(),
-            stabilizers: group.stabilizers().to_vec(),
-            destabilizers: Vec::new(),
-            logical_zs: Vec::new(),
-            logical_xs: Vec::new(),
-            distance: None,
-        }
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if the stabilizer generators are linearly dependent.
+    pub fn from_stabilizer_group(group: &pecos_quantum::PauliStabilizerGroup) -> Result<Self> {
+        Self::from_stabilizers(group.num_qubits(), group.stabilizers().to_vec())
     }
 
     /// Creates a `StabilizerCodeSpec` from a [`StabilizerCode`](crate::StabilizerCode),
@@ -813,7 +859,8 @@ impl StabilizerCodeSpec {
         code: &crate::StabilizerCode,
     ) -> std::result::Result<Self, crate::LogicalDiscoveryError> {
         let mut spec =
-            Self::from_stabilizers(code.num_qubits(), code.group().stabilizers().to_vec());
+            Self::from_stabilizers(code.num_qubits(), code.group().stabilizers().to_vec())
+                .map_err(|_| crate::LogicalDiscoveryError::StabilizersNotIndependent)?;
         spec.discover_logicals()?;
         Ok(spec)
     }
@@ -959,6 +1006,76 @@ impl StabilizerCodeSpecBuilder {
         self
     }
 
+    /// Appends X-type and Z-type stabilizers from role-neutral CSS matrices.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if either matrix width differs from the builder width,
+    /// or if the first non-orthogonal X/Z row pair is found.
+    pub fn checks_from_css(
+        mut self,
+        x_stabilizers: &ParityCheckMatrix,
+        z_stabilizers: &ParityCheckMatrix,
+    ) -> Result<Self> {
+        if x_stabilizers.num_qubits() != self.num_qubits {
+            return Err(StabilizerCodeSpecError::MatrixWidthMismatch {
+                matrix: "X parity-check",
+                expected: self.num_qubits,
+                actual: x_stabilizers.num_qubits(),
+            });
+        }
+        if z_stabilizers.num_qubits() != self.num_qubits {
+            return Err(StabilizerCodeSpecError::MatrixWidthMismatch {
+                matrix: "Z parity-check",
+                expected: self.num_qubits,
+                actual: z_stabilizers.num_qubits(),
+            });
+        }
+
+        let overlaps = x_stabilizers
+            .matrix()
+            .mul(&z_stabilizers.matrix().transpose());
+        for x_row in 0..overlaps.num_rows() {
+            for z_row in 0..overlaps.num_cols() {
+                if overlaps.get(x_row, z_row) == 1 {
+                    return Err(StabilizerCodeSpecError::CssRowsNotOrthogonal { x_row, z_row });
+                }
+            }
+        }
+
+        self.stabilizers.extend(x_stabilizers.to_x_stabilizers());
+        self.stabilizers.extend(z_stabilizers.to_z_stabilizers());
+        Ok(self)
+    }
+
+    /// Appends mutually commuting stabilizers from symplectic rows.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if the matrix width differs from the builder width, or
+    /// if the first anticommuting row pair is found.
+    pub fn checks_from_symplectic(mut self, matrix: &SymplecticMatrix) -> Result<Self> {
+        if matrix.num_qubits() != self.num_qubits {
+            return Err(StabilizerCodeSpecError::MatrixWidthMismatch {
+                matrix: "Symplectic",
+                expected: self.num_qubits,
+                actual: matrix.num_qubits(),
+            });
+        }
+
+        let stabilizers = matrix.to_positive_paulis();
+        for i in 0..stabilizers.len() {
+            for j in (i + 1)..stabilizers.len() {
+                if !stabilizers[i].commutes_with(&stabilizers[j]) {
+                    return Err(StabilizerCodeSpecError::StabilizersAnticommute(i, j));
+                }
+            }
+        }
+
+        self.stabilizers.extend(stabilizers);
+        Ok(self)
+    }
+
     /// Adds a logical Z operator from a `PauliString` directly.
     #[must_use]
     pub fn logical_z_pauli(mut self, pauli: PauliString) -> Self {
@@ -1047,7 +1164,8 @@ impl StabilizerCodeSpecBuilder {
     pub fn build_with_discovered_logicals(
         self,
     ) -> std::result::Result<StabilizerCodeSpec, crate::LogicalDiscoveryError> {
-        let mut code = StabilizerCodeSpec::from_stabilizers(self.num_qubits, self.stabilizers);
+        let mut code = StabilizerCodeSpec::from_stabilizers(self.num_qubits, self.stabilizers)
+            .map_err(|_| crate::LogicalDiscoveryError::StabilizersNotIndependent)?;
         code.discover_logicals()?;
         Ok(code)
     }
@@ -1056,7 +1174,7 @@ impl StabilizerCodeSpecBuilder {
 #[cfg(test)]
 mod tests {
     use super::*;
-    use pecos_core::Pauli;
+    use pecos_core::{Pauli, Zs};
 
     /// Helper to create a `PauliString` from a simple specification.
     fn pauli_string(paulis: &[(Pauli, usize)]) -> PauliString {
@@ -1067,6 +1185,30 @@ mod tests {
         )
     }
 
+    #[test]
+    fn display_includes_generators_and_paired_logicals_in_order() {
+        let code = StabilizerCodeSpecBuilder::new(3)
+            .check(Zs([0, 1]))
+            .check(Zs([1, 2]))
+            .build_with_discovered_logicals()
+            .unwrap();
+
+        let rendered = code.to_string();
+        assert!(rendered.starts_with("[[3, 1]]\nStabilizer generators:\n"));
+
+        let stabilizer_section = rendered.find("Stabilizer generators:").unwrap();
+        let destabilizer_section = rendered.find("Destabilizer generators:").unwrap();
+        let logical_section = rendered.find("Logical operators:").unwrap();
+        assert!(stabilizer_section < destabilizer_section);
+        assert!(destabilizer_section < logical_section);
+
+        for operator in code.stabilizers().iter().chain(code.destabilizers()) {
+            assert!(rendered.contains(&operator.to_dense_str(Some(3))));
+        }
+        assert!(rendered.contains("Z1: "));
+        assert!(rendered.contains("X1: "));
+    }
+
     #[test]
     fn test_three_qubit_bit_flip_code() {
         // 3-qubit bit flip code: [[3, 1, 1]]
@@ -1114,7 +1256,7 @@ mod tests {
         let stab1 = pauli_string(&[(Pauli::X, 0)]);
         let stab2 = pauli_string(&[(Pauli::Z, 0)]);
 
-        let code = StabilizerCodeSpec::from_stabilizers(1, vec![stab1, stab2]);
+        let code = StabilizerCodeSpec::from_stabilizers(1, vec![stab1, stab2]).unwrap();
 
         let result = code.verify_stabilizers_commute();
         assert!(matches!(
@@ -1123,6 +1265,36 @@ mod tests {
         ));
     }
 
+    #[test]
+    fn new_rejects_dependent_stabilizers() {
+        let result = StabilizerCodeSpec::new(
+            3,
+            vec![Zs([0, 1]), Zs([1, 2]), Zs([0, 2])],
+            Vec::new(),
+            Vec::new(),
+        );
+
+        assert!(matches!(
+            &result,
+            Err(StabilizerCodeSpecError::DependentStabilizers { rank: 2, count: 3 })
+        ));
+        assert_eq!(
+            result.unwrap_err().to_string(),
+            "Stabilizer generators are dependent: rank 2, count 3"
+        );
+    }
+
+    #[test]
+    fn from_stabilizers_rejects_dependent_stabilizers() {
+        let result =
+            StabilizerCodeSpec::from_stabilizers(3, vec![Zs([0, 1]), Zs([1, 2]), Zs([0, 2])]);
+
+        assert!(matches!(
+            result,
+            Err(StabilizerCodeSpecError::DependentStabilizers { rank: 2, count: 3 })
+        ));
+    }
+
     #[test]
     fn test_logical_pair_must_anticommute() {
         // Create a code where logical X and Z commute (invalid)
@@ -1145,7 +1317,7 @@ mod tests {
         // 3-qubit bit flip code
         let stab1 = pauli_string(&[(Pauli::Z, 0), (Pauli::Z, 1)]);
         let stab2 = pauli_string(&[(Pauli::Z, 1), (Pauli::Z, 2)]);
-        let code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]);
+        let code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]).unwrap();
 
         // X error on qubit 0 should trigger stabilizer 0 only
         let x0 = pauli_string(&[(Pauli::X, 0)]);
@@ -1168,7 +1340,7 @@ mod tests {
     fn test_code_parameters_string() {
         let stab1 = pauli_string(&[(Pauli::Z, 0), (Pauli::Z, 1)]);
         let stab2 = pauli_string(&[(Pauli::Z, 1), (Pauli::Z, 2)]);
-        let mut code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]);
+        let mut code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]).unwrap();
 
         assert_eq!(code.code_parameters(), "[[3, 1, ?]]");
 
@@ -1260,7 +1432,7 @@ mod tests {
         // 3-qubit bit flip code
         let stab1 = pauli_string(&[(Pauli::Z, 0), (Pauli::Z, 1)]);
         let stab2 = pauli_string(&[(Pauli::Z, 1), (Pauli::Z, 2)]);
-        let code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]);
+        let code = StabilizerCodeSpec::from_stabilizers(3, vec![stab1, stab2]).unwrap();
         let index = code.build_stabilizer_index();
 
         // X error on qubit 0 should trigger stabilizer 0 only
@@ -1441,6 +1613,81 @@ mod tests {
     // Builder tests
     // ========================================================================
 
+    #[test]
+    fn builder_css_checks_validate_orthogonality() {
+        let x = ParityCheckMatrix::from_dense(vec![vec![1, 0]]).unwrap();
+        let z = ParityCheckMatrix::from_dense(vec![vec![1, 0]]).unwrap();
+
+        let result = StabilizerCodeSpecBuilder::new(2).checks_from_css(&x, &z);
+        assert!(matches!(
+            &result,
+            Err(StabilizerCodeSpecError::CssRowsNotOrthogonal { x_row: 0, z_row: 0 })
+        ));
+        assert_eq!(
+            result.unwrap_err().to_string(),
+            "CSS X row 0 and Z row 0 are not orthogonal"
+        );
+    }
+
+    #[test]
+    fn builder_css_checks_construct_steane_code() {
+        let h = ParityCheckMatrix::from_dense(vec![
+            vec![1, 0, 1, 0, 1, 0, 1],
+            vec![0, 1, 1, 0, 0, 1, 1],
+            vec![0, 0, 0, 1, 1, 1, 1],
+        ])
+        .unwrap();
+
+        let code = StabilizerCodeSpecBuilder::new(7)
+            .checks_from_css(&h, &h)
+            .unwrap()
+            .build_with_discovered_logicals()
+            .unwrap();
+
+        assert_eq!(code.num_logical_qubits(), 1);
+        assert_eq!(
+            crate::calculate_distance(&code, &crate::DistanceSearchConfig::default())
+                .unwrap()
+                .distance,
+            3
+        );
+    }
+
+    #[test]
+    fn builder_symplectic_checks_construct_five_qubit_code() {
+        let matrix = SymplecticMatrix::from_dense(vec![
+            vec![1, 0, 0, 1, 0, 0, 1, 1, 0, 0],
+            vec![0, 1, 0, 0, 1, 0, 0, 1, 1, 0],
+            vec![1, 0, 1, 0, 0, 0, 0, 0, 1, 1],
+            vec![0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
+        ])
+        .unwrap();
+
+        let code = StabilizerCodeSpecBuilder::new(5)
+            .checks_from_symplectic(&matrix)
+            .unwrap()
+            .build_with_discovered_logicals()
+            .unwrap();
+
+        assert_eq!(
+            crate::calculate_distance(&code, &crate::DistanceSearchConfig::default())
+                .unwrap()
+                .distance,
+            3
+        );
+    }
+
+    #[test]
+    fn builder_symplectic_checks_reject_anticommuting_rows() {
+        let matrix = SymplecticMatrix::from_dense(vec![vec![1, 0], vec![0, 1]]).unwrap();
+
+        let result = StabilizerCodeSpecBuilder::new(1).checks_from_symplectic(&matrix);
+        assert!(matches!(
+            result,
+            Err(StabilizerCodeSpecError::StabilizersAnticommute(0, 1))
+        ));
+    }
+
     #[test]
     fn test_builder_three_qubit_bit_flip() {
         use pecos_core::{Xs, Zs};
@@ -1568,7 +1815,8 @@ mod tests {
                 Zs([0, 1]), // ZZI
                 Zs([1, 2]), // IZZ
             ],
-        );
+        )
+        .unwrap();
 
         assert!(!code.has_logicals());
 
@@ -1631,7 +1879,7 @@ mod tests {
     #[test]
     fn test_from_stabilizer_group() {
         let steane = crate::StabilizerCode::steane();
-        let code = StabilizerCodeSpec::from_stabilizer_group(steane.group());
+        let code = StabilizerCodeSpec::from_stabilizer_group(steane.group()).unwrap();
 
         assert_eq!(code.num_qubits(), 7);
         assert_eq!(code.num_stabilizers(), 6);
@@ -1654,13 +1902,25 @@ mod tests {
             .unwrap();
 
         let group = original.to_stabilizer_group().unwrap();
-        let roundtripped = StabilizerCodeSpec::from_stabilizer_group(&group);
+        let roundtripped = StabilizerCodeSpec::from_stabilizer_group(&group).unwrap();
 
         assert_eq!(roundtripped.num_qubits(), original.num_qubits());
         assert_eq!(roundtripped.num_stabilizers(), original.num_stabilizers());
         assert!(roundtripped.verify_stabilizers_commute().is_ok());
     }
 
+    #[test]
+    fn from_stabilizer_group_rejects_dependent_stabilizers() {
+        let group =
+            pecos_quantum::PauliStabilizerGroup::new(vec![Zs([0, 1]), Zs([1, 2]), Zs([0, 2])])
+                .unwrap();
+
+        assert!(matches!(
+            StabilizerCodeSpec::from_stabilizer_group(&group),
+            Err(StabilizerCodeSpecError::DependentStabilizers { rank: 2, count: 3 })
+        ));
+    }
+
     #[test]
     fn test_stabilizer_group_algebraic_analysis() {
         use pecos_core::pauli::*;
diff --git a/crates/pecos-quantum/src/lib.rs b/crates/pecos-quantum/src/lib.rs
index 3a59ea966..031c70683 100644
--- a/crates/pecos-quantum/src/lib.rs
+++ b/crates/pecos-quantum/src/lib.rs
@@ -74,6 +74,7 @@ pub mod pauli_group;
 pub mod pauli_sequence;
 pub mod pauli_set;
 pub mod stabilizer_group;
+pub mod symplectic_matrix;
 mod tick_circuit;
 pub mod unitary_matrix;
 
@@ -127,6 +128,7 @@ pub use pauli_group::{PauliGroup, PauliGroupError};
 pub use pauli_sequence::{F2Matrix, PauliSequence};
 pub use pauli_set::PauliSet;
 pub use stabilizer_group::{PauliStabilizerGroup, PauliStabilizerGroupError};
+pub use symplectic_matrix::{SymplecticMatrix, SymplecticMatrixError};
 
 // Re-export HUGR types when the feature is enabled
 #[cfg(feature = "hugr")]
diff --git a/crates/pecos-quantum/src/symplectic_matrix.rs b/crates/pecos-quantum/src/symplectic_matrix.rs
new file mode 100644
index 000000000..605865137
--- /dev/null
+++ b/crates/pecos-quantum/src/symplectic_matrix.rs
@@ -0,0 +1,348 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+//! Validated binary symplectic matrices for Pauli operators.
+
+use crate::{F2Matrix, PauliSequence};
+use pecos_core::{Pauli, PauliString, QuarterPhase, QubitId};
+use std::fmt;
+
+/// Errors that can occur when constructing a [`SymplecticMatrix`].
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum SymplecticMatrixError {
+    /// No rows were supplied, so the matrix width cannot be inferred.
+    EmptyRows,
+    /// A row has a different width from the first row.
+    RaggedRows {
+        /// Index of the mismatched row.
+        row: usize,
+        /// Width inferred from the first row.
+        expected: usize,
+        /// Actual width of the mismatched row.
+        actual: usize,
+    },
+    /// A dense entry was not binary.
+    InvalidEntry {
+        /// Row containing the invalid entry.
+        row: usize,
+        /// Column containing the invalid entry.
+        column: usize,
+        /// Invalid value.
+        value: u8,
+    },
+    /// A symplectic matrix must have equally sized X and Z blocks.
+    OddColumnCount {
+        /// Number of supplied columns.
+        columns: usize,
+    },
+    /// A Pauli operator acts beyond the requested explicit width.
+    OperatorExceedsWidth {
+        /// Index of the offending operator.
+        row: usize,
+        /// Offending qubit index.
+        qubit: usize,
+        /// Requested number of qubits.
+        num_qubits: usize,
+    },
+}
+
+impl fmt::Display for SymplecticMatrixError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Self::EmptyRows => write!(
+                f,
+                "cannot infer symplectic matrix width from empty input; use SymplecticMatrix::zeros"
+            ),
+            Self::RaggedRows {
+                row,
+                expected,
+                actual,
+            } => write!(
+                f,
+                "symplectic matrix row {row} has {actual} columns, expected {expected}"
+            ),
+            Self::InvalidEntry { row, column, value } => write!(
+                f,
+                "symplectic matrix entry at row {row}, column {column} is {value}, expected 0 or 1"
+            ),
+            Self::OddColumnCount { columns } => write!(
+                f,
+                "symplectic matrix has {columns} columns, expected an even column count"
+            ),
+            Self::OperatorExceedsWidth {
+                row,
+                qubit,
+                num_qubits,
+            } => write!(
+                f,
+                "Pauli operator at row {row} acts on qubit {qubit}, outside the explicit width of {num_qubits} qubits"
+            ),
+        }
+    }
+}
+
+impl std::error::Error for SymplecticMatrixError {}
+
+/// A binary symplectic matrix whose rows are Pauli operators.
+///
+/// For `n` qubits, columns are ordered as
+/// `[x_0, ..., x_{n-1} | z_0, ..., z_{n-1}]`.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub struct SymplecticMatrix {
+    matrix: F2Matrix,
+}
+
+impl SymplecticMatrix {
+    /// Constructs a validated symplectic matrix from dense binary rows.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error for empty input, ragged rows, non-binary entries, or
+    /// an odd number of columns.
+    pub fn from_dense(rows: Vec<Vec<u8>>) -> Result<Self, SymplecticMatrixError> {
+        let Some(first) = rows.first() else {
+            return Err(SymplecticMatrixError::EmptyRows);
+        };
+        let num_cols = first.len();
+        if num_cols % 2 != 0 {
+            return Err(SymplecticMatrixError::OddColumnCount { columns: num_cols });
+        }
+        for (row_index, row) in rows.iter().enumerate() {
+            if row.len() != num_cols {
+                return Err(SymplecticMatrixError::RaggedRows {
+                    row: row_index,
+                    expected: num_cols,
+                    actual: row.len(),
+                });
+            }
+            for (column, &value) in row.iter().enumerate() {
+                if value > 1 {
+                    return Err(SymplecticMatrixError::InvalidEntry {
+                        row: row_index,
+                        column,
+                        value,
+                    });
+                }
+            }
+        }
+        Ok(Self {
+            matrix: F2Matrix::from_rows(rows),
+        })
+    }
+
+    /// Constructs an all-zero matrix with an explicit number of qubits.
+    #[must_use]
+    pub fn zeros(num_rows: usize, num_qubits: usize) -> Self {
+        Self {
+            matrix: F2Matrix::zeros(num_rows, 2 * num_qubits),
+        }
+    }
+
+    /// Converts a Pauli sequence to symplectic form at an explicit width.
+    ///
+    /// Pauli phases are deliberately ignored. Use [`to_positive_paulis`](Self::to_positive_paulis)
+    /// to recover operators with phase `+1`.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if any operator acts on a qubit outside `num_qubits`.
+    pub fn from_pauli_sequence_ignoring_phase(
+        sequence: &PauliSequence,
+        num_qubits: usize,
+    ) -> Result<Self, SymplecticMatrixError> {
+        for (row, operator) in sequence.iter().enumerate() {
+            if let Some(qubit) = operator
+                .qubits()
+                .into_iter()
+                .find(|&qubit| qubit >= num_qubits)
+            {
+                return Err(SymplecticMatrixError::OperatorExceedsWidth {
+                    row,
+                    qubit,
+                    num_qubits,
+                });
+            }
+        }
+
+        let inferred_num_qubits = sequence.num_qubits();
+        let inferred = sequence.to_symplectic_matrix();
+        let mut matrix = F2Matrix::zeros(sequence.len(), 2 * num_qubits);
+        for row in 0..sequence.len() {
+            for qubit in 0..inferred_num_qubits {
+                matrix.set(row, qubit, inferred.get(row, qubit));
+                matrix.set(
+                    row,
+                    num_qubits + qubit,
+                    inferred.get(row, inferred_num_qubits + qubit),
+                );
+            }
+        }
+        Ok(Self { matrix })
+    }
+
+    /// Converts each row to a phase-`+1` Pauli operator.
+    ///
+    /// Symplectic matrices contain no sign or quarter-phase information, so
+    /// every returned operator necessarily has positive phase.
+    #[must_use]
+    pub fn to_positive_paulis(&self) -> Vec<PauliString> {
+        let num_qubits = self.num_qubits();
+        (0..self.num_rows())
+            .map(|row| {
+                let mut paulis = Vec::new();
+                for qubit in 0..num_qubits {
+                    let x = self.matrix.get(row, qubit);
+                    let z = self.matrix.get(row, num_qubits + qubit);
+                    let pauli = match (x, z) {
+                        (1, 0) => Some(Pauli::X),
+                        (0, 1) => Some(Pauli::Z),
+                        (1, 1) => Some(Pauli::Y),
+                        _ => None,
+                    };
+                    if let Some(pauli) = pauli {
+                        paulis.push((pauli, QubitId::new(qubit)));
+                    }
+                }
+                PauliString::with_phase_and_paulis(QuarterPhase::PlusOne, paulis)
+            })
+            .collect()
+    }
+
+    /// Returns the number of matrix rows.
+    #[must_use]
+    pub fn num_rows(&self) -> usize {
+        self.matrix.num_rows()
+    }
+
+    /// Returns the number of represented qubits.
+    #[must_use]
+    pub fn num_qubits(&self) -> usize {
+        self.matrix.num_cols() / 2
+    }
+
+    /// Returns a copy of the X block.
+    #[must_use]
+    pub fn x_block(&self) -> F2Matrix {
+        let num_qubits = self.num_qubits();
+        let mut block = F2Matrix::zeros(self.num_rows(), num_qubits);
+        for row in 0..self.num_rows() {
+            for qubit in 0..num_qubits {
+                block.set(row, qubit, self.matrix.get(row, qubit));
+            }
+        }
+        block
+    }
+
+    /// Returns a copy of the Z block.
+    #[must_use]
+    pub fn z_block(&self) -> F2Matrix {
+        let num_qubits = self.num_qubits();
+        let mut block = F2Matrix::zeros(self.num_rows(), num_qubits);
+        for row in 0..self.num_rows() {
+            for qubit in 0..num_qubits {
+                block.set(row, qubit, self.matrix.get(row, num_qubits + qubit));
+            }
+        }
+        block
+    }
+
+    /// Returns the rank over GF(2).
+    #[must_use]
+    pub fn rank(&self) -> usize {
+        self.matrix.row_reduce().1.len()
+    }
+
+    /// Returns dense copies of all rows.
+    #[must_use]
+    pub fn rows(&self) -> Vec<Vec<u8>> {
+        self.matrix.rows()
+    }
+
+    /// Returns a dense copy of one row, or `None` if the index is out of range.
+    #[must_use]
+    pub fn row(&self, index: usize) -> Option<Vec<u8>> {
+        (index < self.num_rows()).then(|| self.matrix.row(index))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn pauli_round_trip_ignores_phase_and_preserves_explicit_width() {
+        let negative_y: PauliString = "-Y0".parse().unwrap();
+        let imaginary_xz: PauliString = "+i X1 Z3".parse().unwrap();
+        let sequence = PauliSequence::new(vec![negative_y, imaginary_xz]);
+
+        let matrix = SymplecticMatrix::from_pauli_sequence_ignoring_phase(&sequence, 5).unwrap();
+        assert_eq!(matrix.num_rows(), 2);
+        assert_eq!(matrix.num_qubits(), 5);
+        assert_eq!(matrix.x_block().rows()[0], vec![1, 0, 0, 0, 0]);
+        assert_eq!(matrix.z_block().rows()[0], vec![1, 0, 0, 0, 0]);
+
+        let positive = matrix.to_positive_paulis();
+        assert_eq!(positive[0].get_phase(), QuarterPhase::PlusOne);
+        assert_eq!(positive[1].get_phase(), QuarterPhase::PlusOne);
+        assert_eq!(positive[0].to_dense_str(Some(5)), "+YIIII");
+        assert_eq!(positive[1].to_dense_str(Some(5)), "+IXIZI");
+    }
+
+    #[test]
+    fn empty_dense_input_requires_explicit_zeros_constructor() {
+        assert_eq!(
+            SymplecticMatrix::from_dense(Vec::new()).unwrap_err(),
+            SymplecticMatrixError::EmptyRows
+        );
+        assert_eq!(SymplecticMatrix::zeros(0, 7).num_qubits(), 7);
+    }
+
+    #[test]
+    fn dense_input_rejects_ragged_rows() {
+        assert_eq!(
+            SymplecticMatrix::from_dense(vec![vec![1, 0], vec![1]]).unwrap_err(),
+            SymplecticMatrixError::RaggedRows {
+                row: 1,
+                expected: 2,
+                actual: 1,
+            }
+        );
+    }
+
+    #[test]
+    fn dense_input_rejects_non_binary_entries() {
+        assert_eq!(
+            SymplecticMatrix::from_dense(vec![vec![0, 2]]).unwrap_err(),
+            SymplecticMatrixError::InvalidEntry {
+                row: 0,
+                column: 1,
+                value: 2,
+            }
+        );
+    }
+
+    #[test]
+    fn pauli_sequence_rejects_operator_beyond_explicit_width() {
+        let sequence = PauliSequence::new(vec![PauliString::x(9)]);
+
+        assert_eq!(
+            SymplecticMatrix::from_pauli_sequence_ignoring_phase(&sequence, 3).unwrap_err(),
+            SymplecticMatrixError::OperatorExceedsWidth {
+                row: 0,
+                qubit: 9,
+                num_qubits: 3,
+            }
+        );
+    }
+}
diff --git a/docs/user-guide/stabilizer-code-verification.md b/docs/user-guide/stabilizer-code-verification.md
new file mode 100644
index 000000000..e2040476f
--- /dev/null
+++ b/docs/user-guide/stabilizer-code-verification.md
@@ -0,0 +1,388 @@
+# Stabilizer-Code Verification
+
+This guide covers designing, verifying, and analyzing stabilizer codes with the
+Rust-backed types in `pecos.quantum`. The workflow starts with Pauli checks,
+discovers a compatible logical basis, and searches for low-weight logical
+operators.
+
+## What You'll Learn
+
+- Building a `StabilizerCodeSpec` from Pauli checks
+- Diagnosing anticommuting and dependent generators
+- Discovering logical operators and calculating code distance
+- Searching a range of low-weight logical operators
+- Importing CSS and symplectic check matrices
+- Choosing between the two exact distance methods
+
+```hidden-python
+import re
+
+import numpy as np
+from pecos.quantum import (
+    ParityCheckMatrix,
+    StabilizerCode,
+    StabilizerCodeSpec,
+    SymplecticMatrix,
+    X,
+    Xs,
+    Y,
+    Ys,
+    Z,
+    Zs,
+    pauli_string,
+)
+
+
+def add_checks(builder, checks):
+    for check in checks:
+        builder.check(check)
+    return builder
+
+
+def original_checks():
+    return [
+        Xs([3, 4, 7, 8]),
+        Xs([5, 6, 7, 9]),
+        Zs([2, 4, 5, 7]),
+        Zs([7, 8, 9]),
+        Zs([0, 1]) * Y(2),
+        pauli_string("X0 X2 Z3 Y4"),
+        pauli_string("X1 X2 Z6 Y5"),
+    ]
+
+
+def fixed_checks():
+    checks = original_checks()
+    checks[4] = Zs([0, 1, 2])
+    return checks
+
+
+def final_checks():
+    return [
+        Zs([2, 4, 5, 7]),
+        Xs([3, 4, 7, 8]),
+        Xs([5, 6, 7, 9]),
+        pauli_string("X0 X2 Z3 Y4"),
+        pauli_string("X1 X2 Z6 Y5"),
+        Zs([0, 1, 2]),
+        Xs([0, 1]),
+        Zs([3, 8]),
+        Zs([6, 9]),
+    ]
+```
+
+## Overview
+
+`StabilizerCodeSpec.builder(num_qubits)` collects stabilizer checks and,
+optionally, explicitly chosen logical operators. A check is a `PauliString`.
+The single-qubit `X`, `Y`, and `Z` constructors compose with `&`; `Xs`, `Ys`,
+and `Zs` construct one Pauli type on several qubits at once. Multiplication
+with `*` provides Pauli multiplication, while `pauli_string` parses sparse
+text.
+
+```python
+first = Xs([3, 4, 7, 8])
+assert first == X(3) & X(4) & X(7) & X(8)
+
+mixed = Zs([0, 1]) * Y(2)
+assert mixed == Z(0) & Z(1) & Y(2)
+assert Ys([1, 5]) == Y(1) & Y(5)
+assert pauli_string("X1 X2 Z6 Y5") == X(1) & X(2) & Z(6) & Y(5)
+
+builder = StabilizerCodeSpec.builder(10)
+builder.check(first)
+builder.check(mixed)
+```
+
+Use `build()` when only count and independence validation is needed,
+`build_verified()` to validate a supplied stabilizer and logical basis, or
+`build_with_discovered_logicals()` to verify the checks and discover paired
+logical operators and destabilizers.
+
+## Developing a Ten-Qubit Code
+
+Consider a ten-qubit design with seven proposed checks. The fifth check has a
+`Y` on qubit 2:
+
+```python
+checks = original_checks()
+builder = add_checks(StabilizerCodeSpec.builder(10), checks)
+
+try:
+    builder.build_verified()
+except ValueError as error:
+    message = str(error)
+else:
+    raise AssertionError("the original checks should not verify")
+
+pair = re.search(r"generators (\d+) and (\d+) anticommute", message)
+assert pair is not None
+first_index, second_index = (int(index) for index in pair.groups())
+assert checks[first_index].anticommutes_with(checks[second_index])
+print(message)
+```
+
+```text
+Stabilizer generators 2 and 4 anticommute
+```
+
+The indices identify the offending entries in insertion order. Replacing that
+mixed check with `Zs([0, 1, 2])` produces a valid `[[10, 3]]` code. Building
+with discovered logicals also supplies the destabilizer and paired-logical
+generators displayed by `print(spec)`:
+
+```python
+builder = add_checks(StabilizerCodeSpec.builder(10), fixed_checks())
+spec = builder.build_with_discovered_logicals()
+result = spec.distance()
+
+assert spec.num_logical_qubits == 3
+assert result is not None
+assert result.distance == 2
+assert result.min_weight_operator.weight() == 2
+print(spec)
+print(result)
+```
+
+```text
+[[10, 3]]
+Stabilizer generators:
++IIIXXIIXXI
++IIIIIXXXIX
++IIZIZZIZII
++IIIIIIIZZZ
++ZZZIIIIIII
++XIXZYIIIII
++IXXIIYZIII
+Destabilizer generators:
++IIIZIIIIII
++IZIIIZIIII
++ZIIIXIIIII
++IIIIIIIIXI
++ZIXIXIIIII
++ZIIIIIIIII
++IZIIIIIIII
+Logical operators:
+Z1: +IZIIIZZIII
+X1: +IZIIIIXIII
+Z2: +IZIZIZIZII
+X2: +ZIIIXIIXXI
+Z3: +IZIIIZIIIZ
+X3: +IIIIIIIIXX
+
+DistanceResult(distance=2, min_weight_operator=X_0 X_1)
+```
+
+The parameters line is `[[n, k]]`; the distance result adds the minimum
+logical weight and one operator attaining it. Adding checks that detect the
+weight-two logicals, while removing the original `Zs([7, 8, 9])` check, gives
+the final nine-check design:
+
+```python
+builder = add_checks(StabilizerCodeSpec.builder(10), final_checks())
+spec = builder.build_with_discovered_logicals()
+result = spec.distance()
+
+assert spec.num_logical_qubits == 1
+assert result is not None
+assert result.distance == 3
+assert result.min_weight_operator.weight() == 3
+assert str(spec).splitlines()[0] == "[[10, 1]]"
+print(str(spec).splitlines()[0])
+print(result)
+```
+
+```text
+[[10, 1]]
+DistanceResult(distance=3, min_weight_operator=X_0 X_2 X_7)
+```
+
+This is a `[[10, 1, 3]]` code: it encodes one logical qubit into ten physical
+qubits and has distance three.
+
+## Exploring Low-Weight Logicals
+
+`min_weight_logicals()` returns every logical operator found at the minimum
+weight. Each `LogicalOperatorInfo` records the Pauli operator, its weight, and
+which chosen logical generators it is equivalent to modulo stabilizers.
+`equivalence_string()` formats that last field compactly.
+
+`shortest_logicals(delta)` continues through `delta` weights above the minimum.
+For the five-qubit code there are 30 weight-three logical operators. No
+weight-four logicals exist, and `delta=2` exposes another 18 at weight five:
+
+```python
+spec = StabilizerCodeSpec.from_stabilizer_code(StabilizerCode.five_qubit())
+minimum = spec.min_weight_logicals()
+spectrum = spec.shortest_logicals(delta=2)
+
+assert len(minimum) == 30
+assert [info.operator for info in spectrum[: len(minimum)]] == [
+    info.operator for info in minimum
+]
+assert {info.weight for info in spectrum} == {3, 5}
+assert sum(info.weight == 5 for info in spectrum) == 18
+assert len(spectrum) == 48
+assert all(info.equivalence_string() for info in spectrum)
+
+first = minimum[0]
+print(first.operator, first.weight, first.equivalence_string())
+```
+
+```text
+X_0 Y_1 X_2 3 X0*Z0
+```
+
+Only genuine logical operators are returned; stabilizers and detected
+operators are excluded even when their weights fall inside the requested
+range.
+
+## Matrix Input
+
+For CSS codes, `ParityCheckMatrix` represents a role-neutral binary
+checks-by-qubits matrix. `checks_from_css(x_stabilizers, z_stabilizers)` chooses
+the role: rows in the first matrix become X-type stabilizers, and rows in the
+second become Z-type stabilizers.
+
+### CSS Parity-Check Matrices
+
+The Steane code uses the classical Hamming parity-check matrix for both
+blocks. Nested Python sequences and NumPy integer arrays are accepted:
+
+```python
+hamming_h = [
+    [1, 0, 1, 0, 1, 0, 1],
+    [0, 1, 1, 0, 0, 1, 1],
+    [0, 0, 0, 1, 1, 1, 1],
+]
+
+plain = ParityCheckMatrix(hamming_h)
+int64_matrix = ParityCheckMatrix(np.asarray(hamming_h, dtype=np.int64))
+uint8_matrix = ParityCheckMatrix(np.asarray(hamming_h, dtype=np.uint8))
+assert plain.rows() == int64_matrix.rows() == uint8_matrix.rows() == hamming_h
+
+builder = StabilizerCodeSpec.builder(7)
+builder.checks_from_css(int64_matrix, uint8_matrix)
+steane = builder.build_with_discovered_logicals()
+result = steane.distance(css=True)
+
+assert steane.num_logical_qubits == 1
+assert result is not None
+assert result.distance == StabilizerCode.steane().distance() == 3
+```
+
+The builder checks CSS orthogonality before appending rows. It reports the
+first X-row and Z-row pair with an odd overlap. The code-spec constructor then
+protects the independent-generator invariant and reports both rank and count:
+
+```python
+builder = StabilizerCodeSpec.builder(2)
+try:
+    builder.checks_from_css(
+        ParityCheckMatrix([[1, 0]]),
+        ParityCheckMatrix([[1, 0]]),
+    )
+except ValueError as error:
+    orthogonality_message = str(error)
+else:
+    raise AssertionError("non-orthogonal CSS rows should be rejected")
+assert "X row 0 and Z row 0" in orthogonality_message
+
+dependent = ParityCheckMatrix(
+    [
+        [1, 1, 0],
+        [0, 1, 1],
+        [1, 0, 1],
+    ]
+)
+builder = StabilizerCodeSpec.builder(3)
+builder.checks_from_css(dependent, ParityCheckMatrix.zeros(0, 3))
+try:
+    builder.build()
+except ValueError as error:
+    dependence_message = str(error)
+else:
+    raise AssertionError("dependent stabilizers should be rejected")
+assert "rank 2, count 3" in dependence_message
+```
+
+`ParityCheckMatrix.zeros(0, n)` carries the width that an empty nested list
+cannot express. It is useful for a code with only one stabilizer type:
+
+```python
+x_stabilizers = ParityCheckMatrix([[1, 1]])
+z_stabilizers = ParityCheckMatrix.zeros(0, 2)
+assert z_stabilizers.rows() == []
+assert z_stabilizers.num_qubits() == 2
+
+builder = StabilizerCodeSpec.builder(2)
+builder.checks_from_css(x_stabilizers, z_stabilizers)
+spec = builder.build_with_discovered_logicals()
+assert spec.stabilizers == x_stabilizers.to_x_stabilizers()
+assert spec.num_logical_qubits == 1
+```
+
+### Symplectic Matrices
+
+`SymplecticMatrix` stores each Pauli row as `[X block | Z block]`. A set bit in
+both blocks represents `Y`; phase information is not present, so
+`to_positive_paulis()` always returns positive-phase operators.
+
+These are the four stabilizer rows of the five-qubit code:
+
+```python
+rows = [
+    [1, 0, 0, 1, 0, 0, 1, 1, 0, 0],
+    [0, 1, 0, 0, 1, 0, 0, 1, 1, 0],
+    [1, 0, 1, 0, 0, 0, 0, 0, 1, 1],
+    [0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
+]
+matrix = SymplecticMatrix.from_dense(rows)
+assert matrix.x_block() == [row[:5] for row in rows]
+assert matrix.z_block() == [row[5:] for row in rows]
+
+builder = StabilizerCodeSpec.builder(5)
+builder.checks_from_symplectic(matrix)
+five_qubit = builder.build_with_discovered_logicals()
+result = five_qubit.distance()
+
+assert matrix.to_positive_paulis() == five_qubit.stabilizers
+assert result is not None
+assert result.distance == 3
+```
+
+As with CSS ingestion, width mismatches and anticommuting rows raise
+`ValueError` before a spec is built.
+
+## Choosing a Distance Method
+
+Two exact distance calculations serve different regimes:
+
+| Method | Search strategy | Best use |
+|--------|-----------------|----------|
+| `StabilizerCode.distance()` | Enumerates stabilizer/logical cosets | Tiny codes with small generator counts |
+| `StabilizerCodeSpec.distance()` | Enumerates Paulis by increasing weight | Codes whose distance is small relative to their length |
+
+The coset method is a useful oracle for tiny built-in codes. The spec method
+supports `max_weight` as a search budget and `verbose=True` to print
+`Checking weight N...` progress to standard error. It returns `None` if no
+logical operator is found within the budget:
+
+```python
+code = StabilizerCode.five_qubit()
+spec = StabilizerCodeSpec.from_stabilizer_code(code)
+
+result = spec.distance()
+assert result is not None
+assert result.distance == code.distance() == 3
+assert spec.distance(max_weight=2, verbose=True) is None
+```
+
+The same `max_weight`, `css`, and `verbose` controls are available on
+`min_weight_logicals()` and `shortest_logicals()`.
+
+## Next Steps
+
+- **[Pauli Algebra and QEC in Python](python-pauli-qec.md)** - Work with Pauli strings, sequences, and stabilizer groups
+- **[Stabilizer Codes](stabilizer-codes.md)** - Understand the Rust stabilizer-code model
+- **[QEC Geometry](qec-geometry.md)** - Describe layouts and check supports for code families
diff --git a/mkdocs.yml b/mkdocs.yml
index b58629fb6..c57c08b11 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -67,6 +67,7 @@ nav:
   - Quantum Information Primitives: user-guide/quantum-info.md
   - Stabilizer Codes: user-guide/stabilizer-codes.md
   - Pauli Algebra and QEC in Python: user-guide/python-pauli-qec.md
+  - Stabilizer-Code Verification: user-guide/stabilizer-code-verification.md
   - Fault Tolerance Analysis: user-guide/fault-tolerance.md
   - Fault Catalog Tutorial: user-guide/fault-catalog.md
   - QEC Geometry: user-guide/qec-geometry.md
diff --git a/python/pecos-rslib/pecos_rslib.pyi b/python/pecos-rslib/pecos_rslib.pyi
index f5a792e9f..313c683a1 100644
--- a/python/pecos-rslib/pecos_rslib.pyi
+++ b/python/pecos-rslib/pecos_rslib.pyi
@@ -1060,14 +1060,150 @@ class PauliString:
     def __eq__(self, other: object) -> bool: ...
 
 def X(qubit: int) -> PauliString: ...
+def Xs(qubits: Sequence[int]) -> PauliString: ...
 def Y(qubit: int) -> PauliString: ...
+def Ys(qubits: Sequence[int]) -> PauliString: ...
 def Z(qubit: int) -> PauliString: ...
+def Zs(qubits: Sequence[int]) -> PauliString: ...
 
 class PauliStabilizerGroup:
     """A group of commuting Pauli operators with real phases."""
 
     ...
 
+class StabilizerCode:
+    """A stabilizer group with an explicit physical-qubit count."""
+
+    def __init__(self, group: PauliStabilizerGroup, num_qubits: int | None = None) -> None: ...
+    @staticmethod
+    def repetition(n: int) -> StabilizerCode: ...
+    @staticmethod
+    def steane() -> StabilizerCode: ...
+    @staticmethod
+    def five_qubit() -> StabilizerCode: ...
+    @staticmethod
+    def shor() -> StabilizerCode: ...
+    @staticmethod
+    def four_two_two() -> StabilizerCode: ...
+    @staticmethod
+    def toric(l: int) -> StabilizerCode: ...
+    def num_qubits(self) -> int: ...
+    def num_logical_qubits(self) -> int: ...
+    def code_parameters(self) -> str: ...
+    def distance(self) -> int | None: ...
+    def syndrome(self, error: PauliString) -> list[bool]: ...
+    def logical_operators(self) -> list[PauliString]: ...
+    def group(self) -> PauliStabilizerGroup: ...
+
+class DistanceResult:
+    """A code distance and one minimum-weight logical operator."""
+
+    @property
+    def distance(self) -> int: ...
+    @property
+    def min_weight_operator(self) -> PauliString: ...
+    def __repr__(self) -> str: ...
+
+class LogicalOperatorInfo:
+    """A minimum-weight logical operator and its logical equivalence."""
+
+    @property
+    def operator(self) -> PauliString: ...
+    @property
+    def weight(self) -> int: ...
+    @property
+    def equivalent_logicals(self) -> list[tuple[str, int]]: ...
+    def equivalence_string(self) -> str: ...
+    def __repr__(self) -> str: ...
+
+class ParityCheckMatrix:
+    """A role-neutral binary parity-check matrix."""
+
+    def __init__(self, rows: Sequence[Sequence[int]]) -> None: ...
+    @classmethod
+    def from_dense(cls, rows: Sequence[Sequence[int]]) -> ParityCheckMatrix: ...
+    @classmethod
+    def zeros(cls, num_checks: int, num_qubits: int) -> ParityCheckMatrix: ...
+    def num_checks(self) -> int: ...
+    def num_qubits(self) -> int: ...
+    def rank(self) -> int: ...
+    def rows(self) -> list[list[int]]: ...
+    def to_x_stabilizers(self) -> list[PauliString]: ...
+    def to_z_stabilizers(self) -> list[PauliString]: ...
+    def __repr__(self) -> str: ...
+
+class SymplecticMatrix:
+    """A binary symplectic matrix whose rows represent Pauli operators."""
+
+    def __init__(self, rows: Sequence[Sequence[int]]) -> None: ...
+    @classmethod
+    def from_dense(cls, rows: Sequence[Sequence[int]]) -> SymplecticMatrix: ...
+    @classmethod
+    def zeros(cls, num_rows: int, num_qubits: int) -> SymplecticMatrix: ...
+    def num_rows(self) -> int: ...
+    def num_qubits(self) -> int: ...
+    def rank(self) -> int: ...
+    def rows(self) -> list[list[int]]: ...
+    def x_block(self) -> list[list[int]]: ...
+    def z_block(self) -> list[list[int]]: ...
+    def to_positive_paulis(self) -> list[PauliString]: ...
+    def __repr__(self) -> str: ...
+
+class StabilizerCodeSpecBuilder:
+    """Mutable Python wrapper around the consuming Rust specification builder."""
+
+    def check(self, op: PauliString) -> None: ...
+    def checks_from_css(self, x_stabilizers: ParityCheckMatrix, z_stabilizers: ParityCheckMatrix) -> None: ...
+    def checks_from_symplectic(self, matrix: SymplecticMatrix) -> None: ...
+    def logical_z(self, op: PauliString) -> None: ...
+    def logical_x(self, op: PauliString) -> None: ...
+    def build(self) -> StabilizerCodeSpec: ...
+    def build_verified(self) -> StabilizerCodeSpec: ...
+    def build_with_discovered_logicals(self) -> StabilizerCodeSpec: ...
+
+class StabilizerCodeSpec:
+    """A complete stabilizer-code specification with paired logical operators."""
+
+    def __init__(
+        self,
+        num_qubits: int,
+        stabilizers: list[PauliString],
+        logical_zs: list[PauliString],
+        logical_xs: list[PauliString],
+    ) -> None: ...
+    @staticmethod
+    def builder(num_qubits: int) -> StabilizerCodeSpecBuilder: ...
+    @classmethod
+    def from_stabilizer_code(cls, code: StabilizerCode) -> StabilizerCodeSpec: ...
+    @property
+    def num_qubits(self) -> int: ...
+    @property
+    def num_logical_qubits(self) -> int: ...
+    @property
+    def stabilizers(self) -> list[PauliString]: ...
+    @property
+    def destabilizers(self) -> list[PauliString]: ...
+    @property
+    def logical_zs(self) -> list[PauliString]: ...
+    @property
+    def logical_xs(self) -> list[PauliString]: ...
+    def verify(self) -> None: ...
+    def distance(
+        self, max_weight: int | None = None, css: bool = False, verbose: bool = False
+    ) -> DistanceResult | None: ...
+    def min_weight_logicals(
+        self, max_weight: int | None = None, css: bool = False, verbose: bool = False
+    ) -> list[LogicalOperatorInfo]: ...
+    def shortest_logicals(
+        self,
+        delta: int = 0,
+        max_weight: int | None = None,
+        css: bool = False,
+        verbose: bool = False,
+    ) -> list[LogicalOperatorInfo]: ...
+    def __str__(self) -> str: ...
+    def __repr__(self) -> str: ...
+
 class PauliSequence:
     """Ordered sequence of Pauli operators with symplectic analysis."""
 
diff --git a/python/pecos-rslib/src/code_matrix_bindings.rs b/python/pecos-rslib/src/code_matrix_bindings.rs
new file mode 100644
index 000000000..cdc8578ed
--- /dev/null
+++ b/python/pecos-rslib/src/code_matrix_bindings.rs
@@ -0,0 +1,200 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
+// in compliance with the License. You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software distributed under the License
+// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
+// or implied. See the License for the specific language governing permissions and limitations under
+// the License.
+
+//! Python bindings for binary QEC code matrices.
+
+use pecos_qec::ParityCheckMatrix as RustParityCheckMatrix;
+use pecos_quantum::SymplecticMatrix as RustSymplecticMatrix;
+use pyo3::prelude::*;
+use pyo3::types::PyType;
+
+use crate::pauli_bindings::PauliString;
+
+fn validated_binary_rows(rows: Vec<Vec<i64>>, name: &str) -> PyResult<Vec<Vec<u8>>> {
+    rows.into_iter()
+        .enumerate()
+        .map(|(row_index, row)| {
+            row.into_iter()
+                .map(|value| match value {
+                    0 | 1 => Ok(u8::from(value == 1)),
+                    _ => Err(pyo3::exceptions::PyValueError::new_err(format!(
+                        "{name} row {row_index} contains invalid value {value}; expected 0 or 1"
+                    ))),
+                })
+                .collect()
+        })
+        .collect()
+}
+
+fn python_binary_rows(rows: Vec<Vec<u8>>) -> Vec<Vec<usize>> {
+    rows.into_iter()
+        .map(|row| row.into_iter().map(usize::from).collect())
+        .collect()
+}
+
+/// A role-neutral binary parity-check matrix.
+#[pyclass(name = "ParityCheckMatrix", module = "pecos_rslib", from_py_object)]
+#[derive(Clone, Debug)]
+pub struct PyParityCheckMatrix {
+    pub(crate) inner: RustParityCheckMatrix,
+}
+
+#[pymethods]
+impl PyParityCheckMatrix {
+    #[new]
+    fn new(rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        Self::from_rows(rows)
+    }
+
+    #[classmethod]
+    fn from_dense(_cls: &Bound<'_, PyType>, rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        Self::from_rows(rows)
+    }
+
+    #[classmethod]
+    fn zeros(_cls: &Bound<'_, PyType>, num_checks: usize, num_qubits: usize) -> Self {
+        Self {
+            inner: RustParityCheckMatrix::zeros(num_checks, num_qubits),
+        }
+    }
+
+    fn num_checks(&self) -> usize {
+        self.inner.num_checks()
+    }
+
+    fn num_qubits(&self) -> usize {
+        self.inner.num_qubits()
+    }
+
+    fn rank(&self) -> usize {
+        self.inner.rank()
+    }
+
+    fn rows(&self) -> Vec<Vec<usize>> {
+        python_binary_rows(self.inner.rows())
+    }
+
+    fn to_x_stabilizers(&self) -> Vec<PauliString> {
+        self.inner
+            .to_x_stabilizers()
+            .into_iter()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    fn to_z_stabilizers(&self) -> Vec<PauliString> {
+        self.inner
+            .to_z_stabilizers()
+            .into_iter()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    fn __repr__(&self) -> String {
+        format!(
+            "ParityCheckMatrix(shape=({}, {}))",
+            self.inner.num_checks(),
+            self.inner.num_qubits()
+        )
+    }
+}
+
+impl PyParityCheckMatrix {
+    fn from_rows(rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        let rows = validated_binary_rows(rows, "parity-check matrix")?;
+        let inner = RustParityCheckMatrix::from_dense(rows)
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(Self { inner })
+    }
+}
+
+/// A binary symplectic matrix whose rows represent Pauli operators.
+#[pyclass(name = "SymplecticMatrix", module = "pecos_rslib", from_py_object)]
+#[derive(Clone, Debug)]
+pub struct PySymplecticMatrix {
+    pub(crate) inner: RustSymplecticMatrix,
+}
+
+#[pymethods]
+impl PySymplecticMatrix {
+    #[new]
+    fn new(rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        Self::from_rows(rows)
+    }
+
+    #[classmethod]
+    fn from_dense(_cls: &Bound<'_, PyType>, rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        Self::from_rows(rows)
+    }
+
+    #[classmethod]
+    fn zeros(_cls: &Bound<'_, PyType>, num_rows: usize, num_qubits: usize) -> Self {
+        Self {
+            inner: RustSymplecticMatrix::zeros(num_rows, num_qubits),
+        }
+    }
+
+    fn num_rows(&self) -> usize {
+        self.inner.num_rows()
+    }
+
+    fn num_qubits(&self) -> usize {
+        self.inner.num_qubits()
+    }
+
+    fn rank(&self) -> usize {
+        self.inner.rank()
+    }
+
+    fn rows(&self) -> Vec<Vec<usize>> {
+        python_binary_rows(self.inner.rows())
+    }
+
+    fn x_block(&self) -> Vec<Vec<usize>> {
+        python_binary_rows(self.inner.x_block().rows())
+    }
+
+    fn z_block(&self) -> Vec<Vec<usize>> {
+        python_binary_rows(self.inner.z_block().rows())
+    }
+
+    fn to_positive_paulis(&self) -> Vec<PauliString> {
+        self.inner
+            .to_positive_paulis()
+            .into_iter()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    fn __repr__(&self) -> String {
+        format!(
+            "SymplecticMatrix(shape=({}, {}))",
+            self.inner.num_rows(),
+            self.inner.num_qubits()
+        )
+    }
+}
+
+impl PySymplecticMatrix {
+    fn from_rows(rows: Vec<Vec<i64>>) -> PyResult<Self> {
+        let rows = validated_binary_rows(rows, "symplectic matrix")?;
+        let inner = RustSymplecticMatrix::from_dense(rows)
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(Self { inner })
+    }
+}
+
+pub fn register_code_matrix_types(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_class::<PyParityCheckMatrix>()?;
+    m.add_class::<PySymplecticMatrix>()?;
+    Ok(())
+}
diff --git a/python/pecos-rslib/src/lib.rs b/python/pecos-rslib/src/lib.rs
index 0de139841..f2931443e 100644
--- a/python/pecos-rslib/src/lib.rs
+++ b/python/pecos-rslib/src/lib.rs
@@ -37,6 +37,7 @@ mod bit_int_bindings;
 mod bit_uint_bindings;
 mod byte_message_bindings;
 mod clifford_rep_bindings;
+mod code_matrix_bindings;
 mod coin_toss_bindings;
 mod dag_circuit_bindings;
 mod decoder_bindings;
@@ -71,6 +72,7 @@ mod sparse_stab_engine_bindings;
 mod stab_bindings;
 mod stab_vec_bindings;
 mod stabilizer_code_bindings;
+mod stabilizer_code_spec_bindings;
 mod stabilizer_group_bindings;
 mod state_vec_bindings;
 mod state_vec_engine_bindings;
@@ -302,7 +304,9 @@ fn pecos_rslib(_py: Python<'_>, m: &Bound<'_, PyModule>) -> PyResult<()> {
 
     // Register stabilizer group, Pauli sequence, and Clifford types
     stabilizer_group_bindings::register_stabilizer_group_types(m)?;
+    code_matrix_bindings::register_code_matrix_types(m)?;
     stabilizer_code_bindings::register_stabilizer_code_types(m)?;
+    stabilizer_code_spec_bindings::register_stabilizer_code_spec_types(m)?;
     pauli_sequence_bindings::register_pauli_sequence_types(m)?;
     clifford_rep_bindings::register_clifford_types(m)?;
 
diff --git a/python/pecos-rslib/src/pauli_bindings.rs b/python/pecos-rslib/src/pauli_bindings.rs
index 07a685d4c..9356207fb 100644
--- a/python/pecos-rslib/src/pauli_bindings.rs
+++ b/python/pecos-rslib/src/pauli_bindings.rs
@@ -24,6 +24,7 @@ use std::hash::{Hash, Hasher};
 use crate::prelude::{
     Pauli as RustPauli, PauliOperator, PauliString as RustPauliString, QuarterPhase, QubitId,
 };
+use pecos_core::{Xs as RustXs, Ys as RustYs, Zs as RustZs};
 use pyo3::prelude::*;
 
 /// Single-qubit Pauli operator (I, X, Y, Z)
@@ -312,7 +313,7 @@ impl PauliString {
     }
 
     /// String representation
-    fn __str__(&self) -> String {
+    pub(crate) fn __str__(&self) -> String {
         // Build string representation
         let phase_str = match self.inner.get_phase() {
             QuarterPhase::PlusOne => "",
@@ -595,6 +596,27 @@ pub fn Z(qubit: usize) -> PauliString {
     }
 }
 
+/// Create a multi-qubit X `PauliString`: `Xs([0, 2, 5])`.
+#[pyfunction]
+#[allow(non_snake_case)]
+pub fn Xs(qubits: Vec<usize>) -> PauliString {
+    PauliString::from_rust(RustXs(qubits))
+}
+
+/// Create a multi-qubit Y `PauliString`: `Ys([0, 2, 5])`.
+#[pyfunction]
+#[allow(non_snake_case)]
+pub fn Ys(qubits: Vec<usize>) -> PauliString {
+    PauliString::from_rust(RustYs(qubits))
+}
+
+/// Create a multi-qubit Z `PauliString`: `Zs([0, 2, 5])`.
+#[pyfunction]
+#[allow(non_snake_case)]
+pub fn Zs(qubits: Vec<usize>) -> PauliString {
+    PauliString::from_rust(RustZs(qubits))
+}
+
 /// Register Pauli types with Python module
 pub fn register_pauli_types(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_class::<Pauli>()?;
@@ -602,5 +624,8 @@ pub fn register_pauli_types(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(pyo3::wrap_pyfunction!(X, m)?)?;
     m.add_function(pyo3::wrap_pyfunction!(Y, m)?)?;
     m.add_function(pyo3::wrap_pyfunction!(Z, m)?)?;
+    m.add_function(pyo3::wrap_pyfunction!(Xs, m)?)?;
+    m.add_function(pyo3::wrap_pyfunction!(Ys, m)?)?;
+    m.add_function(pyo3::wrap_pyfunction!(Zs, m)?)?;
     Ok(())
 }
diff --git a/python/pecos-rslib/src/stabilizer_code_bindings.rs b/python/pecos-rslib/src/stabilizer_code_bindings.rs
index d6a0f2781..40545e5e8 100644
--- a/python/pecos-rslib/src/stabilizer_code_bindings.rs
+++ b/python/pecos-rslib/src/stabilizer_code_bindings.rs
@@ -41,7 +41,7 @@ use crate::stabilizer_group_bindings::PyPauliStabilizerGroup;
 #[pyclass(name = "StabilizerCode", module = "pecos_rslib", from_py_object)]
 #[derive(Debug, Clone)]
 pub struct PyStabilizerCode {
-    inner: RustCode,
+    pub(crate) inner: RustCode,
 }
 
 unsafe impl Send for PyStabilizerCode {}
diff --git a/python/pecos-rslib/src/stabilizer_code_spec_bindings.rs b/python/pecos-rslib/src/stabilizer_code_spec_bindings.rs
new file mode 100644
index 000000000..b893b6dd5
--- /dev/null
+++ b/python/pecos-rslib/src/stabilizer_code_spec_bindings.rs
@@ -0,0 +1,411 @@
+// Copyright 2026 The PECOS Developers
+//
+// Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
+// in compliance with the License. You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software distributed under the License
+// is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express
+// or implied. See the License for the specific language governing permissions and limitations under
+// the License.
+
+//! Python bindings for stabilizer-code specifications and distance search.
+
+use pecos_qec::{
+    DistanceResult as RustDistanceResult, DistanceSearchConfig,
+    LogicalOperatorInfo as RustLogicalOperatorInfo, StabilizerCodeSpec as RustCodeSpec,
+    StabilizerCodeSpecBuilder as RustCodeSpecBuilder, calculate_distance,
+    find_min_weight_logicals_with_info, find_shortest_logicals,
+};
+use pyo3::prelude::*;
+use pyo3::types::PyType;
+
+use crate::code_matrix_bindings::{PyParityCheckMatrix, PySymplecticMatrix};
+use crate::pauli_bindings::PauliString;
+use crate::stabilizer_code_bindings::PyStabilizerCode;
+
+/// Result of a stabilizer-code distance search.
+#[pyclass(name = "DistanceResult", module = "pecos_rslib", skip_from_py_object)]
+#[derive(Clone, Debug)]
+pub struct PyDistanceResult {
+    inner: RustDistanceResult,
+}
+
+#[pymethods]
+impl PyDistanceResult {
+    /// The code distance.
+    #[getter]
+    fn distance(&self) -> usize {
+        self.inner.distance
+    }
+
+    /// A logical operator achieving the code distance.
+    #[getter]
+    fn min_weight_operator(&self) -> PauliString {
+        PauliString::from_rust(self.inner.min_weight_operator.clone())
+    }
+
+    fn __repr__(&self) -> String {
+        let operator = PauliString::from_rust(self.inner.min_weight_operator.clone());
+        format!(
+            "DistanceResult(distance={}, min_weight_operator={})",
+            self.inner.distance,
+            operator.__str__()
+        )
+    }
+}
+
+impl From<RustDistanceResult> for PyDistanceResult {
+    fn from(inner: RustDistanceResult) -> Self {
+        Self { inner }
+    }
+}
+
+/// A minimum-weight logical operator and its logical equivalence information.
+#[pyclass(
+    name = "LogicalOperatorInfo",
+    module = "pecos_rslib",
+    skip_from_py_object
+)]
+#[derive(Clone, Debug)]
+pub struct PyLogicalOperatorInfo {
+    inner: RustLogicalOperatorInfo,
+}
+
+#[pymethods]
+impl PyLogicalOperatorInfo {
+    /// The physical Pauli operator.
+    #[getter]
+    fn operator(&self) -> PauliString {
+        PauliString::from_rust(self.inner.operator.clone())
+    }
+
+    /// The physical weight of the operator.
+    #[getter]
+    fn weight(&self) -> usize {
+        self.inner.weight
+    }
+
+    /// Logical operations implemented by the operator.
+    #[getter]
+    fn equivalent_logicals(&self) -> Vec<(String, usize)> {
+        self.inner
+            .equivalent_logicals
+            .iter()
+            .map(|(logical_type, index)| (logical_type.to_string(), *index))
+            .collect()
+    }
+
+    /// Return the logical equivalence as a compact string such as ``X0*Z1``.
+    fn equivalence_string(&self) -> String {
+        self.inner.equivalence_string()
+    }
+
+    fn __repr__(&self) -> String {
+        let operator = PauliString::from_rust(self.inner.operator.clone());
+        format!(
+            "LogicalOperatorInfo(operator={}, weight={}, equivalence={})",
+            operator.__str__(),
+            self.inner.weight,
+            self.inner.equivalence_string()
+        )
+    }
+}
+
+impl From<RustLogicalOperatorInfo> for PyLogicalOperatorInfo {
+    fn from(inner: RustLogicalOperatorInfo) -> Self {
+        Self { inner }
+    }
+}
+
+/// Builder for a stabilizer-code specification.
+#[pyclass(
+    name = "StabilizerCodeSpecBuilder",
+    module = "pecos_rslib",
+    skip_from_py_object
+)]
+pub struct PyStabilizerCodeSpecBuilder {
+    inner: Option<RustCodeSpecBuilder>,
+}
+
+impl PyStabilizerCodeSpecBuilder {
+    fn new(num_qubits: usize) -> Self {
+        Self {
+            inner: Some(RustCodeSpecBuilder::new(num_qubits)),
+        }
+    }
+
+    fn take_inner(&mut self) -> PyResult<RustCodeSpecBuilder> {
+        self.inner.take().ok_or_else(|| {
+            pyo3::exceptions::PyRuntimeError::new_err(
+                "StabilizerCodeSpecBuilder has already been consumed",
+            )
+        })
+    }
+}
+
+#[pymethods]
+impl PyStabilizerCodeSpecBuilder {
+    /// Add a stabilizer generator.
+    fn check(&mut self, op: &PauliString) -> PyResult<()> {
+        let builder = self.take_inner()?;
+        self.inner = Some(builder.check(op.to_rust()));
+        Ok(())
+    }
+
+    /// Add X-type and Z-type stabilizers from CSS parity-check matrices.
+    fn checks_from_css(
+        &mut self,
+        x_stabilizers: &PyParityCheckMatrix,
+        z_stabilizers: &PyParityCheckMatrix,
+    ) -> PyResult<()> {
+        let builder = self.take_inner()?;
+        self.inner = Some(
+            builder
+                .checks_from_css(&x_stabilizers.inner, &z_stabilizers.inner)
+                .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?,
+        );
+        Ok(())
+    }
+
+    /// Add stabilizers from the rows of a symplectic matrix.
+    fn checks_from_symplectic(&mut self, matrix: &PySymplecticMatrix) -> PyResult<()> {
+        let builder = self.take_inner()?;
+        self.inner = Some(
+            builder
+                .checks_from_symplectic(&matrix.inner)
+                .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?,
+        );
+        Ok(())
+    }
+
+    /// Add a logical Z operator.
+    fn logical_z(&mut self, op: &PauliString) -> PyResult<()> {
+        let builder = self.take_inner()?;
+        self.inner = Some(builder.logical_z(op.to_rust()));
+        Ok(())
+    }
+
+    /// Add a logical X operator.
+    fn logical_x(&mut self, op: &PauliString) -> PyResult<()> {
+        let builder = self.take_inner()?;
+        self.inner = Some(builder.logical_x(op.to_rust()));
+        Ok(())
+    }
+
+    /// Build with count validation only.
+    fn build(&mut self) -> PyResult<PyStabilizerCodeSpec> {
+        let inner = self
+            .take_inner()?
+            .build()
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(PyStabilizerCodeSpec { inner })
+    }
+
+    /// Build and fully verify all commutation relations.
+    fn build_verified(&mut self) -> PyResult<PyStabilizerCodeSpec> {
+        let inner = self
+            .take_inner()?
+            .build_verified()
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(PyStabilizerCodeSpec { inner })
+    }
+
+    /// Build and automatically discover paired logical operators.
+    fn build_with_discovered_logicals(&mut self) -> PyResult<PyStabilizerCodeSpec> {
+        let inner = self
+            .take_inner()?
+            .build_with_discovered_logicals()
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(PyStabilizerCodeSpec { inner })
+    }
+}
+
+/// A complete stabilizer-code specification with paired logical operators.
+#[pyclass(name = "StabilizerCodeSpec", module = "pecos_rslib", from_py_object)]
+#[derive(Clone, Debug)]
+pub struct PyStabilizerCodeSpec {
+    inner: RustCodeSpec,
+}
+
+#[pymethods]
+impl PyStabilizerCodeSpec {
+    /// Create a builder for a code with the specified number of qubits.
+    #[staticmethod]
+    fn builder(num_qubits: usize) -> PyStabilizerCodeSpecBuilder {
+        PyStabilizerCodeSpecBuilder::new(num_qubits)
+    }
+
+    /// Create a stabilizer-code specification.
+    #[new]
+    fn new(
+        num_qubits: usize,
+        stabilizers: Vec<PauliString>,
+        logical_zs: Vec<PauliString>,
+        logical_xs: Vec<PauliString>,
+    ) -> PyResult<Self> {
+        let stabilizers = stabilizers
+            .into_iter()
+            .map(|pauli| pauli.to_rust())
+            .collect();
+        let logical_zs = logical_zs
+            .into_iter()
+            .map(|pauli| pauli.to_rust())
+            .collect();
+        let logical_xs = logical_xs
+            .into_iter()
+            .map(|pauli| pauli.to_rust())
+            .collect();
+        let inner = RustCodeSpec::new(num_qubits, stabilizers, logical_zs, logical_xs)
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(Self { inner })
+    }
+
+    /// Create a full specification from a ``StabilizerCode``.
+    #[classmethod]
+    fn from_stabilizer_code(_cls: &Bound<'_, PyType>, code: &PyStabilizerCode) -> PyResult<Self> {
+        let inner = RustCodeSpec::from_stabilizer_code(&code.inner)
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))?;
+        Ok(Self { inner })
+    }
+
+    /// Number of physical qubits.
+    #[getter]
+    fn num_qubits(&self) -> usize {
+        self.inner.num_qubits()
+    }
+
+    /// Number of encoded logical qubits.
+    #[getter]
+    fn num_logical_qubits(&self) -> usize {
+        self.inner.num_logical_qubits()
+    }
+
+    /// Stabilizer generators.
+    #[getter]
+    fn stabilizers(&self) -> Vec<PauliString> {
+        self.inner
+            .stabilizers()
+            .iter()
+            .cloned()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    /// Destabilizer generators.
+    #[getter]
+    fn destabilizers(&self) -> Vec<PauliString> {
+        self.inner
+            .destabilizers()
+            .iter()
+            .cloned()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    /// Logical Z operators.
+    #[getter]
+    fn logical_zs(&self) -> Vec<PauliString> {
+        self.inner
+            .logical_zs()
+            .iter()
+            .cloned()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    /// Logical X operators.
+    #[getter]
+    fn logical_xs(&self) -> Vec<PauliString> {
+        self.inner
+            .logical_xs()
+            .iter()
+            .cloned()
+            .map(PauliString::from_rust)
+            .collect()
+    }
+
+    /// Verify all stabilizer and logical commutation relations.
+    fn verify(&self) -> PyResult<()> {
+        self.inner
+            .verify()
+            .map_err(|error| pyo3::exceptions::PyValueError::new_err(error.to_string()))
+    }
+
+    /// Find the code distance and one minimum-weight logical operator.
+    #[pyo3(signature = (max_weight=None, css=false, verbose=false))]
+    fn distance(
+        &self,
+        max_weight: Option<usize>,
+        css: bool,
+        verbose: bool,
+    ) -> Option<PyDistanceResult> {
+        let config = DistanceSearchConfig {
+            max_weight,
+            css_only: css,
+            verbose,
+        };
+        calculate_distance(&self.inner, &config).map(PyDistanceResult::from)
+    }
+
+    /// Find all logical operators at the minimum weight searched.
+    #[pyo3(signature = (max_weight=None, css=false, verbose=false))]
+    fn min_weight_logicals(
+        &self,
+        max_weight: Option<usize>,
+        css: bool,
+        verbose: bool,
+    ) -> Vec<PyLogicalOperatorInfo> {
+        let config = DistanceSearchConfig {
+            max_weight,
+            css_only: css,
+            verbose,
+        };
+        find_min_weight_logicals_with_info(&self.inner, &config)
+            .into_iter()
+            .map(PyLogicalOperatorInfo::from)
+            .collect()
+    }
+
+    /// Find logical operators through ``delta`` weights above the minimum.
+    #[pyo3(signature = (delta=0, max_weight=None, css=false, verbose=false))]
+    fn shortest_logicals(
+        &self,
+        delta: usize,
+        max_weight: Option<usize>,
+        css: bool,
+        verbose: bool,
+    ) -> Vec<PyLogicalOperatorInfo> {
+        let config = DistanceSearchConfig {
+            max_weight,
+            css_only: css,
+            verbose,
+        };
+        find_shortest_logicals(&self.inner, &config, delta)
+            .into_iter()
+            .map(PyLogicalOperatorInfo::from)
+            .collect()
+    }
+
+    fn __str__(&self) -> String {
+        self.inner.to_string()
+    }
+
+    fn __repr__(&self) -> String {
+        format!(
+            "StabilizerCodeSpec([[{}, {}]])",
+            self.inner.num_qubits(),
+            self.inner.num_logical_qubits()
+        )
+    }
+}
+
+/// Register stabilizer-code specification and distance-result types.
+pub fn register_stabilizer_code_spec_types(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_class::<PyStabilizerCodeSpec>()?;
+    m.add_class::<PyStabilizerCodeSpecBuilder>()?;
+    m.add_class::<PyDistanceResult>()?;
+    m.add_class::<PyLogicalOperatorInfo>()?;
+    Ok(())
+}
diff --git a/python/quantum-pecos/src/pecos/quantum/__init__.py b/python/quantum-pecos/src/pecos/quantum/__init__.py
index 6cd99b993..ba00c26c1 100644
--- a/python/quantum-pecos/src/pecos/quantum/__init__.py
+++ b/python/quantum-pecos/src/pecos/quantum/__init__.py
@@ -108,6 +108,7 @@
         SZ,
         SZZ,
         CliffordRep,
+        DistanceResult,
         F,
         F2dg,
         F3dg,
@@ -118,21 +119,29 @@
         GateDefBuilder,
         GateRegistry,
         H,
+        LogicalOperatorInfo,
+        ParityCheckMatrix,
         Pauli,
         PauliSequence,
         PauliStabilizerGroup,
         PauliString,
         StabilizerCode,
+        StabilizerCodeSpec,
+        StabilizerCodeSpecBuilder,
         SXdg,
         SXXdg,
         SYdg,
+        SymplecticMatrix,
         SYYdg,
         SZdg,
         SZZdg,
         TableauWrapper,
         X,
+        Xs,
         Y,
+        Ys,
         Z,
+        Zs,
         adjust_tableau_string,
         sparse_stab,
     )
@@ -290,6 +299,7 @@ def pauli_string(
     "CliffordRep",
     "DagCircuit",
     "DagCircuitWouldCycleError",
+    "DistanceResult",
     "F",
     "F2dg",
     "F3dg",
@@ -305,6 +315,8 @@ def pauli_string(
     "HostedGateRecord",
     "HostedOperationBinding",
     "HugrConversionError",
+    "LogicalOperatorInfo",
+    "ParityCheckMatrix",
     "Pauli",
     "PauliSequence",
     "PauliStabilizerGroup",
@@ -318,6 +330,9 @@ def pauli_string(
     "SZZdg",
     "SZdg",
     "StabilizerCode",
+    "StabilizerCodeSpec",
+    "StabilizerCodeSpecBuilder",
+    "SymplecticMatrix",
     "TableauWrapper",
     "Tick",
     "TickCircuit",
@@ -325,8 +340,11 @@ def pauli_string(
     "TickMeasureHandle",
     "TickPrepHandle",
     "X",
+    "Xs",
     "Y",
+    "Ys",
     "Z",
+    "Zs",
     "adjust_tableau_string",
     "commute",
     "gate_groups",
diff --git a/python/quantum-pecos/src/pecos/tools/fault_tolerance_checks.py b/python/quantum-pecos/src/pecos/tools/fault_tolerance_checks.py
deleted file mode 100644
index fc1fa13fa..000000000
--- a/python/quantum-pecos/src/pecos/tools/fault_tolerance_checks.py
+++ /dev/null
@@ -1,457 +0,0 @@
-"""Fault tolerance verification for quantum error correction."""
-
-# Copyright 2018 The PECOS Developers
-# Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract
-# DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software.
-#
-# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with
-# the License.You may obtain a copy of the License at
-#
-#     https://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an
-# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the
-# specific language governing permissions and limitations under the License.
-from __future__ import annotations
-
-import itertools as it
-from itertools import combinations, product
-from typing import TYPE_CHECKING, TypeVar
-
-import pecos as pc
-from pecos.analysis.stabilizer_funcs import circ2set, find_stab, op_commutes, remove_stab
-from pecos.circuits import LogicalCircuit, QuantumCircuit
-from pecos.decoders import MWPM2D
-from pecos.engines.circuit_runners import Standard
-from pecos.noise.parent_class_error_gen import ErrorCircuits
-from pecos.simulators import SparseStabPy
-
-if TYPE_CHECKING:
-    from collections.abc import Iterable
-
-    from pecos.protocols import Decoder, QECCProtocol, SimulatorProtocol
-
-T = TypeVar("T")
-
-
-def powerset(
-    iterable: Iterable[T],
-    bound: int | None = None,
-) -> it.chain[tuple[T, ...]]:
-    """Returns the power set of an iterable."""
-    powerlist = list(iterable)
-    if bound is None:
-        bound = len(powerlist)
-    return it.chain.from_iterable(it.combinations(powerlist, t) for t in range(bound + 1))
-
-
-def t_errors_check(
-    qecc: QECCProtocol,
-    logical_gate: QuantumCircuit | LogicalCircuit | None = None,
-    syn_extract: QuantumCircuit | LogicalCircuit | None = None,
-    decoder: Decoder | None = None,
-    t_weight: int | None = None,
-    error_set: Iterable[tuple[set[int], set[int]]] | None = None,
-    *,
-    verbose: bool = True,
-    data_errors: bool = True,
-    ancilla_errors: bool = False,
-) -> tuple[bool, int]:
-    """Check exRec conditions for fault-free error correction or logical gate.
-
-    This checks that the exRec conditions for a fault-free error correction (EC) or logical gate (Ga) as described in
-    arXiv:quant-ph/0504218.
-
-    For fault-free EC, weight <= t errors in produce no errors out.
-
-    For fault-free Ga, weight <= t errors in produce weight <= errors out.
-
-
-    Fault-free EC:
-                     ------------------
-    error wt <= t -> |EC (fault free) | -> no errors => No syndrome in subsequent fault-free EC (+ no logical faults)
-                     ------------------
-
-
-    Fault-free Ga:
-                     ------------------
-    error wt <= t -> |Ga (fault free) | ->error wt <= t   => A following fault-free EC + Recovery will result in a state
-                     ------------------
-    with no logical fault.
-
-
-    Args:
-    ----
-        qecc: The quantum error correcting code instance.
-        logical_gate(QuantumCircuit): The logical gate circuit to test (None for error correction only).
-        syn_extract(QuantumCircuit): The syndrome extraction circuit to use.
-        decoder: The decoder instance for error correction.
-        t_weight: The maximum weight of errors to check (typically pc.floor((distance-1)/2)).
-        error_set: Custom set of errors to check (if None, all Pauli errors are checked).
-        verbose: If True, prints detailed information about failures.
-        data_errors: If True, includes errors on data qubits.
-        ancilla_errors: If True, includes errors on ancilla qubits.
-
-    Returns:
-    -------
-        tuple (bool, int): The bool is whether the check is passed. The int is the weight of error last checked. If the
-        bool is True then int == t_weight. If bool == False, int == weight of error that caused a logical error.
-
-    """
-    qudit_set = set()
-
-    if data_errors:
-        qudit_set.update(qecc.data_qudit_set)
-
-    if ancilla_errors:
-        qudit_set.update(qecc.ancilla_qudit_set)
-
-    if t_weight is None:
-        t_weight = pc.floor((qecc.distance - 1) / 2)
-
-    if error_set is None:
-        error_set = {"X", "Y", "Z"}
-
-    circ_sim = Standard()
-
-    # init |0> circuit
-    initzero = LogicalCircuit(suppress_warning=True)
-    initzero.append(qecc.gate("ideal init |0>"))
-
-    # init |+> circuit
-    initplus = LogicalCircuit(suppress_warning=True)
-    initplus.append(qecc.gate("ideal init |+>"))
-
-    if syn_extract is not None and logical_gate is not None:
-        msg = "Both syn_extract and logical_gate cannot be set (not None)."
-        raise Exception(msg)
-
-    if syn_extract is None:
-        # Syndrome extraction
-        syn_extract = LogicalCircuit(suppress_warning=True)
-        syn_extract.append(qecc.gate("I", num_syn_extract=1, forced_outcome=1))
-
-    logic = syn_extract if logical_gate is None else logical_gate
-
-    logical_ops_zero = qecc.instruction("instr_init_zero").logical_stabs[0]
-    logical_ops_plus = qecc.instruction("instr_init_plus").logical_stabs[0]
-
-    if decoder is None:
-        decoder = MWPM2D(qecc)
-
-    for qubit_comb in powerset(qudit_set):
-        if len(qubit_comb) > t_weight:
-            break
-
-        error_combinations = product(error_set, repeat=len(qubit_comb))
-
-        for error_comb in error_combinations:
-            error_circ = QuantumCircuit(1)
-            errors = ErrorCircuits()
-
-            errors.simple_add(0, 0, 0, before_errors=error_circ)
-
-            for e, q in zip(error_comb, qubit_comb, strict=False):
-                error_circ.update(e, {q})
-
-            state_zero = SparseStabPy(qecc.num_qudits)
-            state_plus = SparseStabPy(qecc.num_qudits)
-
-            circ_sim.run(state_zero, initzero)
-            circ_sim.run(state_plus, initplus)
-
-            output, _ = circ_sim.run(state_zero, logic, error_circuits=errors)
-            circ_sim.run(state_plus, logic, error_circuits=errors)
-
-            syn = output.simplified(last=True)
-
-            if syn:
-                # Recovery operation
-                recovery = decoder.decode(syn)
-                circ_sim.run(state_zero, recovery)
-                circ_sim.run(state_plus, recovery)
-
-            sign_zero = state_zero.logical_sign(*logical_ops_zero)
-            sign_plus = state_plus.logical_sign(*logical_ops_plus)
-
-            if sign_zero or sign_plus:
-                if verbose:
-                    print(errors)
-                return False, len(error_comb)
-
-            if logical_gate is None:  # The following is only required for EC.
-                # Any remaining syndromes?
-                output, _ = circ_sim.run(state_zero, syn_extract)
-                syn = output.simplified(last=True)
-
-                if syn:
-                    if verbose:
-                        print(f"syndromes = {syn}")
-                        print(errors)
-                    return False, len(error_comb)
-
-    return True, int(t_weight)
-
-
-def fault_check(
-    qecc: QECCProtocol,
-    logical_gate: QuantumCircuit | LogicalCircuit | None = None,
-    decoder: Decoder | None = None,
-    t_weight: int | None = None,
-    error_set: Iterable[tuple[set[int], set[int]]] | None = None,
-    *,
-    verbose: bool = True,
-    data_errors: bool = True,
-    ancilla_errors: bool = False,
-) -> tuple[bool, int]:
-    """Check exRec conditions for faulty error correction or logical gate.
-
-    This checks that the exRec conditions for a faulty error correction (EC) or logical gate (Ga) as described in
-    arXiv:quant-ph/0504218.
-
-    For fault-free EC, weight <= t errors in produce no errors out.
-
-    For fault-free Ga, weight <= t errors in produce weight <= errors out.
-
-
-    Fault-free EC:
-                     ------------------
-    error wt <= t -> |EC (fault free) | -> no errors => No syndrome in subsequent fault-free EC (+ no logical faults)
-                     ------------------
-
-
-    Fault-free Ga:
-                     ------------------
-    error wt <= t -> |Ga (fault free) | ->error wt <= t   => A following fault-free EC + Recovery will result in a state
-                     ------------------
-    with no logical fault.
-
-
-    Args:
-    ----
-        qecc: The quantum error correcting code instance.
-        logical_gate(QuantumCircuit): The logical gate circuit to test (None for error correction only).
-        decoder: The decoder instance for error correction.
-        t_weight: The maximum weight of errors to check (typically pc.floor((distance-1)/2)).
-        error_set: Custom set of errors to check (if None, all Pauli errors are checked).
-        verbose: If True, prints detailed information about failures.
-        data_errors: If True, includes errors on data qubits.
-        ancilla_errors: If True, includes errors on ancilla qubits.
-
-    Returns:
-    -------
-        tuple (bool, int): The bool is whether the check is passed. The int is the weight of error last checked. If the
-        bool is True then int == t_weight. If bool == False, int == weight of error that caused a logical error.
-
-    """
-    qudit_set = set()
-
-    if data_errors:
-        qudit_set.update(qecc.data_qudit_set)
-
-    if ancilla_errors:
-        qudit_set.update(qecc.ancilla_qudit_set)
-
-    if t_weight is None:
-        t_weight = pc.floor((qecc.distance - 1) / 2)
-
-    if error_set is None:
-        error_set = {"X", "Y", "Z"}
-
-    circ_sim = Standard()
-
-    # init |0> circuit
-    initzero = LogicalCircuit(suppress_warning=True)
-    initzero.append(qecc.gate("ideal init |0>"))
-
-    # init |+> circuit
-    initplus = LogicalCircuit(suppress_warning=True)
-    initplus.append(qecc.gate("ideal init |+>"))
-
-    if logical_gate is None:
-        # Syndrome extraction
-        syn_extract = LogicalCircuit(suppress_warning=True)
-        syn_extract.append(qecc.gate("I", num_syn_extract=1, forced_outcome=1))
-        logic = syn_extract
-    else:
-        logic = logical_gate
-
-    logical_ops_zero = qecc.instruction("instr_init_zero").logical_stabs[0]
-    logical_ops_plus = qecc.instruction("instr_init_plus").logical_stabs[0]
-
-    if decoder is None:
-        decoder = MWPM2D(qecc)
-
-    for qubit_comb in powerset(qudit_set):
-        if len(qubit_comb) > t_weight:
-            break
-
-        error_combinations = product(error_set, repeat=len(qubit_comb))
-
-        for error_comb in error_combinations:
-            error_circ = QuantumCircuit(1)
-            errors = ErrorCircuits()
-
-            errors.simple_add(0, 0, 0, before_errors=error_circ)
-
-            for e, q in zip(error_comb, qubit_comb, strict=False):
-                error_circ.update(e, {q})
-
-            state_zero = SparseStabPy(qecc.num_qudits)
-            state_plus = SparseStabPy(qecc.num_qudits)
-
-            circ_sim.run(state_zero, initzero)
-            circ_sim.run(state_plus, initplus)
-
-            output, _ = circ_sim.run(state_zero, logic, error_circuits=errors)
-            circ_sim.run(state_plus, logic, error_circuits=errors)
-
-            syn = output.simplified(last=True)
-
-            if syn:
-                # Recovery operation
-                recovery = decoder.decode(syn)
-                circ_sim.run(state_zero, recovery)
-                circ_sim.run(state_plus, recovery)
-
-            sign_zero = state_zero.logical_sign(*logical_ops_zero)
-            sign_plus = state_plus.logical_sign(*logical_ops_plus)
-
-            if sign_zero or sign_plus:
-                if verbose:
-                    print(errors)
-                return False, len(error_comb)
-
-    return True, int(t_weight)
-
-
-def distance_check(
-    qecc: QECCProtocol,
-    mode: str | None = None,
-    dist_mode: str | None = None,
-) -> int:
-    """Determines the distance of the code by looking for the smallest logical errors.
-
-    Args:
-    ----
-        qecc: The quantum error correcting code instance.
-        mode: The mode for distance checking ('X', 'x', 'Z', 'z', or None for automatic).
-        dist_mode: The specific distance checking mode to use (if None, uses default based on mode).
-
-    Returns:
-    -------
-        Tuple (bool, int). The bool is whether the check is passed. The int is the weight of error last checked. If the
-        bool is True then int == t_weight. If bool == False, int == weight of error that caused a logical error.
-
-    """
-    qudit_set = qecc.data_qudit_set
-
-    circ_sim = Standard()
-    state = SparseStabPy(qecc.num_qudits)
-
-    ideal_initlogic = LogicalCircuit(suppress_warning=True)
-    ideal_initlogic.append(qecc.gate("ideal init |0>"))
-
-    circ_sim.run(state, ideal_initlogic)
-
-    logical_op, delogical_op = qecc.instruction("instr_init_zero").logical_stabs[0]
-
-    destab_xs, destab_zs = circ2set(delogical_op.items(params=False))
-    stab_xs, stab_zs = circ2set(logical_op.items(params=False))
-
-    remove_stab(state, stab_xs, stab_zs, destab_xs, destab_zs)
-
-    if dist_mode is None:
-        if mode in {"X", "x"}:
-            print("x")
-            return dist_mode_x(state, qudit_set)
-        if mode in {"Z", "z"}:
-            print("z")
-            return dist_mode_z(state, qudit_set)
-        if mode == "power":
-            return dist_mode_powerset(state, qudit_set)
-        return dist_mode_smallest(state, qudit_set)
-
-    return dist_mode(state, qudit_set)
-
-
-def dist_mode_powerset(state: SimulatorProtocol, qudit_set: set[int]) -> str | bool:
-    """Check for logical errors using powerset of all possible X and Z errors.
-
-    Args:
-    state: Stabilizer state to check.
-    qudit_set: Set of qudit/qubit indices to check for errors.
-    """
-    for x_errors in powerset(qudit_set):
-        for z_errors in powerset(qudit_set):
-            if op_commutes(x_errors, z_errors, state.stabs) and not find_stab(
-                state,
-                x_errors,
-                z_errors,
-            ):
-                return f"Logical error found: Xs - {x_errors} Zs - {z_errors}"
-
-    return False
-
-
-def dist_mode_smallest(state: SimulatorProtocol, qudit_set: set[int]) -> str | bool:
-    """Find smallest logical error by checking errors in increasing size.
-
-    Args:
-    ----
-        state: Stabilizer state to check.
-        qudit_set: Set of qudit/qubit indices to check for errors.
-    """
-    for lenq in range(len(qudit_set) + 1):
-        for qs in combinations(qudit_set, lenq):
-            if op_commutes(qs, qs, state.stabs) and not find_stab(state, qs, qs):
-                return f"Logical error found: Xs - {qs} Zs - {qs}"
-
-            for qs2 in powerset(qudit_set, len(qs) - 1):
-                if op_commutes(qs2, qs, state.stabs) and not find_stab(state, qs2, qs):
-                    return f"Logical error found: Xs - {qs2} Zs - {qs}"
-
-                if op_commutes(qs, qs2, state.stabs) and not find_stab(state, qs, qs2):
-                    return f"Logical error found: Xs - {qs} Zs - {qs2}"
-
-    return False
-
-
-def dist_mode_x(state: SimulatorProtocol, qudit_set: set[int]) -> str | bool:
-    """Check for X-type logical errors only.
-
-    Args:
-    ----
-        state: Stabilizer state to check.
-        qudit_set: Set of qudit/qubit indices to check for errors.
-    """
-    z_errors = ()
-    for x_errors in powerset(qudit_set):
-        if op_commutes(x_errors, z_errors, state.stabs) and not find_stab(
-            state,
-            x_errors,
-            z_errors,
-        ):
-            return f"Logical error found: Xs - {x_errors} Zs - {z_errors}"
-
-    return False
-
-
-def dist_mode_z(state: SimulatorProtocol, qudit_set: set[int]) -> str | bool:
-    """Check for Z-type logical errors only.
-
-    Args:
-    ----
-        state: Stabilizer state to check.
-        qudit_set: Set of qudit/qubit indices to check for errors.
-    """
-    x_errors = ()
-    for z_errors in powerset(qudit_set):
-        if op_commutes(x_errors, z_errors, state.stabs) and not find_stab(
-            state,
-            x_errors,
-            z_errors,
-        ):
-            return f"Logical error found: Xs - {x_errors} Zs - {z_errors}"
-
-    return False
diff --git a/python/quantum-pecos/src/pecos/tools/stabilizer_verification.py b/python/quantum-pecos/src/pecos/tools/stabilizer_verification.py
deleted file mode 100644
index 0feb0974f..000000000
--- a/python/quantum-pecos/src/pecos/tools/stabilizer_verification.py
+++ /dev/null
@@ -1,1117 +0,0 @@
-"""Stabilizer verification tools for quantum error correction.
-
-This module provides utilities for verifying stabilizer codes and analyzing
-their properties, including stabilizer group verification, code distance
-calculation, and logical operator validation.
-"""
-
-# Copyright 2018 The PECOS Developers
-# Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS). Under the terms of Contract
-# DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software.
-#
-# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with
-# the License.You may obtain a copy of the License at
-#
-#     https://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an
-# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the
-# specific language governing permissions and limitations under the License.
-
-from __future__ import annotations
-
-from itertools import combinations, product
-from typing import TYPE_CHECKING
-
-import pecos as pc
-from pecos.circuits import QuantumCircuit
-
-if TYPE_CHECKING:
-    from collections.abc import Generator, Sequence
-
-    from pecos.protocols import SimulatorProtocol
-    from pecos.typing import LogicalOpInfo, StabilizerVerificationResult
-
-# TODO: NEED TO ADD SIGN TRACKING TO DESTABILIZERS TO GET THE RIGHT SIGN FOR LOGICAL Xs
-
-
-class VerifyStabilizers:
-    """Used to define a stabilizer QECC."""
-
-    def __init__(self) -> None:
-        """Initialize the VerifyStabilizers instance.
-
-        Sets up the circuit simulator and initializes empty data structures
-        for stabilizer checks, logical operators, and qubit tracking.
-        """
-        self.circ_sim = pc.simulators.SparseStabPy
-
-        self.checks = []
-        self.logical_zs = []
-        self.logical_xs = []
-        self.logical_zs_defined = []  # User chosen logical Zs TODO: ...
-        self.logical_xs_defined = []  # User chosen logical Xs TODO: ...
-        self.logical_zs_reference = {}
-        self.logical_xs_reference = {}
-
-        self.data_qubits = set()
-        self.ancilla_qubits = set()
-        self.circuit = None
-
-        self.check_row_x = None
-        self.check_row_z = None
-        self.check_col_x = None
-        self.check_col_z = None
-
-        # stabilizer ids:
-        self.data_gens = None
-        self.ancilla_gens = None
-        self.logical_gens = None
-
-        # distance
-        self.dist = None
-
-        self.state = None
-
-    def check(self, paulis: str | Sequence[str], qubits: Sequence[int]) -> None:
-        """Check if the given Pauli operators stabilize the state.
-
-        Args:
-            paulis: Sequence of Pauli operators to check (e.g., ['X', 'Z', 'Y']).
-            qubits: Sequence of qubit indices corresponding to the Pauli operators.
-
-        Returns: None
-
-        """
-        if not qubits:
-            msg = "No qubit ids given."
-            raise Exception(msg)
-
-        check_string = "check(" + str(paulis) + ", " + str(qubits) + ")"
-
-        if isinstance(paulis, str):
-            if paulis not in {"X", "x", "Y", "y", "Z", "z"}:
-                msg = 'Paulis should be "X", "Y" or "Z"!'
-                raise Exception(msg)
-
-            paulis_new = [paulis for _ in qubits]
-            paulis = paulis_new
-
-        if not isinstance(paulis, str) and len(paulis) != len(qubits):
-            msg = "Number of Paulis and qubits do not match!!!"
-            raise Exception(msg)
-
-        self.checks.append((paulis, qubits, check_string))
-        self.data_qubits.update(qubits)
-
-    def logicalz(self, paulis: str | Sequence[str], qubits: Sequence[int]) -> None:
-        """Used to define logical Z.
-
-        Args:
-        ----
-            paulis: Either a single Pauli string ('X', 'Y', or 'Z') or a list of Pauli operators.
-            qubits: List of qubit indices where the logical Z operator acts.
-        """
-        if not qubits:
-            msg = "No qubit ids given."
-            raise Exception(msg)
-
-        logical_string = "check(" + str(paulis) + ", " + str(qubits) + ")"
-
-        if isinstance(paulis, str):
-            if paulis not in {"X", "x", "Y", "y", "Z", "z"}:
-                msg = 'Paulis should be "X", "Y" or "Z"!'
-                raise Exception(msg)
-
-            paulis_new = [paulis for _ in qubits]
-            paulis = paulis_new
-
-        if not isinstance(paulis, str) and len(paulis) != len(qubits):
-            msg = "Number of Paulis and qubits do not match!!!"
-            raise Exception(msg)
-
-        self.logical_zs.append((paulis, qubits, logical_string))
-
-    def logicalx(self, paulis: str | Sequence[str], qubits: Sequence[int]) -> None:
-        """Used to define logical X.
-
-        Args:
-        ----
-            paulis: Either a single Pauli string ('X', 'Y', or 'Z') or a list of Pauli operators.
-            qubits: List of qubit indices where the logical X operator acts.
-        """
-        if not qubits:
-            msg = "No qubit ids given."
-            raise Exception(msg)
-
-        logical_string = "check(" + str(paulis) + ", " + str(qubits) + ")"
-
-        if isinstance(paulis, str):
-            if paulis not in {"X", "x", "Y", "y", "Z", "z"}:
-                msg = 'Paulis should be "X", "Y" or "Z"!'
-                raise Exception(msg)
-
-            paulis_new = [paulis for _ in qubits]
-            paulis = paulis_new
-
-        if not isinstance(paulis, str) and len(paulis) != len(qubits):
-            msg = "Number of Paulis and qubits do not match!!!"
-            raise Exception(msg)
-
-        self.logical_xs.append((paulis, qubits, logical_string))
-
-    def num_logical_qubits(self) -> int:
-        """Calculate the number of logical qubits in the stabilizer code.
-
-        Returns:
-            Number of logical qubits (data qubits minus stabilizer checks).
-        """
-        return len(self.data_qubits) - len(self.checks)
-
-    def generators(
-        self,
-        *,
-        print_y: bool = True,
-        verbose: bool = True,
-    ) -> tuple[
-        list[dict[str, set[int]]],
-        list[dict[str, set[int]]],
-        list[str],
-        list[str],
-    ]:
-        """Evaluates the stabilizer generators that have been supplied via the `check` method.
-
-        Args:
-        ----
-            print_y: If True, includes Y operators in the output (otherwise converts to X and Z).
-            verbose: If True, prints detailed information about the generators.
-        """
-        if self.circuit is None:
-            msg = "Must compile circuits first!"
-            raise Exception(msg)
-
-        state = self.state
-        z, x, stab_strings, destab_strings = self.get_info(
-            state,
-            print_y=print_y,
-            verbose=verbose,
-        )
-
-        return z, x, stab_strings, destab_strings
-
-    def _check_all_labels(self) -> None:
-        """This checks to see that all the consecutive qubit ids have been used and none are missing."""
-        qubit_labels = set()
-
-        checks = self.checks
-        for check in checks:
-            _, qs, _ = check
-            qubit_labels.update(qs)
-
-        largest_labels = max(qubit_labels)
-        labels_should_have = set(range(largest_labels + 1))
-
-        dont_have = labels_should_have - qubit_labels
-
-        if dont_have:
-            msg = f"Qubit ids missing: {dont_have}"
-            raise Exception(msg)
-
-    def _check2rowcol(self) -> None:
-        """Creates row and column matrices."""
-        checks = self.checks
-
-        num_checks = len(checks)
-        row_x = [set() for _ in range(num_checks)]
-        row_z = [set() for _ in range(num_checks)]
-        col_x = [set() for _ in range(self.num_data_qubits)]
-        col_z = [set() for _ in range(self.num_data_qubits)]
-
-        for stab_id, check in enumerate(checks):
-            ps, qs, _ = check
-
-            for p, q in zip(ps, qs, strict=False):
-                if p in {"X", "x"}:
-                    row_x[stab_id].add(q)
-                    col_x[q].add(stab_id)
-
-                elif p in {"Z", "z"}:
-                    row_z[stab_id].add(q)
-                    col_z[q].add(stab_id)
-
-                elif p in {"Y", "y"}:
-                    row_x[stab_id].add(q)
-                    row_z[stab_id].add(q)
-                    col_x[q].add(stab_id)
-                    col_z[q].add(stab_id)
-
-        self.check_row_x = row_x
-        self.check_row_z = row_z
-        self.check_col_x = col_x
-        self.check_col_z = col_z
-
-    def _check_commute(self) -> bool:
-        """Checks to see that all the stabilizer generators commute.
-
-        Returns:
-            Returns bool value if all the checks commute or not.
-        """
-        row_x = self.check_row_x
-        row_z = self.check_row_z
-        col_x = self.check_col_x
-        col_z = self.check_col_z
-
-        for stab_id in range(len(self.checks)):
-            anti_zs = set()
-            for q in row_x[stab_id]:
-                anti_zs ^= col_z[q]
-
-            anti_xs = set()
-            for q in row_z[stab_id]:
-                anti_xs ^= col_x[q]
-
-            anti = anti_xs ^ anti_zs
-            anti.discard(stab_id)
-
-            if anti:
-                print("\nChecks anticommute!")
-                print("\nCheck:")
-                for s in anti:
-                    print(self.checks[s][2])
-                print("\nanticommutes with:")
-                print(self.checks[stab_id][2])
-
-                msg = "Checks anticommute!"
-                raise Exception(msg)
-        return True
-
-    def compile(self) -> None:
-        """Checks commutation relations and creates a circuit to measure the checks."""
-        if self.circuit:
-            msg = "Measurement encoding-circuit has already been compiled!"
-            raise Exception(msg)
-
-        # Check the qubit ids
-        self._check_all_labels()
-
-        # Create row and column matrices.
-        self._check2rowcol()
-
-        # Checks that all the stabilizer generators (checks) commute.
-        self._check_commute()
-
-        # Create check circuits:
-        # ----------------------
-        ancilla_qubits = set()
-        qc = QuantumCircuit()
-
-        ancilla_id = sorted(self.data_qubits)[-1]
-
-        for ps, qs, _ in self.checks:
-            ancilla_id += 1
-            ancilla_qubits.add(ancilla_id)
-            qc.append("init |+>", {ancilla_id})
-
-            for p, q in zip(ps, qs, strict=False):
-                symbol = None
-
-                if p in {"X", "x"}:
-                    symbol = "CNOT"
-                elif p in {"Z", "z"}:
-                    symbol = "CZ"
-                elif p in {"Y", "y"}:
-                    symbol = "CY"
-
-                qc.append(symbol, {(ancilla_id, q)})
-
-            qc.append("measure X", {ancilla_id}, random_outcome=0)
-
-        self.ancilla_qubits = ancilla_qubits
-
-        self.circuit = qc
-
-        # Run circuits
-        # ------------
-        # Separate the checks, logical stabilizers, and ancilla stabilizers.
-        circuit = self.circuit
-        state = pc.simulators.SparseStabPy(self.num_qubits)
-        state.run_circuit(circuit)
-        self.get_info(state, verbose=False)
-        self.state = state
-
-        self._verify_checks()
-        self._check_logical_commute()
-
-    def _check_logical_commute(self) -> None:
-        logical_z_col_x = [set() for _ in range(self.num_data_qubits)]
-        logical_z_col_z = [set() for _ in range(self.num_data_qubits)]
-        logical_z_row_x = [set() for _ in range(len(self.logical_zs))]
-        logical_z_row_z = [set() for _ in range(len(self.logical_zs))]
-
-        logical_x_col_x = [set() for _ in range(self.num_data_qubits)]
-        logical_x_col_z = [set() for _ in range(self.num_data_qubits)]
-        logical_x_row_x = [set() for _ in range(len(self.logical_xs))]
-        logical_x_row_z = [set() for _ in range(len(self.logical_xs))]
-
-        for i, (ps, qs, _) in enumerate(self.logical_zs):
-            for p, q in zip(ps, qs, strict=False):
-                if p in {"X", "Y"}:
-                    logical_z_col_x[q].add(i)
-                    logical_z_row_x[i].add(q)
-
-                if p in {"Z", "Y"}:
-                    logical_z_col_z[q].add(i)
-                    logical_z_row_z[i].add(q)
-
-        for i, (ps, qs, _) in enumerate(self.logical_xs):
-            for p, q in zip(ps, qs, strict=False):
-                if p in {"X", "Y"}:
-                    logical_x_col_x[q].add(i)
-                    logical_x_row_x[i].add(q)
-
-                if p in {"Z", "Y"}:
-                    logical_x_col_z[q].add(i)
-                    logical_x_row_z[i].add(q)
-
-        for s in range(len(self.logical_zs)):
-            anti_zs = set()
-            for q in logical_z_row_x[s]:
-                anti_zs ^= logical_z_col_z[q]
-
-            anti_xs = set()
-            for q in logical_z_row_z[s]:
-                anti_xs ^= logical_z_col_x[q]
-
-            anti = anti_xs ^ anti_zs
-            anti.discard(s)
-
-            if anti:
-                print("\nLogical Zs anticommute!")
-                print("\nLogical Zs:")
-                for i in anti:
-                    print(self.logical_zs[i][2])
-                print("\nanticommutes with:")
-                print(self.logical_zs[s][2])
-
-                msg = "Logical Zs anticommute!"
-                raise Exception(msg)
-
-        for s in range(len(self.logical_xs)):
-            anti_zs = set()
-            for q in logical_x_row_x[s]:
-                anti_zs ^= logical_x_col_z[q]
-
-            anti_xs = set()
-            for q in logical_x_row_z[s]:
-                anti_xs ^= logical_x_col_x[q]
-
-            anti = anti_xs ^ anti_zs
-            anti.discard(s)
-
-            if anti:
-                print("\nLogical Xs anticommute!")
-                print("\nLogical Xs:")
-                for i in anti:
-                    print(self.logical_xs[i][2])
-                print("\nanticommutes with:")
-                print(self.logical_xs[s][2])
-
-                msg = "Logical Xs anticommute!"
-                raise Exception(msg)
-
-        # So far checked that all the logical Zs and logical Xs commute with themselves...
-        # - Next check that they commute with the stabilizers...
-        # - Then find if the there are anti-commuting pairs of logical Zs and Xs
-        # - Then search for the logical operators and refactor...
-        # This step might require switching logical Xs and Zs... Might be a bit complicated... as we can modify
-        # "logical Xs" with destabilizers and swap those... So need to do a search through all stabilizers and
-        # destabilizers and then determine if the required multiplication is valid with fix stabiliziers and whatever
-        # has been fixed for the logical operators...
-
-    def _verify_checks(self) -> bool:
-        # Stabilizers:
-        checks = []
-
-        for strings, qids, _ in self.checks:
-            check_dict = {}
-            for pauli, q in zip(strings, qids, strict=False):
-                if pauli == "X":
-                    qset = check_dict.setdefault("X", set())
-                elif pauli == "Z":
-                    qset = check_dict.setdefault("Z", set())
-                else:
-                    qset = check_dict.setdefault("Y", set())
-                qset.add(q)
-            checks.append(check_dict)
-
-        checks2 = []
-        for i in range(len(self.checks)):
-            xs = self.check_row_x[i]
-            zs = self.check_row_z[i]
-
-            stab_dict = {}
-
-            if xs - zs:
-                stab_dict["X"] = xs - zs
-
-            if zs - xs:
-                stab_dict["Z"] = zs - xs
-
-            if xs & zs:
-                stab_dict["Y"] = xs & zs
-
-            checks2.append(stab_dict)
-
-        if checks != checks2:
-            print(
-                "WARNING: PECOS didn't refactor the stabilizers into the checks supplied!",
-            )
-
-        return checks != checks2
-
-    def eval(self, *, verbose: bool = False) -> StabilizerVerificationResult:
-        """Evaluate the stabilizer code verification.
-
-        Args:
-            verbose: Whether to print detailed output during evaluation.
-
-        Returns:
-            Verification result containing success status and error details.
-        """
-        if self.circuit is None:
-            self.compile()
-
-        z, x, _, destab_strings = self.generators(verbose=verbose)
-
-        if self.dist is None:
-            self.distance(verbose=verbose)
-
-        # Stabilizers:
-        checks = []
-
-        for strings, qids, _ in self.checks:
-            check_dict = {}
-            for pauli, q in zip(strings, qids, strict=False):
-                if pauli == "X":
-                    qset = check_dict.setdefault("X", set())
-                elif pauli == "Z":
-                    qset = check_dict.setdefault("Z", set())
-                else:
-                    qset = check_dict.setdefault("Y", set())
-                qset.add(q)
-            checks.append(check_dict)
-
-        # Destabilizers:
-        destabs = []
-        for i in self.data_gens:
-            destab_dict = {}
-            for j in range(len(destab_strings[i]) - self.num_ancilla_qubits):
-                pauli = destab_strings[i][j]
-                if pauli == "X":
-                    qset = destab_dict.setdefault("X", set())
-                elif pauli == "Z":
-                    qset = destab_dict.setdefault("Z", set())
-                elif pauli == "Y":
-                    qset = destab_dict.setdefault("Y", set())
-                else:
-                    continue
-                qset.add(j - 2)
-            destabs.append(destab_dict)
-
-        output_dict = {
-            "num_datas": self.num_data_qubits,
-            "num_logical_qubits": self.num_logical_qubits(),
-            "distance": self.dist,
-            "[[n, k, d]]": f"[[{self.num_data_qubits}, {self.num_logical_qubits()}, {self.dist}]]",
-            "checks": checks,
-            "destabilizers": destabs,
-            "logical_xs": x,
-            "logical_zs": z,
-        }
-
-        self.logical_xs_reference = {}
-        self.logical_zs_reference = {}
-
-        for i, xi in enumerate(x):
-            self.logical_xs_reference["X" + str(i)] = xi
-
-        for i, zi in enumerate(z):
-            self.logical_zs_reference["Z" + str(i)] = zi
-
-        return output_dict
-
-    @property
-    def num_data_qubits(self) -> int:
-        """Get the number of data qubits.
-
-        Returns:
-            Number of data qubits in the stabilizer code.
-        """
-        return len(self.data_qubits)
-
-    @property
-    def num_ancilla_qubits(self) -> int:
-        """Get the number of ancilla qubits.
-
-        Returns:
-            Number of ancilla qubits in the stabilizer code.
-        """
-        return len(self.ancilla_qubits)
-
-    @property
-    def num_qubits(self) -> int:
-        """Get the total number of qubits.
-
-        Returns:
-            Total number of qubits (data + ancilla).
-        """
-        return len(self.data_qubits) + len(self.ancilla_qubits)
-
-    def refactor(self, state: SimulatorProtocol) -> None:
-        """Refactor the stabilizer state to match the expected generators.
-
-        Args:
-            state: Simulator state to refactor.
-        """
-        found_stab_ids = set()
-
-        refactor_things = list(self.checks)
-        refactor_things.extend(self.logical_zs)
-        # TODO: NEED TO REFACTOR THE DESTABILIZER OF LOGICAL Z TO GET THE RIGHT LOGICAL X.....
-
-        for ps, qs, _ in refactor_things:
-            xs = set()
-            zs = set()
-
-            for p, q in zip(ps, qs, strict=False):
-                if p in {"X", "x"}:
-                    xs.add(q)
-                elif p in {"Z", "z"}:
-                    zs.add(q)
-                elif p in {"Y", "y"}:
-                    xs.add(q)
-                    zs.add(q)
-
-            try:
-                found, stab_id = state.refactor(
-                    xs,
-                    zs,
-                    choose=0,
-                    protected=found_stab_ids,
-                )
-            except IndexError:
-                xonly = xs - zs
-                zonly = zs - xs
-                ys = xs & zs
-                msg = f"IndexError.\nThe stabilizer {{'X': {xonly}, 'Y': {ys}, 'Z': {zonly}}} is likely redundant!"
-                raise Exception(msg) from IndexError
-
-            found_stab_ids.add(stab_id)
-
-            if not found:
-                msg = "Could not find check:"
-                raise Exception(msg, (ps, qs))
-
-        for q in self.ancilla_qubits:
-            found, stab_id = state.refactor(
-                {q},
-                set(),
-                choose=-1,
-                protected=found_stab_ids,
-            )
-            found_stab_ids.add(stab_id)
-
-            if not found:
-                msg = f"Could not find ancilla {q}"
-                raise Exception(msg)
-
-    def get_check_ancilla(
-        self,
-    ) -> tuple[list[tuple[set[int], set[int]]], list[tuple[set[int], set[int]]]]:
-        """Get check and ancilla operator sets.
-
-        Returns:
-            Tuple containing lists of (X set, Z set) tuples for checks and ancillas.
-        """
-        check_tuples = []
-        ancilla_tuples = []
-
-        for ps, qs, _ in self.checks:
-            xs = set()
-            zs = set()
-
-            for p, q in zip(ps, qs, strict=False):
-                if p in {"X", "x"}:
-                    xs.add(q)
-                elif p in {"Z", "z"}:
-                    zs.add(q)
-                elif p in {"Y", "y"}:
-                    xs.add(q)
-                    zs.add(q)
-
-            check_tuples.append((xs, zs))
-
-        ancilla_tuples.extend(({q}, set()) for q in self.ancilla_qubits)
-
-        return check_tuples, ancilla_tuples
-
-    def get_info(
-        self,
-        state: SimulatorProtocol,
-        stop_search: int = 1000,
-        *,
-        verbose: bool = True,
-        print_y: bool = False,
-    ) -> tuple[
-        list[dict[str, set[int]]],
-        list[dict[str, set[int]]],
-        list[str],
-        list[str],
-    ]:
-        """Get stabilizer information from the quantum state.
-
-        Args:
-            state: Simulator state to analyze.
-            stop_search: Maximum number of refactoring attempts.
-            verbose: Whether to print detailed information.
-            print_y: Whether to include Y operators in output.
-
-        Returns:
-            Tuple of logical Z operators, logical X operators, stabilizer strings, destabilizer strings.
-        """
-        if self.circuit is None:
-            return Exception("Must run `compile()` first!")
-
-        self.refactor(state)
-        stab_strs, destab_strs = state.print_stabs(
-            verbose=False,
-            print_y=print_y,
-            print_destabs=True,
-        )
-
-        num_ancillas = len(self.ancilla_qubits)
-
-        num_logical = self.num_logical_qubits()
-        num_checks = len(self.checks)
-
-        if verbose:
-            print(f"Number of data qubits: {self.num_data_qubits}")
-            print(f"Number of checks: {num_checks}")
-            print(f"Number of logical qubits: {num_logical}")
-
-        check_tuples, ancilla_tuples = self.get_check_ancilla()
-        # determine the gen_id of the checks and logicals
-        check_gens = []
-        logical_gens = []
-        ancilla_gens = []
-
-        missing_checks = list(check_tuples)
-        missing_ancillas = list(ancilla_tuples)
-        notmatched_gens = list(range(state.num_qubits))
-
-        found_all = False
-        search_count = 0
-        while not found_all:
-            if verbose:
-                print("----")
-            for g, gtuple in enumerate(
-                zip(state.stabs.row_x, state.stabs.row_z, strict=False),
-            ):
-                if gtuple in missing_checks:
-                    missing_checks.remove(gtuple)
-                    try:
-                        notmatched_gens.remove(g)
-                    except ValueError:
-                        msg = f"list.remove(x): x not in list.\nThe stabilizer {gtuple!s} is likely redundant!"
-                        raise Exception(
-                            msg,
-                        ) from ValueError
-
-                    check_gens.append(g)
-                elif gtuple in missing_ancillas:
-                    missing_ancillas.remove(gtuple)
-                    notmatched_gens.remove(g)
-                    ancilla_gens.append(g)
-
-            if len(notmatched_gens) == num_logical:
-                logical_gens = notmatched_gens
-                found_all = True
-            else:
-                for xs, zs in missing_checks:
-                    state.refactor(xs, zs, choose=0, prefer=notmatched_gens)
-                state.print_stabs(verbose=False, print_y=print_y, print_destabs=True)
-
-            if search_count == stop_search:
-                msg = "Can not refactor properly!"
-                raise Exception(msg)
-            search_count += 1
-
-        self.data_gens = set(check_gens)
-        self.ancilla_gens = set(ancilla_gens)
-        self.logical_gens = set(logical_gens)
-
-        if verbose:
-            if len(check_gens) != num_checks:
-                print("Found:", check_gens)
-                print("Want:", check_tuples)
-                msg = f"Did not find the correct number of stabilizer generators. {len(check_gens)}/{num_checks}"
-                raise Exception(msg)
-
-            if len(logical_gens) != num_logical:
-                print("Found:", logical_gens)
-                msg = f"Did not find the correct number of logical generators. {len(logical_gens)}/{num_logical}"
-                raise Exception(msg)
-
-            print("\nStabilizer generators:")
-            for gen in check_gens:
-                print(stab_strs[gen][: len(stab_strs[gen]) - num_ancillas])
-
-            print("\nDestabilizer generators:")
-            for gen in check_gens:
-                print(destab_strs[gen][: len(destab_strs[gen]) - num_ancillas])
-
-            print("\nLogical operators:")
-
-            for i, gen in enumerate(logical_gens):
-                print(f"\n. Logical Z #{i + 1!s}:")
-                print(stab_strs[gen][: len(stab_strs[gen]) - num_ancillas])
-                print(f". Logical X #{i + 1!s}:")
-                print(destab_strs[gen][: len(destab_strs[gen]) - num_ancillas])
-
-        logical_z_strings = []
-        logical_x_strings = []
-
-        for gen in logical_gens:
-            z_string = stab_strs[gen][: len(stab_strs[gen]) - num_ancillas]
-            x_string = destab_strs[gen][: len(destab_strs[gen]) - num_ancillas]
-
-            check_dict = {}
-            for q, pauli in enumerate(z_string):
-                if pauli == "X":
-                    qset = check_dict.setdefault("X", set())
-                elif pauli == "Z":
-                    qset = check_dict.setdefault("Z", set())
-                elif pauli == "Y":
-                    qset = check_dict.setdefault("Y", set())
-                else:
-                    continue
-                qset.add(q - 2)
-            logical_z_strings.append(check_dict)
-
-            check_dict = {}
-            for q, pauli in enumerate(x_string):
-                if pauli == "X":
-                    qset = check_dict.setdefault("X", set())
-                elif pauli == "Z":
-                    qset = check_dict.setdefault("Z", set())
-                elif pauli == "Y":
-                    qset = check_dict.setdefault("Y", set())
-                else:
-                    continue
-
-                qset.add(q - 2)
-            logical_x_strings.append(check_dict)
-
-        return logical_z_strings, logical_x_strings, stab_strs, destab_strs
-
-    def distance(
-        self,
-        *,
-        css: bool = False,
-        verbose: bool = True,
-    ) -> tuple[set[int], set[int]] | None:
-        """Checks the distance of the code."""
-        if self.circuit is None:
-            msg = "Must compile circuits first!"
-            raise Exception(msg)
-
-        qudit_set = self.data_qubits
-
-        state = self.state
-        found = self._dist_mode_smallest(state, qudit_set, css=css, verbose=verbose)
-
-        if verbose and found:
-            xs, zs = found
-            distance = len(xs | zs)
-
-            print(
-                f"\nThis is a [[{self.num_data_qubits}, {self.num_logical_qubits()}, {distance}]] code.",
-            )
-
-        if not found:
-            print(
-                "No logical errors found... Checks might describe a stabilizer state.",
-            )
-            return None
-        xs, zs = found
-        self.dist = len(xs | zs)
-        return found
-
-    def _dist_mode_smallest(
-        self,
-        state: SimulatorProtocol,
-        qudit_set: set[int],
-        *,
-        css: bool = False,
-        verbose: bool = True,
-        start_len: int | None = None,
-        end_len: int | None = None,
-        list_ops: bool = False,
-    ) -> Generator[tuple[set, set], None, None]:
-        """Determine if a logical error can be found by starting with the smallest weight errors.
-
-        Args:
-        ----
-            state: The quantum state to check for logical errors.
-            qudit_set: Set of qudit indices to consider for errors.
-            css: If True, only checks CSS (Calderbank-Shor-Steane) type errors.
-            verbose: If True, prints progress and found logical operators.
-            start_len: Starting weight of errors to check (default: 1).
-            end_len: Maximum weight of errors to check (default: number of qudits).
-            list_ops: If True, returns a list of all logical operators found.
-
-        """
-        ops = []
-
-        if start_len is None:
-            start_len = 1
-
-        if end_len is None:
-            end_len = len(qudit_set)
-
-        for lenq in range(start_len, end_len + 1):
-            if verbose:
-                print(f"Checking Paulis of weight {lenq}...")
-
-            for xs, zs in self.gen_errors(qudit_set, lenq, lenq, css=css):
-                if self._is_logical_error(state, xs, zs):
-                    if verbose:
-                        print(f"Logical operator found: Xs - {xs} Zs - {zs}")
-
-                    if list_ops:
-                        ops.append({"X": xs, "Z": zs})
-                    else:
-                        return xs, zs
-
-        return ops
-
-    def gen_errors(
-        self,
-        qubits: set[int] | Sequence[int],
-        min_errors: int = 1,
-        *,
-        max_errors: bool | int = False,
-        css: bool = False,
-    ) -> Generator[tuple[set[int], set[int]], None, None]:
-        """Generate error patterns for testing stabilizer codes.
-
-        Args:
-        ----
-            qubits (set of int): Set of qubit indices to generate errors on.
-            min_errors (int): Minimum number of errors to generate.
-            max_errors (bool, int): Maximum number of errors to generate. False for no limit.
-            css (bool): If True, generate only CSS-compatible errors (X and Z only).
-
-        """
-        paulis = ("X", "Z", "Y")
-
-        num_qubits = len(qubits)
-
-        for i in range(min_errors, num_qubits + 1):
-            if max_errors and i > max_errors:
-                break
-
-            xs = next(product(("X",), repeat=i))
-            zs = next(product(("Z",), repeat=i))
-
-            xzs = [xs, zs]
-
-            for b in combinations(qubits, i):
-                for ps in xzs:
-                    x_set = set()
-                    z_set = set()
-                    for p, q in zip(ps, b, strict=False):
-                        if p == "X":
-                            x_set.add(q)
-                        else:
-                            z_set.add(q)
-                    yield x_set, z_set
-
-            if not css:
-                for a in product(paulis, repeat=i):
-                    if a in {xs, zs}:
-                        continue
-
-                    for b in combinations(qubits, i):
-                        x_set = set()
-                        z_set = set()
-                        for p, q in zip(a, b, strict=False):
-                            if p == "X":
-                                x_set.add(q)
-                            elif p == "Z":
-                                z_set.add(q)
-                            else:
-                                x_set.add(q)
-                                z_set.add(q)
-                        yield x_set, z_set
-
-    def _is_logical_error(
-        self,
-        state: SimulatorProtocol,
-        xs: set[int],
-        zs: set[int],
-    ) -> bool:
-        # A trivial error anticommutes with the checks. (Might or might not anticommute with the logical stabilizers)
-        # A logical error commutes with the checks and is not a product of checks.
-
-        # Does the error anticommute with the checks?
-        x_anticoms = set()
-        z_anticoms = set()
-        for q in xs:
-            x_anticoms ^= state.stabs.col_z[q]
-
-        for q in zs:
-            z_anticoms ^= state.stabs.col_x[q]
-
-        anticoms = x_anticoms ^ z_anticoms
-        anticom_logical_zs = self.logical_gens & anticoms
-        anticoms -= self.logical_gens
-
-        if anticoms:
-            return False
-        if anticom_logical_zs:
-            # So the error commutes with all the stabilizers
-            # Did it anticommute with any logical Z operations? If so... It is a product of logical Xs!
-            # (and possibly other things)
-            return True
-        # Let's see if the error anticommuted with any logical X operators:
-
-        x_anticoms_destabs = set()
-        z_anticoms_destabs = set()
-
-        for q in xs:
-            x_anticoms_destabs ^= state.destabs.col_z[q]
-
-        for q in zs:
-            z_anticoms_destabs ^= state.destabs.col_x[q]
-
-        anticoms_destabs = x_anticoms_destabs ^ z_anticoms_destabs
-        anticom_logical_xs = self.logical_gens & anticoms_destabs
-
-        # The error is a product of logical Zs
-        return bool(anticom_logical_xs)
-
-    def shortest_logicals(
-        self,
-        start_weight: int | None = None,
-        delta: int = 0,
-        *,
-        verbose: bool = True,
-        css: bool = False,
-    ) -> tuple[
-        list[LogicalOpInfo],
-        dict[str, dict[str, set[int]]],
-        dict[str, dict[str, set[int]]],
-    ]:
-        """Find the shortest logical operators.
-
-        Args:
-            start_weight (int): Weight of operators to begin searching.
-            delta (int): Method will look for all logical ops with weight =< minimum weight + `delta`.
-            verbose (bool): If True, print progress information during the search.
-            css (bool): If True, restrict search to CSS-compatible operators (X and Z only).
-
-        Returns:
-        -------
-            Dictionary of logical ops...
-
-        """
-        # if not self.logical_xs_reference and not self.logical_zs_reference:
-
-        if start_weight is None:
-            start_weight = self.dist if self.dist is not None else 1
-
-        end_weight = start_weight + delta
-
-        if self.circuit is None:
-            msg = "Must compile circuits first!"
-            raise Exception(msg)
-
-        qudit_set = self.data_qubits
-
-        end_weight = min(end_weight, len(qudit_set))
-
-        state = self.state
-        found = self._dist_mode_smallest(
-            state,
-            qudit_set,
-            css=css,
-            verbose=False,
-            start_len=start_weight,
-            end_len=end_weight,
-            list_ops=True,
-        )
-
-        xs_labels = sorted(self.logical_xs_reference.keys())
-        zs_labels = sorted(self.logical_zs_reference.keys())
-
-        oplist = []
-
-        if found:
-            for paulis in found:
-                op_product = []
-                for xi, op_label in enumerate(xs_labels):
-                    if self.op_anticommute(paulis, self.logical_xs_reference[op_label]):
-                        op_product.append(zs_labels[xi])
-
-                for zi, op_label in enumerate(zs_labels):
-                    if self.op_anticommute(paulis, self.logical_zs_reference[op_label]):
-                        op_product.append(xs_labels[zi])
-
-                op_product = sorted(op_product)
-
-                oplist.append(
-                    {
-                        "X": paulis["X"],
-                        "Z": paulis["Z"],
-                        "equiv_ops": tuple(op_product),
-                    },
-                )
-
-        if verbose:
-            print("Reference Logical Operators:")
-            print("\nLogical Xs:")
-            for op_label in xs_labels:
-                op = self.logical_xs_reference[op_label]
-                print(op_label, op)
-            print("\nLogical Zs:")
-            for op_label in zs_labels:
-                op = self.logical_zs_reference[op_label]
-                print(op_label, op)
-
-            print("\nLogical Ops Found:\n")
-            for foundop in oplist:
-                print(
-                    "X - {} Z - {} Equiv Ops - {}".format(
-                        foundop["X"],
-                        foundop["Z"],
-                        foundop["equiv_ops"],
-                    ),
-                )
-
-        return oplist, self.logical_xs_reference, self.logical_zs_reference
-
-    @staticmethod
-    def op_anticommute(op1: dict[str, set[int]], op2: dict[str, set[int]]) -> bool:
-        """Check if two Pauli operators anticommute.
-
-        Args:
-            op1: First Pauli operator as dictionary with X, Y, Z keys and qubit sets.
-            op2: Second Pauli operator as dictionary with X, Y, Z keys and qubit sets.
-
-        Returns:
-            True if the operators anticommute, False otherwise.
-        """
-        return bool(
-            (len(op1.get("X", set()) & op2.get("Z", set())) + len(op2.get("X", set()) & op1.get("Z", set()))) % 2,
-        )
diff --git a/python/quantum-pecos/tests/pecos/test_stabilizer_code_spec_bindings.py b/python/quantum-pecos/tests/pecos/test_stabilizer_code_spec_bindings.py
new file mode 100644
index 000000000..0d3b6d539
--- /dev/null
+++ b/python/quantum-pecos/tests/pecos/test_stabilizer_code_spec_bindings.py
@@ -0,0 +1,315 @@
+# Copyright 2026 The PECOS Developers
+#
+# Licensed under the Apache License, Version 2.0
+
+from collections.abc import Callable
+
+import numpy as np
+import pytest
+from pecos.quantum import (
+    DistanceResult,
+    LogicalOperatorInfo,
+    ParityCheckMatrix,
+    PauliString,
+    StabilizerCode,
+    StabilizerCodeSpec,
+    SymplecticMatrix,
+    Zs,
+)
+
+_HAMMING_H = [
+    [1, 0, 1, 0, 1, 0, 1],
+    [0, 1, 1, 0, 0, 1, 1],
+    [0, 0, 0, 1, 1, 1, 1],
+]
+
+
+def _five_qubit_spec() -> StabilizerCodeSpec:
+    return StabilizerCodeSpec(
+        5,
+        [
+            PauliString.from_dense_str("XZZXI"),
+            PauliString.from_dense_str("IXZZX"),
+            PauliString.from_dense_str("XIXZZ"),
+            PauliString.from_dense_str("ZXIXZ"),
+        ],
+        [PauliString.from_dense_str("ZZZZZ")],
+        [PauliString.from_dense_str("XXXXX")],
+    )
+
+
+def _repetition_spec() -> StabilizerCodeSpec:
+    return StabilizerCodeSpec(
+        3,
+        [
+            PauliString.from_dense_str("ZZI"),
+            PauliString.from_dense_str("IZZ"),
+        ],
+        [PauliString.from_dense_str("ZZZ")],
+        [PauliString.from_dense_str("XXX")],
+    )
+
+
+def test_five_qubit_hand_built_spec_finds_genuine_weight_three_logical() -> None:
+    spec = _five_qubit_spec()
+    spec.verify()
+
+    assert spec.num_qubits == 5
+    assert spec.num_logical_qubits == 1
+    assert len(spec.stabilizers) == 4
+    assert spec.logical_zs == [PauliString.from_dense_str("ZZZZZ")]
+    assert spec.logical_xs == [PauliString.from_dense_str("XXXXX")]
+
+    result = spec.distance()
+
+    assert isinstance(result, DistanceResult)
+    assert result.distance == 3
+    assert result.min_weight_operator.weight() == 3
+    assert "distance=" in repr(result)
+    assert str(result.min_weight_operator) in repr(result)
+    assert StabilizerCode.five_qubit().syndrome(result.min_weight_operator) == [False] * 4
+
+
+def test_steane_css_and_general_searches_both_find_distance_three() -> None:
+    spec = StabilizerCodeSpec.from_stabilizer_code(StabilizerCode.steane())
+
+    general = spec.distance()
+    css = spec.distance(css=True)
+
+    assert general is not None
+    assert css is not None
+    assert general.distance == css.distance == 3
+
+
+def test_repetition_min_weight_logicals_expose_equivalence_information() -> None:
+    spec = _repetition_spec()
+    result = spec.distance()
+    logicals = spec.min_weight_logicals()
+
+    assert result is not None
+    assert result.distance == 1
+    assert logicals
+    assert all(isinstance(info, LogicalOperatorInfo) for info in logicals)
+    assert all(info.weight == info.operator.weight() == 1 for info in logicals)
+    assert all(info.equivalent_logicals == [("Z", 0)] for info in logicals)
+    assert {info.equivalence_string() for info in logicals} <= {"X0", "Z0"}
+    assert all(str(info.operator) in repr(info) for info in logicals)
+    assert all(info.equivalence_string() in repr(info) for info in logicals)
+
+
+def test_verbose_distance_searches_write_progress_to_stderr(
+    capfd: pytest.CaptureFixture[str],
+) -> None:
+    spec = _repetition_spec()
+
+    assert spec.distance(verbose=True) is not None
+    assert "Checking weight" in capfd.readouterr().err
+
+    assert spec.min_weight_logicals(verbose=True)
+    assert "Checking weight" in capfd.readouterr().err
+
+
+def test_from_stabilizer_code_steane_round_trip_finds_distance_three() -> None:
+    code = StabilizerCode.steane()
+    spec = StabilizerCodeSpec.from_stabilizer_code(code)
+
+    spec.verify()
+    result = spec.distance()
+
+    assert spec.num_qubits == code.num_qubits()
+    assert spec.num_logical_qubits == code.num_logical_qubits()
+    assert result is not None
+    assert result.distance == 3
+
+
+def test_max_weight_below_true_distance_returns_no_results() -> None:
+    spec = _five_qubit_spec()
+
+    assert spec.distance(max_weight=2) is None
+    assert spec.min_weight_logicals(max_weight=2) == []
+    assert spec.shortest_logicals(delta=1, max_weight=2) == []
+
+
+def test_five_qubit_shortest_logicals_include_requested_weight_range() -> None:
+    spec = _five_qubit_spec()
+    minimum = spec.min_weight_logicals()
+    delta_zero = spec.shortest_logicals()
+    delta_one = spec.shortest_logicals(delta=1)
+    delta_two = spec.shortest_logicals(delta=2)
+
+    minimum_operators = [info.operator for info in minimum]
+    delta_zero_operators = [info.operator for info in delta_zero]
+    delta_one_operators = [info.operator for info in delta_one]
+    delta_two_operators = [info.operator for info in delta_two]
+
+    assert len(minimum_operators) == 30
+    assert delta_zero_operators == minimum_operators
+    assert delta_one_operators == minimum_operators
+    assert delta_two_operators[: len(minimum_operators)] == minimum_operators
+    assert sum(info.weight == 5 for info in delta_two) == 18
+    assert {info.weight for info in delta_two} == {3, 5}
+    assert all(StabilizerCode.five_qubit().syndrome(info.operator) == [False] * 4 for info in delta_two)
+
+
+def test_repetition_shortest_logicals_exclude_non_logicals_in_range() -> None:
+    logicals = _repetition_spec().shortest_logicals(delta=1)
+
+    assert logicals
+    assert {info.weight for info in logicals} == {1}
+
+
+@pytest.mark.parametrize(
+    "constructor",
+    [StabilizerCode.steane, StabilizerCode.five_qubit, StabilizerCode.shor],
+)
+def test_spec_distance_matches_stabilizer_code_oracle(
+    constructor: Callable[[], StabilizerCode],
+) -> None:
+    code = constructor()
+    result = StabilizerCodeSpec.from_stabilizer_code(code).distance()
+
+    assert result is not None
+    assert result.distance == code.distance()
+
+
+def test_noncommuting_stabilizers_raise_python_exception_on_verify() -> None:
+    spec = StabilizerCodeSpec(
+        1,
+        [PauliString.from_dense_str("X"), PauliString.from_dense_str("Z")],
+        [],
+        [],
+    )
+
+    with pytest.raises(ValueError, match="Stabilizer generators 0 and 1 anticommute"):
+        spec.verify()
+
+
+def test_constructor_errors_are_python_exceptions() -> None:
+    with pytest.raises(ValueError, match="Number of logical X and Z operators must match"):
+        StabilizerCodeSpec(
+            1,
+            [],
+            [PauliString.from_dense_str("Z")],
+            [],
+        )
+
+
+@pytest.mark.parametrize(
+    "rows",
+    [
+        _HAMMING_H,
+        np.asarray(_HAMMING_H, dtype=np.int64),
+        np.asarray(_HAMMING_H, dtype=np.uint8),
+    ],
+    ids=["lists", "numpy-int64", "numpy-uint8"],
+)
+def test_parity_check_matrix_builds_steane_code_from_dense_inputs(rows: object) -> None:
+    matrix = ParityCheckMatrix(rows)
+    builder = StabilizerCodeSpec.builder(7)
+    builder.checks_from_css(matrix, matrix)
+    spec = builder.build_with_discovered_logicals()
+    result = spec.distance()
+
+    assert matrix.num_checks() == matrix.rank() == 3
+    assert matrix.num_qubits() == 7
+    assert matrix.rows() == _HAMMING_H
+    assert len(matrix.to_x_stabilizers()) == len(matrix.to_z_stabilizers()) == 3
+    assert repr(matrix) == "ParityCheckMatrix(shape=(3, 7))"
+    assert result is not None
+    assert result.distance == StabilizerCode.steane().distance() == 3
+
+
+def test_symplectic_matrix_builds_five_qubit_code() -> None:
+    rows = [
+        [1, 0, 0, 1, 0, 0, 1, 1, 0, 0],
+        [0, 1, 0, 0, 1, 0, 0, 1, 1, 0],
+        [1, 0, 1, 0, 0, 0, 0, 0, 1, 1],
+        [0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
+    ]
+    matrix = SymplecticMatrix.from_dense(rows)
+    builder = StabilizerCodeSpec.builder(5)
+    builder.checks_from_symplectic(matrix)
+    spec = builder.build_with_discovered_logicals()
+    result = spec.distance()
+
+    assert matrix.num_rows() == matrix.rank() == 4
+    assert matrix.num_qubits() == 5
+    assert matrix.rows() == rows
+    assert matrix.x_block() == [row[:5] for row in rows]
+    assert matrix.z_block() == [row[5:] for row in rows]
+    assert matrix.to_positive_paulis() == spec.stabilizers
+    assert repr(matrix) == "SymplecticMatrix(shape=(4, 5))"
+    assert result is not None
+    assert result.distance == 3
+
+
+def test_code_matrix_entry_and_shape_errors_are_value_errors() -> None:
+    with pytest.raises(ValueError, match=r"row 0.*value 2"):
+        ParityCheckMatrix([[0, 2]])
+    with pytest.raises(ValueError, match=r"row 0.*value -1"):
+        SymplecticMatrix([[0, -1]])
+    with pytest.raises(ValueError, match=r"row 1.*columns"):
+        ParityCheckMatrix([[1, 0], [1]])
+    with pytest.raises(ValueError, match="even column count"):
+        SymplecticMatrix.from_dense([[1, 0, 1]])
+
+
+def test_code_matrix_builder_validation_errors_preserve_diagnostics() -> None:
+    width_mismatch = StabilizerCodeSpec.builder(2)
+    with pytest.raises(ValueError, match=r"3 qubits, expected 2"):
+        width_mismatch.checks_from_css(
+            ParityCheckMatrix([[1, 0, 0]]),
+            ParityCheckMatrix.zeros(0, 2),
+        )
+
+    nonorthogonal = StabilizerCodeSpec.builder(2)
+    with pytest.raises(ValueError, match=r"X row 0 and Z row 0"):
+        nonorthogonal.checks_from_css(
+            ParityCheckMatrix([[1, 0]]),
+            ParityCheckMatrix([[1, 0]]),
+        )
+
+    symplectic_width_mismatch = StabilizerCodeSpec.builder(2)
+    with pytest.raises(ValueError, match=r"3 qubits, expected 2"):
+        symplectic_width_mismatch.checks_from_symplectic(SymplecticMatrix.zeros(0, 3))
+
+
+def test_spec_constructor_rejects_dependent_stabilizers() -> None:
+    with pytest.raises(ValueError, match=r"rank 2, count 3"):
+        StabilizerCodeSpec(
+            3,
+            [Zs([0, 1]), Zs([1, 2]), Zs([0, 2])],
+            [],
+            [],
+        )
+
+
+def test_zero_row_parity_check_matrix_preserves_width_for_x_only_code() -> None:
+    x_checks = ParityCheckMatrix([[1, 1]])
+    z_checks = ParityCheckMatrix.zeros(0, 2)
+    builder = StabilizerCodeSpec.builder(2)
+    builder.checks_from_css(x_checks, z_checks)
+    spec = builder.build_with_discovered_logicals()
+
+    assert z_checks.num_checks() == 0
+    assert z_checks.num_qubits() == 2
+    assert z_checks.rows() == []
+    assert spec.num_logical_qubits == 1
+    assert spec.stabilizers == x_checks.to_x_stabilizers()
+
+
+def test_quantum_namespace_exports_distance_search_types() -> None:
+    import pecos.quantum as quantum
+
+    assert quantum.StabilizerCodeSpec is StabilizerCodeSpec
+    assert quantum.DistanceResult is DistanceResult
+    assert quantum.LogicalOperatorInfo is LogicalOperatorInfo
+    assert quantum.ParityCheckMatrix is ParityCheckMatrix
+    assert quantum.SymplecticMatrix is SymplecticMatrix
+    assert {
+        "StabilizerCodeSpec",
+        "DistanceResult",
+        "LogicalOperatorInfo",
+        "ParityCheckMatrix",
+        "SymplecticMatrix",
+    } <= set(quantum.__all__)
diff --git a/python/quantum-pecos/tests/pecos/test_stabilizer_code_verification_bindings.py b/python/quantum-pecos/tests/pecos/test_stabilizer_code_verification_bindings.py
new file mode 100644
index 000000000..e57fee001
--- /dev/null
+++ b/python/quantum-pecos/tests/pecos/test_stabilizer_code_verification_bindings.py
@@ -0,0 +1,139 @@
+# Copyright 2026 The PECOS Developers
+#
+# Licensed under the Apache License, Version 2.0
+
+import re
+
+import pytest
+from pecos.quantum import (
+    PauliString,
+    StabilizerCodeSpec,
+    StabilizerCodeSpecBuilder,
+    X,
+    Xs,
+    Y,
+    Ys,
+    Z,
+    Zs,
+    pauli_string,
+)
+
+
+def _original_checks() -> list[PauliString]:
+    return [
+        Xs([3, 4, 7, 8]),
+        Xs([5, 6, 7, 9]),
+        Zs([2, 4, 5, 7]),
+        Zs([7, 8, 9]),
+        Zs([0, 1]) * Y(2),
+        pauli_string("X0 X2 Z3 Y4"),
+        pauli_string("X1 X2 Z6 Y5"),
+    ]
+
+
+def _fixed_checks() -> list[PauliString]:
+    checks = _original_checks()
+    checks[4] = Zs([0, 1, 2])
+    return checks
+
+
+def _final_checks() -> list[PauliString]:
+    return [
+        Zs([2, 4, 5, 7]),
+        Xs([3, 4, 7, 8]),
+        Xs([5, 6, 7, 9]),
+        pauli_string("X0 X2 Z3 Y4"),
+        pauli_string("X1 X2 Z6 Y5"),
+        Zs([0, 1, 2]),
+        Xs([0, 1]),
+        Zs([3, 8]),
+        Zs([6, 9]),
+    ]
+
+
+def _builder_with_checks(checks: list[PauliString]) -> StabilizerCodeSpecBuilder:
+    builder = StabilizerCodeSpec.builder(10)
+    for check in checks:
+        builder.check(check)
+    return builder
+
+
+def test_original_doc_checks_report_a_real_anticommuting_pair() -> None:
+    checks = _original_checks()
+
+    with pytest.raises(
+        ValueError,
+        match=r"Stabilizer generators \d+ and \d+ anticommute",
+    ) as exc_info:
+        _builder_with_checks(checks).build_verified()
+
+    pair = re.search(r"generators (\d+) and (\d+) anticommute", str(exc_info.value))
+    assert pair is not None
+    first, second = (int(index) for index in pair.groups())
+    assert checks[first].anticommutes_with(checks[second])
+
+    with pytest.raises(ValueError, match="Stabilizers do not all commute with each other"):
+        _builder_with_checks(checks).build_with_discovered_logicals()
+
+
+def test_fixed_doc_checks_build_a_distance_two_code() -> None:
+    spec = _builder_with_checks(_fixed_checks()).build_with_discovered_logicals()
+    result = spec.distance()
+
+    assert spec.num_logical_qubits == 3
+    assert result is not None
+    assert result.distance == 2
+    assert result.min_weight_operator.weight() == 2
+
+
+def test_final_doc_checks_build_a_distance_three_code() -> None:
+    spec = _builder_with_checks(_final_checks()).build_with_discovered_logicals()
+    result = spec.distance()
+
+    assert spec.num_logical_qubits == 1
+    assert len(spec.destabilizers) == 9
+    assert result is not None
+    assert result.distance == 3
+
+
+def test_multi_qubit_pauli_helpers_match_single_qubit_composition() -> None:
+    assert Xs([0, 2, 5]) == X(0) & X(2) & X(5)
+    assert Ys((1, 3)) == Y(1) & Y(3)
+    assert Zs([]) == PauliString.I()
+    assert Zs(range(3)) == Z(0) & Z(1) & Z(2)
+
+
+def test_builder_is_consumed_and_validates_logical_counts() -> None:
+    builder = StabilizerCodeSpec.builder(1)
+    spec = builder.build()
+
+    assert spec.num_logical_qubits == 1
+    with pytest.raises(RuntimeError, match="already been consumed"):
+        builder.check(Z(0))
+
+    builder = StabilizerCodeSpec.builder(1)
+    builder.logical_z(Z(0))
+    builder.logical_x(X(0))
+    spec = builder.build_verified()
+
+    assert spec.num_logical_qubits == 1
+    assert spec.logical_zs == [Z(0)]
+    assert spec.logical_xs == [X(0)]
+
+    mismatched = StabilizerCodeSpec.builder(1)
+    mismatched.logical_x(X(0))
+    with pytest.raises(ValueError, match="Number of logical X and Z operators must match"):
+        mismatched.build()
+
+
+def test_string_summary_lists_final_code_generators() -> None:
+    spec = _builder_with_checks(_final_checks()).build_with_discovered_logicals()
+    summary = str(spec)
+
+    assert "[[10, 1]]" in summary
+    assert "Stabilizer generators:" in summary
+    assert "Destabilizer generators:" in summary
+    assert "Z1:" in summary
+    assert "X1:" in summary
+    assert repr(spec) == "StabilizerCodeSpec([[10, 1]])"
+    assert all(operator.to_dense_str(10) in summary for operator in spec.stabilizers)