e831e513-b95a-4211-b11a-d686a4b91dfb/tutorial-video

"""Quantum error-correction utilities inspired by
Calderbank, Rains, Shor and Sloane (1996).

This module provides lightweight, dependency-free tools that capture the
algebraic framework described in *Quantum Error Correction and Orthogonal
Geometry*.  The focus is on the extraspecial 2-group of real Pauli
operators, the associated quadratic and symplectic forms, and the
construction of stabilizer (additive) quantum codes from totally isotropic
subspaces.

The implementation is intentionally self-contained and aims to be both
educational and practical:

* The :class:`Pauli` class models elements of the extraspecial 2-group ``E``
  using their binary (symplectic) representation ``(a | b)``.
* Clifford generators listed in the paper appear here as explicit binary
  symplectic transformations that map Pauli operators to Pauli operators.
* The :class:`StabilizerCode` class encapsulates stabilizer codes arising
  from totally isotropic subspaces, exposes their parameters ``[[n, k, d]]``,
  and can estimate the distance by exhaustive search.
* Ready-to-use constructors replicate the canonical five-qubit code
  (Example 1) and the CSS construction based on nested classical codes
  (Example 2).

Example
-------

>>> code = five_qubit_code()
>>> code.parameters  # doctest: +NORMALIZE_WHITESPACE
{'n': 5, 'k': 1, 'd': 3}
>>> code.is_totally_isotropic()
True
>>> logical_x, logical_z = code.logicals
>>> logical_x * logical_z == logical_z * logical_x
False

The example mirrors the group-theoretic presentation in the paper: the
stabilizer is a 4-dimensional totally singular subspace of ``E`` whose
simultaneous eigenspace encodes one protected qubit.
"""

from __future__ import annotations

from dataclasses import dataclass
from functools import reduce
from itertools import combinations, product
from typing import Iterable, List, Optional, Sequence, Tuple, Union

BinaryVector = Tuple[int, ...]

__all__ = [
    "Pauli",
    "symplectic_inner",
    "quadratic_form",
    "is_totally_isotropic",
    "StabilizerCode",
    "css_code",
    "five_qubit_code",
    "steane_code",
]


# ---------------------------------------------------------------------------
# Binary linear algebra helpers
# ---------------------------------------------------------------------------

def _ensure_binary_vector(vector: Sequence[int]) -> BinaryVector:
    """Validate and convert a sequence into a tuple of bits."""
    bits = tuple(int(v) & 1 for v in vector)
    if len(bits) != len(vector):
        raise AssertionError("Internal error while normalising binary vector.")
    return bits


def _parse_binary_string(bitstring: str) -> BinaryVector:
    """Parse a 0/1 string into a binary vector."""
    cleaned = bitstring.replace(" ", "")
    if not cleaned:
        return tuple()
    if any(ch not in "01" for ch in cleaned):
        raise ValueError(f"Invalid binary string: {bitstring!r}")
    return tuple(int(ch) for ch in cleaned)


def gf2_rank(matrix: Sequence[BinaryVector]) -> int:
    """Return the rank of a binary matrix using Gauss elimination in GF(2)."""
    if not matrix:
        return 0
    rows = [list(row) for row in matrix]
    n_rows, n_cols = len(rows), len(rows[0])
    rank = 0
    for col in range(n_cols):
        pivot = None
        for row in range(rank, n_rows):
            if rows[row][col]:
                pivot = row
                break
        if pivot is None:
            continue
        rows[rank], rows[pivot] = rows[pivot], rows[rank]
        for row in range(n_rows):
            if row != rank and rows[row][col]:
                rows[row] = [
                    (rows[row][c] ^ rows[rank][c]) for c in range(n_cols)
                ]
        rank += 1
        if rank == n_rows:
            break
    return rank


def gf2_row_space_basis(matrix: Sequence[BinaryVector]) -> Tuple[BinaryVector, ...]:
    """Return a basis of the row space in reduced row-echelon form."""
    if not matrix:
        return tuple()
    rows = [list(row) for row in matrix]
    n_rows, n_cols = len(rows), len(rows[0])
    rank = 0
    for col in range(n_cols):
        pivot = None
        for row in range(rank, n_rows):
            if rows[row][col]:
                pivot = row
                break
        if pivot is None:
            continue
        rows[rank], rows[pivot] = rows[pivot], rows[rank]
        for row in range(n_rows):
            if row != rank and rows[row][col]:
                rows[row] = [
                    (rows[row][c] ^ rows[rank][c]) for c in range(n_cols)
                ]
        rank += 1
        if rank == n_rows:
            break
    return tuple(tuple(row) for row in rows[:rank])


def vector_in_row_space(vector: BinaryVector, basis: Sequence[BinaryVector]) -> bool:
    """Decide whether ``vector`` belongs to the row space spanned by ``basis``."""
    if not basis:
        return not any(vector)
    work = list(vector)
    basis_rows = [list(row) for row in basis]
    pivots = []
    for row in basis_rows:
        pivot = next((idx for idx, entry in enumerate(row) if entry), None)
        if pivot is not None:
            pivots.append(pivot)
    for row, pivot in zip(basis_rows, pivots):
        if work[pivot]:
            work = [w ^ r for w, r in zip(work, row)]
    return not any(work)


# ---------------------------------------------------------------------------
# Pauli operators and extraspecial 2-group utilities
# ---------------------------------------------------------------------------

def symplectic_inner(lhs: Tuple[BinaryVector, BinaryVector], rhs: Tuple[BinaryVector, BinaryVector]) -> int:
    """Compute the symplectic inner product ``(a | b) · (a' | b')`` modulo 2.

    The value equals ``a · b' + a' · b (mod 2)`` and determines whether the
    corresponding Pauli operators commute (value 0) or anticommute (value 1).
    """
    (a, b), (a_p, b_p) = lhs, rhs
    if len(a) != len(a_p) or len(b) != len(b_p):
        raise ValueError("Symplectic product requires equal lengths.")
    dot1 = sum(x & y for x, y in zip(a, b_p))
    dot2 = sum(x & y for x, y in zip(a_p, b))
    return (dot1 + dot2) & 1


def quadratic_form(pair: Tuple[BinaryVector, BinaryVector]) -> int:
    """Return the quadratic form ``Q(a | b) = Σ a_j b_j (mod 2)``."""
    a, b = pair
    if len(a) != len(b):
        raise ValueError("Quadratic form requires vectors of equal length.")
    return sum(x & y for x, y in zip(a, b)) & 1


@dataclass(frozen=True)
class Pauli:
    r"""Element of the extraspecial 2-group ``E`` in binary symplectic form.

    Parameters
    ----------
    x :
        Binary vector describing where the Pauli has an ``X`` component.
    z :
        Binary vector describing where the Pauli has a ``Z`` component.
    phase :
        A bit encoding the global sign ``(-1)^phase``.  The group operation
        only ever produces real phases, matching the real Clifford group
        analysed in the paper (no ``±i`` factors).
    """

    x: BinaryVector
    z: BinaryVector
    phase: int = 0

    def __post_init__(self) -> None:
        if len(self.x) != len(self.z):
            raise ValueError("Binary vectors x and z must have the same length.")
        object.__setattr__(self, "x", _ensure_binary_vector(self.x))
        object.__setattr__(self, "z", _ensure_binary_vector(self.z))
        object.__setattr__(self, "phase", int(self.phase) & 1)

    @property
    def n(self) -> int:
        """Number of qubits."""
        return len(self.x)

    def support(self) -> Tuple[int, ...]:
        """Indices of non-identity tensor factors."""
        return tuple(
            idx for idx, (xi, zi) in enumerate(zip(self.x, self.z)) if xi or zi
        )

    def weight(self) -> int:
        """Hamming weight, i.e. the number of acted-upon qubits."""
        return len(self.support())

    def quadratic(self) -> int:
        """Return the quadratic form ``Q(\\bar{e})``."""
        return quadratic_form((self.x, self.z))

    def symplectic(self, other: "Pauli") -> int:
        """Return the symplectic inner product with ``other``."""
        return symplectic_inner((self.x, self.z), (other.x, other.z))

    def commutes_with(self, other: "Pauli") -> bool:
        """Check whether two Pauli operators commute."""
        return self.symplectic(other) == 0

    def __mul__(self, other: "Pauli") -> "Pauli":
        """Group multiplication in the extraspecial 2-group."""
        if self.n != other.n:
            raise ValueError("Cannot multiply Paulis acting on different spaces.")
        phase = (self.phase + other.phase + self.symplectic(other)) & 1
        x = tuple((xi ^ xj) for xi, xj in zip(self.x, other.x))
        z = tuple((zi ^ zj) for zi, zj in zip(self.z, other.z))
        return Pauli(x, z, phase)

    def __neg__(self) -> "Pauli":
        """Return the Pauli with opposite sign."""
        return Pauli(self.x, self.z, self.phase ^ 1)

    def __str__(self) -> str:
        symbols = []
        for xi, zi in zip(self.x, self.z):
            if xi and zi:
                symbols.append("Y")
            elif xi:
                symbols.append("X")
            elif zi:
                symbols.append("Z")
            else:
                symbols.append("I")
        prefix = "-" if self.phase else ""
        return prefix + "".join(symbols)

    def to_pair(self) -> Tuple[BinaryVector, BinaryVector]:
        """Return the ``(a | b)`` pair."""
        return self.x, self.z


# ---------------------------------------------------------------------------
# Stabilizer codes
# ---------------------------------------------------------------------------

def is_totally_isotropic(paulis: Iterable[Pauli]) -> bool:
    """Check whether a family of Pauli operators forms a totally isotropic set.

    The test verifies both ``Q(\\bar{e}) = 0`` for each element and mutual
    commutativity (zero symplectic inner product) for every pair.
    """
    pauli_list = list(paulis)
    for pauli in pauli_list:
        if pauli.quadratic() != 0:
            return False
    for idx, p_i in enumerate(pauli_list):
        for p_j in pauli_list[idx + 1 :]:
            if p_i.symplectic(p_j):
                return False
    return True


class StabilizerCode:
    """Finite-dimensional stabilizer code derived from a totally isotropic set.

    Parameters
    ----------
    stabilizers :
        A sequence of commuting :class:`Pauli` operators that generate the
        stabilizer group.  They must have trivial quadratic form to match the
        real Clifford algebra presentation of the paper.
    logicals :
        Optional tuple of Pauli operators representing logical operators.  When
        provided, the constructor validates that they commute with the
        stabilizer and exhibit the expected symplectic structure.
    """

    def __init__(
        self,
        stabilizers: Sequence[Pauli],
        logicals: Optional[Sequence[Pauli]] = None,
    ) -> None:
        if not stabilizers:
            raise ValueError("A stabilizer code needs at least one generator.")
        n = stabilizers[0].n
        for pauli in stabilizers[1:]:
            if pauli.n != n:
                raise ValueError("All stabilizers must act on the same number of qubits.")
        if not is_totally_isotropic(stabilizers):
            raise ValueError("Stabilizers must form a totally isotropic subspace.")

        self._stabilizers: Tuple[Pauli, ...] = tuple(stabilizers)
        self._n = n

        matrix = [pauli.x + pauli.z for pauli in self._stabilizers]
        self._stabilizer_basis = gf2_row_space_basis(matrix)
        self._rank = len(self._stabilizer_basis)

        self._logicals: Tuple[Pauli, ...] = tuple()
        if logicals:
            self._logicals = tuple(logicals)
            for pauli in self._logicals:
                if pauli.n != n:
                    raise ValueError("Logical operators must act on the same space.")
                for stab in self._stabilizers:
                    if pauli.symplectic(stab):
                        raise ValueError("Logical operators must commute with the stabilizer.")
            # Basic sanity check: logicals should not lie in the stabilizer row space
            for pauli in self._logicals:
                if vector_in_row_space(pauli.x + pauli.z, self._stabilizer_basis):
                    raise ValueError("Logical operator is contained in the stabilizer group.")
            # If two logicals anticommute, ensure one of each type is present.
            for idx, p_i in enumerate(self._logicals):
                for p_j in self._logicals[idx + 1 :]:
                    if p_i.symplectic(p_j) not in (0, 1):
                        raise AssertionError("Symplectic product must be a single bit.")

    # -- properties -----------------------------------------------------

    @property
    def n(self) -> int:
        """Number of physical qubits."""
        return self._n

    @property
    def k(self) -> int:
        """Number of logical qubits (``n - rank(H)``)."""
        return self.n - self._rank

    @property
    def stabilizers(self) -> Tuple[Pauli, ...]:
        """Return the stabilizer generators."""
        return self._stabilizers

    @property
    def logicals(self) -> Tuple[Pauli, ...]:
        """Return the logical operators (may be empty)."""
        return self._logicals

    @property
    def parameters(self) -> dict:
        """Dictionary with the canonical ``[[n, k, d]]`` code parameters."""
        distance = self.distance()
        return {"n": self.n, "k": self.k, "d": distance}

    # -- structural checks ---------------------------------------------

    def is_totally_isotropic(self) -> bool:
        """Re-check that the stabilizer generators are totally isotropic."""
        return is_totally_isotropic(self._stabilizers)

    def contains(self, pauli: Pauli) -> bool:
        """Decide whether a Pauli operator belongs to the stabilizer group."""
        if pauli.n != self.n:
            return False
        return vector_in_row_space(pauli.x + pauli.z, self._stabilizer_basis)

    # -- distance estimation -------------------------------------------

    def distance(self, max_weight: Optional[int] = None) -> Optional[int]:
        """Return the minimum weight of a logical operator (code distance).

        The routine performs an exhaustive search up to ``max_weight``.  If the
        maximum weight is not specified, it defaults to the number of qubits.
        The function returns ``None`` when a decisive value cannot be found
        within the requested range.
        """
        limit = self.n if max_weight is None else min(self.n, max_weight)
        if limit <= 0:
            return None

        for weight in range(1, limit + 1):
            for candidate in _pauli_operators_of_weight(self.n, weight):
                if not all(candidate.commutes_with(stab) for stab in self._stabilizers):
                    continue
                if self.contains(candidate):
                    continue
                return weight
        return None

    # -- Stim circuit synthesis ---------------------------------------

    def to_stim_circuit(
        self,
        rounds: int = 1,
        noise_probability: float = 0.0,
        logical_basis: Optional[Sequence[Pauli]] = None,
        backend: str = "auto",
    ) -> Union[str, "stim.Circuit"]:
        """Return a Stim circuit describing repeated stabilizer measurements.

        The circuit that is generated mirrors the stabilizer measurement
        protocol discussed in :mod:`stim` tutorials: during each cycle every
        stabilizer generator is measured using an ``MPP`` instruction.  When
        multiple rounds are requested, ``DETECTOR`` instructions compare the
        outcome of a stabilizer with its value in the previous round, enabling
        detection of time-like error chains.  Optional single-qubit
        depolarising noise can be injected before each round.

        Parameters
        ----------
        rounds:
            Number of stabilizer measurement rounds (``>= 1``).  Each round
            corresponds to a complete quantum error-correction (QEC) cycle.
        noise_probability:
            Probability for a single-qubit depolarising fault applied to each
            data qubit ahead of every round.  Set to zero to omit noise.
        logical_basis:
            Iterable of Pauli operators to be read out as logical observables
            after the syndrome rounds.  When omitted, the code's logical
            operators provided at construction time are used (if any).
        backend:
            Controls the returned representation.  Accepted values are
            ``"auto"`` (default), ``"text"``, and ``"stim"``.  When ``stim`` is
            requested and the optional dependency is unavailable, a
            :class:`ModuleNotFoundError` is raised.  In ``"auto"`` mode the
            method produces a :class:`stim.Circuit` object when Stim is
            available, otherwise it falls back to the textual representation.

        Returns
        -------
        Union[str, stim.Circuit]
            Either a textual Stim circuit description or a compiled
            :class:`stim.Circuit` instance, depending on ``backend`` and the
            availability of the optional dependency.
        """
        if rounds <= 0:
            raise ValueError("Number of rounds must be strictly positive.")
        if noise_probability < 0 or noise_probability >= 1:
            raise ValueError("Noise probability must be in the interval [0, 1).")

        backend = backend.lower()
        if backend not in {"auto", "stim", "text"}:
            raise ValueError("Backend must be one of 'auto', 'stim', or 'text'.")

        if logical_basis is None:
            logicals = self._logicals
        else:
            logicals = tuple(logical_basis)
            for pauli in logicals:
                if pauli.n != self.n:
                    raise ValueError("Logical basis elements must act on the same space.")

        stabilizer_count = len(self._stabilizers)
        if stabilizer_count == 0:
            raise RuntimeError("No stabilizers available to build a circuit.")

        description = _build_stim_text_description(
            n_qubits=self.n,
            stabilizers=self._stabilizers,
            logicals=logicals,
            rounds=rounds,
            noise_probability=noise_probability,
        )

        if backend == "text":
            return description

        try:
            import stim  # type: ignore
        except ModuleNotFoundError:
            if backend == "stim":
                raise
            return description

        return stim.Circuit(description)


# ---------------------------------------------------------------------------
# Constructors mirroring the paper
# ---------------------------------------------------------------------------

PAULI_FROM_LABEL = {
    "I": (0, 0),
    "X": (1, 0),
    "Z": (0, 1),
    "Y": (1, 1),
}


def pauli_from_product(label: str, sign: int = 0) -> Pauli:
    """Create a :class:`Pauli` from a string of tensor factors.

    The string may contain the characters ``I``, ``X``, ``Y``, ``Z`` and is
    case-insensitive.  Whitespace is ignored.  The optional ``sign`` flag
    toggles the global phase ``(-1)^sign``.
    """
    cleaned = label.replace(" ", "").upper()
    x_bits, z_bits = [], []
    for char in cleaned:
        if char not in PAULI_FROM_LABEL:
            raise ValueError(f"Unknown Pauli symbol: {char!r}")
        x_bit, z_bit = PAULI_FROM_LABEL[char]
        x_bits.append(x_bit)
        z_bits.append(z_bit)
    return Pauli(tuple(x_bits), tuple(z_bits), sign)


def pauli_from_binary_pair(pair: str, sign: int = 0) -> Pauli:
    """Create a :class:`Pauli` from a binary pair string of the form ``a|b``."""
    if "|" not in pair:
        raise ValueError("Expected 'a|b' format when parsing binary pair.")
    x_str, z_str = pair.split("|", maxsplit=1)
    x_bits = _parse_binary_string(x_str)
    z_bits = _parse_binary_string(z_str)
    if len(x_bits) != len(z_bits):
        raise ValueError("Binary pair must have equal lengths on both sides.")
    return Pauli(x_bits, z_bits, sign)


def css_code(
    hx: Sequence[BinaryVector],
    hz: Sequence[BinaryVector],
) -> StabilizerCode:
    """Construct a CSS stabilizer code from classical parity-check matrices.

    This implements the standard Calderbank–Shor–Steane (CSS) construction
    described in Example 2 of the paper: we choose classical codes ``C_X`` and
    ``C_Z`` with parity-check matrices ``H_X`` and ``H_Z`` such that
    ``H_X H_Z^T = 0``.  The stabilizer generators are products of ``X`` (from
    rows of ``H_X``) and ``Z`` (from rows of ``H_Z``).
    """
    hx = tuple(_ensure_binary_vector(row) for row in hx)
    hz = tuple(_ensure_binary_vector(row) for row in hz)
    if not hx and not hz:
        raise ValueError("At least one of H_X or H_Z must be non-empty.")

    n = len(hx[0]) if hx else len(hz[0])
    for row in hx:
        if len(row) != n:
            raise ValueError("All rows in H_X must have the same length.")
    for row in hz:
        if len(row) != n:
            raise ValueError("All rows in H_Z must have the same length.")

    # Orthogonality condition H_X H_Z^T = 0
    for row_x in hx:
        for row_z in hz:
            if sum(x & z for x, z in zip(row_x, row_z)) & 1:
                raise ValueError("Parity-check matrices are not mutually orthogonal.")

    stabilizers: List[Pauli] = []
    for row in hx:
        stabilizers.append(Pauli(row, tuple(0 for _ in range(n))))
    for row in hz:
        stabilizers.append(Pauli(tuple(0 for _ in range(n)), row))

    return StabilizerCode(stabilizers)


def five_qubit_code() -> StabilizerCode:
    """Return the cyclic five-qubit code discussed in Example 1.

    The stabilizer is generated by cyclic permutations of ``XZZXI``.  The code
    encodes one qubit, has distance three, and is perfect with respect to a
    single-qubit Pauli error model.
    """
    generators = [
        pauli_from_product("X Z Z X I"),
        pauli_from_product("I X Z Z X"),
        pauli_from_product("X I X Z Z"),
        pauli_from_product("Z X I X Z"),
    ]
    logical_x = pauli_from_product("X X X X X")
    logical_z = pauli_from_product("Z Z Z Z Z")
    return StabilizerCode(generators, logicals=(logical_x, logical_z))


def steane_code() -> StabilizerCode:
    """Return the seven-qubit Steane CSS code from the [7,4,3] Hamming code."""
    # Parity check matrix for the classical [7,4,3] Hamming code.
    hx = (
        (1, 0, 1, 0, 1, 0, 1),
        (0, 1, 1, 0, 0, 1, 1),
        (0, 0, 0, 1, 1, 1, 1),
    )
    hz = hx  # Self-dual CSS construction
    return css_code(hx, hz)


# ---------------------------------------------------------------------------
# Internal helpers
# ---------------------------------------------------------------------------

def _pauli_operators_of_weight(n: int, weight: int) -> Iterable[Pauli]:
    """Generate Paulis on ``n`` qubits with a given weight."""
    if weight < 0 or weight > n:
        return ()
    axes = ("X", "Y", "Z")
    for positions in combinations(range(n), weight):
        for tensor_factors in product(axes, repeat=weight):
            label = ["I"] * n
            for pos, axis in zip(positions, tensor_factors):
                label[pos] = axis
            yield pauli_from_product("".join(label))


def _pauli_to_mpp_string(pauli: Pauli) -> str:
    """Return the Stim ``MPP`` string describing ``pauli``."""
    targets: List[str] = []
    for index, (x_bit, z_bit) in enumerate(zip(pauli.x, pauli.z)):
        if x_bit == 0 and z_bit == 0:
            continue
        if x_bit and z_bit:
            label = "Y"
        elif x_bit:
            label = "X"
        else:
            label = "Z"
        targets.append(f"{label}{index}")
    if not targets:
        raise ValueError("Cannot form an MPP instruction from the identity operator.")
    if pauli.phase:
        targets[0] = "!" + targets[0]
    return "*".join(targets)


def _rec_label(offset: int) -> str:
    """Return a Stim ``rec`` label offset steps into the past."""
    if offset < 1:
        raise ValueError("Offset must be a positive integer.")
    if offset == 1:
        return "rec[-1]"
    return f"rec[-{offset}]"


def _build_stim_text_description(
    *,
    n_qubits: int,
    stabilizers: Sequence[Pauli],
    logicals: Sequence[Pauli],
    rounds: int,
    noise_probability: float,
) -> str:
    """Construct a textual Stim circuit description."""
    lines: List[str] = []
    lines.append(f"# Stim circuit auto-generated for an [[{n_qubits}, ?, ?]] stabilizer code.")
    lines.append("# Data qubits are arranged linearly along the x-axis.")
    for qubit in range(n_qubits):
        lines.append(f"QUBIT_COORDS({qubit}) {float(qubit)} 0")
    lines.append("TICK")

    all_qubits = " ".join(str(q) for q in range(n_qubits))
    stabilizer_count = len(stabilizers)
    previous_offset = stabilizer_count + 1

    for round_index in range(rounds):
        lines.append(f"# --- QEC round {round_index + 1} ---")
        if noise_probability:
            lines.append(f"DEPOLARIZE1({noise_probability}) {all_qubits}")
            lines.append("TICK")
        for stabilizer_index, stabilizer in enumerate(stabilizers):
            mpp_string = _pauli_to_mpp_string(stabilizer)
            lines.append(f"MPP {mpp_string}")
            detector_x = stabilizer_index + 0.5
            detector_y = round_index + 0.5
            if round_index == 0:
                lines.append(f"DETECTOR {_rec_label(1)} {detector_x} {detector_y}")
            else:
                lines.append(
                    f"DETECTOR {_rec_label(1)} {_rec_label(previous_offset)} {detector_x} {detector_y}"
                )
        lines.append("TICK")

    if logicals:
        lines.append("# --- Logical observable readout ---")
    for observable_index, logical in enumerate(logicals):
        mpp_string = _pauli_to_mpp_string(logical)
        lines.append("TICK")
        lines.append(f"MPP {mpp_string}")
        lines.append(f"OBSERVABLE_INCLUDE {_rec_label(1)} {observable_index} {rounds + 1.0}")

    return "\n".join(lines) + "\n"


def _debug_logical_anticommutation(logicals: Sequence[Pauli]) -> List[Tuple[int, int]]:
    """Return indices of logical operator pairs that anticommute (for testing)."""
    issues: List[Tuple[int, int]] = []
    for idx, lhs in enumerate(logicals):
        for jdx, rhs in enumerate(logicals[idx + 1 :], start=idx + 1):
            if lhs.symplectic(rhs) == 1:
                issues.append((idx, jdx))
    return issues


# ---------------------------------------------------------------------------
# Module entry point for manual experimentation
# ---------------------------------------------------------------------------

def _example_usage() -> None:
    """Simple demonstration when the module is executed as a script."""
    print("Five-qubit code parameters:", five_qubit_code().parameters)
    print("Steane code parameters:", steane_code().parameters)


if __name__ == "__main__":  # pragma: no cover - manual inspection aide
    _example_usage()