diff --git a/src/qldpc/abstract/rings.py b/src/qldpc/abstract/rings.py
index 0869eb1f1..8a56f4cd0 100644
--- a/src/qldpc/abstract/rings.py
+++ b/src/qldpc/abstract/rings.py
@@ -1082,10 +1082,9 @@ def _get_block_howell_form(matrix: galois.FieldArray, *, right: bool = False) ->
     if right and size > 1:
         matrix = matrix.transpose(0, 1, 3, 2)
 
-    # row-reduce as an expanded 2-D matrix and remove all-zero rows
+    # row-reduce as an expanded 2-D matrix, keeping a basis of the row space (no all-zero rows)
     shape = (num_block_rows * size, num_block_cols * size)
-    matrix = matrix.transpose(0, 2, 1, 3).reshape(shape).view(field).row_reduce()
-    matrix = matrix[qldpc.math.first_nonzero_cols(matrix) < matrix.shape[1]].view(field)
+    matrix = matrix.transpose(0, 2, 1, 3).reshape(shape).view(field).row_space()
 
     if size > 1:
         # insert zero rows to shift pivots down so that they always lie on the diagonal of a block
diff --git a/src/qldpc/circuits/encoding.py b/src/qldpc/circuits/encoding.py
index 6dd1c885c..8d64dfcd7 100644
--- a/src/qldpc/circuits/encoding.py
+++ b/src/qldpc/circuits/encoding.py
@@ -226,7 +226,7 @@ def get_state_stabilizers(
 
     # identify stabilizers that are supported entirely on the data qubits of the code
     stabilizers = []
-    for row in matrix.row_reduce():
+    for row in matrix.row_space():
         state_xs = row[cols_state_x]
         state_zs = row[cols_state_z]
         other_xs = row[cols_other_x]
diff --git a/src/qldpc/codes/common.py b/src/qldpc/codes/common.py
index 45a8f4e73..7fdef879a 100644
--- a/src/qldpc/codes/common.py
+++ b/src/qldpc/codes/common.py
@@ -34,6 +34,7 @@
 import scipy.linalg
 import scipy.sparse
 import stim
+from typing_extensions import Self
 
 from qldpc import abstract, decoders, external, math
 from qldpc.math import IntegerArray
@@ -203,9 +204,10 @@ def graph_to_matrix(graph: nx.DiGraph) -> galois.FieldArray:
     def dimension(self) -> int:
         """The number of logical (qu)dits encoded by this code."""
 
-    def forget_distance(self) -> None:
+    def forget_distance(self) -> Self:
         """Forget the known distance of this code."""
         self._distance = None
+        return self
 
 
 ################################################################################
@@ -265,10 +267,9 @@ def canonicalized(self) -> ClassicalCode:
         """The same code with its parity matrix in reduced row echelon form."""
         if self._is_canonicalized:  # pragma: no cover
             return self
-        matrix_rref = self.matrix.row_reduce()
-        matrix_rref = matrix_rref[np.any(matrix_rref, axis=1), :]
-        code = ClassicalCode(matrix_rref, self.field)
-        code._dimension = len(self) - len(matrix_rref)
+        matrix = self.matrix.row_space()
+        code = ClassicalCode(matrix, self.field)
+        code._dimension = len(self) - len(matrix)
         code._distance = self._distance
         code._is_canonicalized = True
         return code
@@ -343,7 +344,7 @@ def generator(self) -> galois.FieldArray:
             self._generator = self.matrix.null_space()
         return self._generator
 
-    def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int]]) -> None:
+    def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int]]) -> Self:
         """Set the generator matrix of this code."""
         generator = np.asanyarray(generator).view(self.field)
         if np.any(self.matrix @ generator.T):
@@ -357,6 +358,7 @@ def set_generator(self, generator: npt.NDArray[np.int_] | Sequence[Sequence[int]
                 f" {required_rank})"
             )
         self._generator = generator
+        return self
 
     def iter_words(self, skip_zero: bool = False) -> Iterator[galois.FieldArray]:
         """Iterate over the code words of this code."""
@@ -891,11 +893,10 @@ def canonicalized(self) -> QuditCode:
         """The same code with its parity matrix in reduced row echelon form."""
         if self._is_canonicalized:  # pragma: no cover
             return self
-        matrix_rref = self.matrix.row_reduce()
-        matrix_rref = matrix_rref[np.any(matrix_rref, axis=1), :]
-        code = QuditCode(matrix_rref, self.field, is_subsystem_code=self._is_subsystem_code)
+        matrix = self.matrix.row_space()
+        code = QuditCode(matrix, self.field, is_subsystem_code=self._is_subsystem_code)
         if not self._is_subsystem_code:
-            code._dimension = len(code) - len(matrix_rref)
+            code._dimension = len(code) - len(matrix)
         code._distance = self._distance
         code._stabilizer_ops = self._stabilizer_ops
         code._gauge_ops = self._gauge_ops
@@ -950,7 +951,10 @@ def maybe_to_css(self) -> QuditCode:
         zs = np.any(matrix_z, axis=1)
         if np.any(xs & zs):
             return self
-        return CSSCode(matrix_x[xs], matrix_z[zs], is_subsystem_code=self._is_subsystem_code)
+        code = CSSCode(matrix_x[xs], matrix_z[zs], is_subsystem_code=self._is_subsystem_code)
+        code._dimension = self._dimension
+        code._distance = self._distance
+        return code
 
     def to_css(self) -> CSSCode:
         """Try to convert this QuditCode into a CSSCode.  Throw an error if we fail."""
@@ -1235,7 +1239,7 @@ def get_standard_form_data(
         parity check matrix (modulo elementary row operations and array reshaping) is
         matrix[:, :, np.argsort(qudit_locs)].  As a sanity check, the following test should pass:
             matrix_2d = matrix[:, :, np.argsort(qudit_locs)].reshape(-1, 2 * len(self))
-            assert np.array_equal(matrix_2d.row_reduce(), self.canonicalized.matrix)
+            assert np.array_equal(matrix_2d.row_space(), self.canonicalized.matrix)
 
         Finally, this method also returns slices (index sets) for all row and column sectors, which
         enables selecting blocks of the parity check matrix with, say, matrix[rows_sx, 1, cols_lz].
@@ -1384,7 +1388,7 @@ def set_logical_ops(
         logical_ops: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the logical operators of this code to the provided logical operators."""
         logical_ops = np.asanyarray(logical_ops).view(self.field)
         if not skip_validation:
@@ -1400,13 +1404,14 @@ def set_logical_ops(
                 raise ValueError("An incorrect number of logical operators was provided")
         self._logical_ops = logical_ops
         self._dimension = len(logical_ops) // 2
+        return self
 
     def set_logical_ops_x(
         self,
         logicals_ops_x: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the X-type logical operators of this code.
 
         Determine suitable Z-type logical operators automatically.  This choice is unique mod
@@ -1431,7 +1436,7 @@ def set_logical_ops_x(
         old_logicals_z = self.get_logical_ops(Pauli.Z, symplectic=True)
         basis_change = np.linalg.inv(math.symplectic_conjugate(old_logicals_z) @ logicals_ops_x.T)
         new_logicals_z = basis_change @ old_logicals_z
-        self.set_logical_ops(
+        return self.set_logical_ops(
             np.vstack([logicals_ops_x, new_logicals_z]), skip_validation=skip_validation
         )
 
@@ -1440,7 +1445,7 @@ def set_logical_ops_z(
         logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the Z-type logical operators of this code.
 
         Determine suitable X-type logical operators automatically.  This choice is unique mod
@@ -1465,7 +1470,7 @@ def set_logical_ops_z(
         old_logicals_x = self.get_logical_ops(Pauli.X, symplectic=True)
         basis_change = np.linalg.inv(old_logicals_x @ math.symplectic_conjugate(logicals_ops_z).T)
         new_logicals_x = basis_change @ old_logicals_x
-        self.set_logical_ops(
+        return self.set_logical_ops(
             np.vstack([new_logicals_x, logicals_ops_z]), skip_validation=skip_validation
         )
 
@@ -1747,6 +1752,8 @@ def transform_ops(ops: galois.FieldArray) -> galois.FieldArray:
             code._stabilizer_ops = transform_ops(self.get_stabilizer_ops())
         if self._gauge_ops is not None:
             code._gauge_ops = transform_ops(self.get_gauge_ops())
+        code._dimension = self._dimension
+        code._distance = self._distance
         return code
 
     def conjugate(self) -> QuditCode:
@@ -1795,6 +1802,7 @@ def transform_ops(ops: galois.FieldArray) -> galois.FieldArray:
         if self._gauge_ops is not None:
             code._gauge_ops = transform_ops(self.get_gauge_ops())
 
+        code._dimension = self._dimension
         return code
 
     @staticmethod
@@ -2198,18 +2206,49 @@ def is_self_dual(self) -> bool:
     def is_swel(self) -> bool:
         """Is this code self-dual with equivalent logicals (SWEL)?
 
-        Specifically, SWEL codes are self-dual codes for which Hadamard-transforming every logical Z
-        operator recovers a dual logical X operator.
+        A SWEL code is a self-dual CSS code for which there exists a basis of logical operators
+        (Lx, Lz) with Lx = Lz; that is, Hadamard-transforming any single-qudit logical Pauli X (in
+        some logical Pauli basis) recovers the associated logical Pauli Z.
 
-        Whether a code is SWEL depends on the choice of logical operator basis.
+        In the case of stabilizer (non-subsystem) codes over a field with characteristic 2 (that
+        is, over GF(2**m), which includes qubits with m = 1), a self-dual CSS code is SWEL iff it
+        has odd block length (Corollary 1 of arXiv:2503.19790), or at least one row of Lx (in any
+        basis) has nonzero self-overlap (Theorem 1 of arXiv:2503.19790).
+
+        This method first checks special cases covered by arXiv:2503.19790, and otherwise directly
+        checks the existence of a logical operator basis (Lx, Lz) with Lx = Lz.
         """
-        return bool(
-            self.is_self_dual
-            and np.array_equal(
-                (logs_z := self.get_logical_ops(Pauli.Z)) @ logs_z.T,
-                np.eye(self.dimension, dtype=int),
-            )
+        if not self.is_self_dual:
+            return False
+        ops_x = self.get_logical_ops(Pauli.X)
+        if not self.is_subsystem_code and self.field.characteristic == 2:
+            return len(self) % 2 == 1 or any(op @ op for op in ops_x)
+        return math.get_orthonormal_basis(ops_x, promise_full_rank=True) is not None
+
+    def get_swel_logical_ops(self) -> galois.FieldArray:
+        """Find a self-dual basis of logical operators for this code: (Lx, Lz) = (L, L).  Return L.
+
+        Raise a ValueError if no such basis exists (see QuditCode.is_swel).
+        """
+        supports = (
+            math.get_orthonormal_basis(self.get_logical_ops(Pauli.X), promise_full_rank=True)
+            if self.is_self_dual
+            else None
         )
+        if supports is None:
+            raise ValueError(
+                "This code has no self-dual logical operator basis; it must be self-dual, and its"
+                " logical operators must admit an orthonormal basis"
+            )
+        return supports
+
+    def set_swel_logical_ops(self) -> Self:
+        """Set the logical operators of this code to those of a self-dual logical operator basis.
+
+        Raise a ValueError if this code is not SWEL (see QuditCode.is_swel).
+        """
+        supports = self.get_swel_logical_ops()
+        return self.set_logical_ops_xz(supports, supports)
 
     @functools.cached_property
     def canonicalized(self) -> CSSCode:
@@ -2490,17 +2529,17 @@ def set_logical_ops_xz(
         logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the logical operators of this code to the provided logical operators."""
         logical_ops = scipy.linalg.block_diag(logicals_ops_x, logicals_ops_z)
-        self.set_logical_ops(logical_ops, skip_validation=skip_validation)
+        return self.set_logical_ops(logical_ops, skip_validation=skip_validation)
 
     def set_logical_ops_x(
         self,
         logicals_ops_x: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the X-type logical operators of this code.
 
         Determine suitable Z-type logical operators automatically.  This choice is unique mod
@@ -2519,14 +2558,16 @@ def set_logical_ops_x(
         logicals_ops_x = np.asanyarray(logicals_ops_x).view(self.field)
         old_logicals_z = self.get_logical_ops(Pauli.Z)
         new_logicals_z = np.linalg.inv(old_logicals_z @ logicals_ops_x.T) @ old_logicals_z
-        self.set_logical_ops_xz(logicals_ops_x, new_logicals_z, skip_validation=skip_validation)
+        return self.set_logical_ops_xz(
+            logicals_ops_x, new_logicals_z, skip_validation=skip_validation
+        )
 
     def set_logical_ops_z(
         self,
         logicals_ops_z: npt.NDArray[np.int_] | Sequence[Sequence[int]],
         *,
         skip_validation: bool = False,
-    ) -> None:
+    ) -> Self:
         """Set the Z-type logical operators of this code.
 
         Determine suitable X-type logical operators automatically.  This choice is unique mod
@@ -2545,7 +2586,9 @@ def set_logical_ops_z(
         logicals_ops_z = np.asanyarray(logicals_ops_z).view(self.field)
         old_logicals_x = self.get_logical_ops(Pauli.X)
         new_logicals_x = np.linalg.inv(old_logicals_x @ logicals_ops_z.T) @ old_logicals_x
-        self.set_logical_ops_xz(new_logicals_x, logicals_ops_z, skip_validation=skip_validation)
+        return self.set_logical_ops_xz(
+            new_logicals_x, logicals_ops_z, skip_validation=skip_validation
+        )
 
     def get_stabilizer_ops(
         self,
@@ -2595,7 +2638,7 @@ def get_gauge_ops(
         return gauge_ops.reshape(-1, 2, len(self))[:, pauli, :].view(self.field)
 
     def dual(self) -> CSSCode:
-        """Dual to this code, which swaps the roles of logical and gauge operators.
+        """The dual of this code, which swaps the roles of logical and gauge operators.
 
         See help(qldpc.codes.QuditCode.dual) for an explanation.
         """
@@ -2863,11 +2906,12 @@ def get_distance_bound_with_decoder(
 
         return min_bound
 
-    def forget_distance(self) -> None:
+    def forget_distance(self) -> Self:
         """Forget the known distance of this code."""
         self._distance_x = self._distance_z = self._distance = None
+        return self
 
-    def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs: Any) -> None:
+    def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs: Any) -> Self:
         """Reduce the weight of a logical operator.
 
         A minimal-weight logical operator is found by enforcing that it has a trivial syndrome, and
@@ -2899,8 +2943,9 @@ def reduce_logical_op(self, pauli: PauliXZ, logical_index: int, **decoder_kwargs
         self._logical_ops.shape = (2, self.dimension, 2, len(self))
         self._logical_ops[pauli, logical_index, pauli, :] = candidate_logical_op
         self._logical_ops.shape = (2 * self.dimension, 2 * len(self))
+        return self
 
-    def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any) -> None:
+    def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any) -> Self:
         """Reduce the weight of all logical operators."""
         assert pauli is None or pauli in PAULIS_XZ
         if pauli is None:
@@ -2909,6 +2954,7 @@ def reduce_logical_ops(self, pauli: PauliXZ | None = None, **decoder_kwargs: Any
         else:
             for logical_index in range(self.dimension):
                 self.reduce_logical_op(pauli, logical_index, **decoder_kwargs)
+        return self
 
     def conjugated(self, qudits: slice | Sequence[int]) -> QuditCode:
         """Apply local Fourier transforms, swapping X-type and Z-type operators.
@@ -2933,6 +2979,10 @@ def conjugate(self) -> CSSCode:
             code._gauge_ops = scipy.linalg.block_diag(
                 self.get_gauge_ops(Pauli.Z), self.get_gauge_ops(Pauli.X)
             ).view(self.field)
+        code._dimension = self._dimension
+        code._distance = self._distance
+        code._distance_x = self._distance_z
+        code._distance_z = self._distance_x
         return code
 
     def deformed(
diff --git a/src/qldpc/codes/common_test.py b/src/qldpc/codes/common_test.py
index 6e9e598c8..07ad813db 100644
--- a/src/qldpc/codes/common_test.py
+++ b/src/qldpc/codes/common_test.py
@@ -599,9 +599,27 @@ def test_css_code(pytestconfig: pytest.Config) -> None:
     assert nx.utils.graphs_equal(subgraphs[0], code.get_graph(Pauli.X))
     assert nx.utils.graphs_equal(subgraphs[1], code.get_graph(Pauli.Z))
 
-    # self-dual codes with equivalent logicals (SWEL)
-    assert codes.SteaneCode().is_swel
-    assert not codes.SurfaceCode(3).is_swel
+
+def test_swel_codes() -> None:
+    """Identify and construct SWEL logical operator bases (see CSSCode.is_swel)."""
+    steane_code = codes.SteaneCode()
+    surface_code = codes.SurfaceCode(3)
+    even_code = codes.CSSCode([[1, 1, 1, 1]], [[1, 1, 1, 1]])
+    assert steane_code.is_swel
+    assert not surface_code.is_swel  # not self-dual
+    assert not even_code.is_swel  # self-dual but not SWEL
+
+    # a self-dual code over GF(3) reaches the general orthonormal-basis test (odd characteristic)
+    qutrit_matrix = [[0, 0, 1, 1, 1], [0, 1, 0, 1, 2]]
+    assert codes.CSSCode(qutrit_matrix, qutrit_matrix, field=3).is_swel
+
+    # we automatically find a self-dual logical operator basis, if it exists
+    code = steane_code.set_swel_logical_ops()
+    assert np.array_equal(code.get_logical_ops(Pauli.X), code.get_logical_ops(Pauli.Z))
+    with pytest.raises(ValueError, match="no self-dual logical operator basis"):
+        surface_code.get_swel_logical_ops()
+    with pytest.raises(ValueError, match="no self-dual logical operator basis"):
+        even_code.get_swel_logical_ops()
 
 
 def test_css_ops(pytestconfig: pytest.Config) -> None:
diff --git a/src/qldpc/math.py b/src/qldpc/math.py
index 44563b54f..2e5f75cb8 100644
--- a/src/qldpc/math.py
+++ b/src/qldpc/math.py
@@ -35,6 +35,24 @@
 DenseIntegerArrayType = TypeVar("DenseIntegerArrayType", galois.FieldArray, npt.NDArray[np.int_])
 
 
+####################################################################################################
+# general math
+
+
+@functools.cache
+def log_choose(n: int, k: int) -> float:
+    """Natural logarithm of (n choose k) = n! / ( k! * (n-k)! )."""
+    return (
+        scipy.special.gammaln(n + 1)
+        - scipy.special.gammaln(k + 1)
+        - scipy.special.gammaln(n - k + 1)
+    )
+
+
+####################################################################################################
+# manipulating Pauli strings and their symplectic vector representations
+
+
 def op_to_string(op: npt.NDArray[np.int_]) -> stim.PauliString:
     """Convert an integer array that represents a Pauli string into a stim.PauliString.
 
@@ -80,14 +98,21 @@ def symplectic_weight(vectors: npt.NDArray[np.int_]) -> int:
     return np.count_nonzero(vectors_x | vectors_z, axis=-1).reshape(vectors.shape[:-1])
 
 
+####################################################################################################
+# matrix helper functions
+
+
 def first_nonzero_cols(
     matrix: npt.NDArray[np.generic] | Sequence[npt.NDArray[np.generic]],
 ) -> npt.NDArray[np.int_]:
     """Get the first nonzero column for every row in a matrix.
 
     If all columns are zero in a particular row, the column value is set to the number of columns.
-    If the array has more than two dimensions, return for each "row" r, the first "column" c for
-    which array[r, c] is not all zero.
+
+    If the input ``matrix`` has more than two dimensions, return an ``output`` vector of length
+    ``matrix.shape[0]``, where ``output[r]`` is the first "column" ``c`` for which
+    ``np.any(matrix[r, c])`` is nonzero; (or set ``output[r] = matrix.shape[1]`` if
+    ``not np.any(matrix[r])``).
     """
     _matrix = np.atleast_2d(np.asarray(matrix))
     if _matrix.size == 0:
@@ -99,19 +124,6 @@ def first_nonzero_cols(
     return first_nonzero_col_index.astype(np.int_, copy=False)
 
 
-def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois.FieldArray:
-    """Construct a dual basis, for which dual_basis @ basis.T = identity_matrix."""
-    if validate and (
-        basis.shape[0] > basis.shape[1] or np.linalg.matrix_rank(basis) != basis.shape[0]
-    ):
-        raise ValueError("A dual basis can only be found for wide matrices of full rank")
-    pivot_cols = first_nonzero_cols(basis.row_reduce())
-    linearly_independent_cols = basis[:, pivot_cols].view(type(basis))
-    dual_basis = np.zeros(basis.shape, dtype=int).view(type(basis))
-    dual_basis[:, pivot_cols] = np.linalg.inv(linearly_independent_cols).T
-    return dual_basis.view(type(basis))
-
-
 def block_matrix(
     blocks: Sequence[Sequence[npt.NDArray[np.generic] | int | object]],
 ) -> npt.NDArray[np.int_]:
@@ -158,11 +170,191 @@ def block_matrix(
     return matrix
 
 
-@functools.cache
-def log_choose(n: int, k: int) -> float:
-    """Natural logarithm of (n choose k) = n! / ( k! * (n-k)! )."""
-    return (
-        scipy.special.gammaln(n + 1)
-        - scipy.special.gammaln(k + 1)
-        - scipy.special.gammaln(n - k + 1)
+####################################################################################################
+# basis builders
+
+
+def get_dual_basis(basis: galois.FieldArray, *, validate: bool = True) -> galois.FieldArray:
+    """Construct a dual basis, for which dual_basis @ basis.T = identity_matrix."""
+    if validate and (
+        basis.shape[0] > basis.shape[1] or np.linalg.matrix_rank(basis) != basis.shape[0]
+    ):
+        raise ValueError("A dual basis can only be found for wide matrices of full rank")
+    pivot_cols = first_nonzero_cols(basis.row_reduce())
+    linearly_independent_cols = basis[:, pivot_cols].view(type(basis))
+    dual_basis = np.zeros(basis.shape, dtype=int).view(type(basis))
+    dual_basis[:, pivot_cols] = np.linalg.inv(linearly_independent_cols).T
+    return dual_basis.view(type(basis))
+
+
+def get_orthonormal_basis(
+    matrix: galois.FieldArray, *, promise_full_rank: bool = False
+) -> galois.FieldArray | None:
+    """An orthonormal basis for the row space of a matrix over a finite field.
+
+    Given a matrix whose rows span a subspace V of GF(q)^n, return a matrix L whose rows are a
+    basis for V with L @ L.T = identity -- that is, a basis of V whose vectors each have
+    self-overlap v @ v = 1 and are pairwise orthogonal under the dot product.  If V has no such
+    basis, return None.
+
+    The rows of the matrix may be linearly dependent; they are first reduced to a basis of V.  Pass
+    promise_full_rank=True to skip this reduction when the rows are already independent.
+
+    An orthonormal basis for V exists if and only if (a) V is nondegenerate, meaning no nonzero
+    vector of V is orthogonal to all of V, and (b), according to the field's characteristic:
+    - over a field of characteristic 2, some vector of V has nonzero self-overlap v @ v;
+    - over a field of odd characteristic, an even number of the vectors in any mutually orthogonal
+      basis of V have a self-overlap with no square root in GF(q).
+
+    The construction is a variant of Gram-Schmidt orthogonalization; the characteristic-2 case
+    follows Algorithm 1 and Lemma 2 of arXiv:2503.19790.
+    """
+    field = type(matrix)
+    dimension = matrix.shape[1]
+    if not promise_full_rank:
+        matrix = matrix.row_space()  # reduce to a basis, discarding linearly dependent rows
+    words = list(matrix)
+    units = (
+        _orthonormalize_char_2(words)
+        if field.characteristic == 2
+        else _orthonormalize_odd(words, field)
     )
+    if units is None:
+        return None
+    return field(units) if units else field.Zeros((0, dimension))
+
+
+def _orthonormalize_char_2(words: list[galois.FieldArray]) -> list[galois.FieldArray] | None:
+    """Try to orthonormalize linearly independent vectors over a field of characteristic 2.
+
+    Reduce the row space to mutually orthogonal "unit" vectors u with u @ u = 1 and "hyperbolic
+    pairs" (b, c) with b @ b = c @ c = 0 and b @ c = 1 (Algorithm 1 of arXiv:2503.19790).  Every
+    element of a characteristic-2 field is a square, so any vector with nonzero self-overlap can be
+    rescaled to a unit vector.  Lemma 2 of arXiv:2503.19790 then rewrites one unit vector and one
+    hyperbolic pair into three unit vectors, eliminating every hyperbolic pair.
+    """
+    units: list[galois.FieldArray] = []  # vectors u with u @ u = 1
+    pairs: list[tuple[galois.FieldArray, galois.FieldArray]] = []  # (b, c) with b @ c = 1
+    while words:
+        index = next((ii for ii, word in enumerate(words) if word @ word), None)
+        if index is not None:
+            # extract a unit vector and orthogonalize other words against it
+            pivot = words.pop(index)
+            unit = pivot / _sqrt(pivot @ pivot)
+            words = [word - (word @ unit) * unit for word in words]
+            units.append(unit)
+        else:
+            # extract a hyperbolic pair and orthogonalize other words against both of its vectors
+            first = words.pop(0)
+            index = next((ii for ii, word in enumerate(words) if first @ word), None)
+            if index is None:
+                return None  # "first" is orthogonal to all of V, so no orthonormal basis exists
+            partner = words.pop(index)
+            partner = partner / (first @ partner)  # rescale so that first @ partner == 1
+            words = [word - (word @ partner) * first - (word @ first) * partner for word in words]
+            pairs.append((first, partner))
+
+    # an orthonormal basis needs at least one unit vector; pure hyperbolic pairs have none
+    if not units and pairs:
+        return None
+
+    # Lemma 2: rewrite a unit vector a and a hyperbolic pair (b, c) as three unit vectors
+    while pairs:
+        vec_b, vec_c = pairs.pop()
+        vec_a = units.pop()
+        units += [vec_a + vec_b + vec_c, vec_a + vec_b, vec_a + vec_c]
+    return units
+
+
+def _orthonormalize_odd(
+    words: list[galois.FieldArray], field: type[galois.FieldArray]
+) -> list[galois.FieldArray] | None:
+    """Try to orthonormalize linearly independent vectors over a field of odd characteristic.
+
+    First, find an orthogonal basis in which each vector has nonzero self-overlap (with Gram-Schmidt
+    orthogonalization): at each step, pick a vector with nonzero self-overlap v @ v and subtract
+    its projection from the others.  Out of the remaining null vectors (with zero self-overlap), if
+    a pair (u, v) have nonzero overlap u @ v, then their sum has nonzero self-overlap (in a field
+    with odd characteristic), (u + v) @ (u + v) = 2 * u @ v, so replace u <- u + v, and
+    orthogonalize remaining vectors against u.  If no remaining null vectors have nonzero overlap,
+    no orthonormal basis exists.
+
+    A vector with nonzero self-overlap can be rescaled to a unit vector (self-overlap 1) exactly
+    when its self-overlap has a square root in the field, by taking v -> v / sqrt(v @ v).  Vectors
+    with non-square self-overlaps cannot be rescaled on their own, but a pair of them may be
+    combined into two unit vectors, allowing an orthonormal basis when the number of such vectors
+    is even.
+    """
+    # diagonalize: build an orthogonal basis of vectors with nonzero self-overlap
+    diagonal: list[tuple[galois.FieldArray, galois.FieldArray]] = []
+    while words:
+        index = next((ii for ii, word in enumerate(words) if word @ word), None)
+        if index is None:
+            # every self-overlap vanishes; sum an overlapping pair to get a nonzero self-overlap
+            index = _combine_null_vectors(words)
+            if index is None:
+                return None  # every overlap vanishes on the remaining space; no basis exists
+        # extract a vector with nonzero overlap, and otrhogonalize the remaining vectors against it
+        pivot = words.pop(index)
+        overlap = pivot @ pivot
+        words = [word - (word @ pivot) / overlap * pivot for word in words]
+        diagonal.append((pivot, overlap))
+
+    # rescale square-norm vectors to unit vectors; collect the non-square-norm vectors
+    units: list[galois.FieldArray] = []
+    non_squares: list[tuple[galois.FieldArray, galois.FieldArray]] = []
+    for pivot, overlap in diagonal:
+        if field.is_square(overlap):
+            units.append(pivot / _sqrt(overlap))
+        else:
+            non_squares.append((pivot, overlap))
+
+    # non-square self-overlaps pair off in twos, so an odd count means no orthonormal basis
+    if len(non_squares) % 2:
+        return None
+
+    # pair non-square-norm vectors into unit vectors, using a fixed non-square element epsilon
+    if non_squares:
+        epsilon = field.primitive_element  # a generator of GF(q)* is always a non-square
+        # solve alpha**2 + beta**2 = 1 / epsilon, so that alpha u + beta w has self-overlap 1
+        alpha, beta = _sum_of_two_squares(field, field(1) / epsilon)
+        for (u_vec, u_norm), (w_vec, w_norm) in zip(non_squares[::2], non_squares[1::2]):
+            # rescale u_vec and w_vec to self-overlap epsilon, then combine into two unit vectors
+            u_vec = u_vec * _sqrt(epsilon / u_norm)
+            w_vec = w_vec * _sqrt(epsilon / w_norm)
+            units += [alpha * u_vec + beta * w_vec, -beta * u_vec + alpha * w_vec]
+    return units
+
+
+def _combine_null_vectors(words: list[galois.FieldArray]) -> int | None:
+    """Add one word to another to obtain a nonzero self-overlap; return its index, or None.
+
+    Assumes every word has zero self-overlap.  In odd characteristic, if some pair has nonzero
+    overlap then their sum has nonzero self-overlap; if no pair does, every overlap is zero.
+    """
+    for ii in range(len(words)):
+        for jj in range(ii + 1, len(words)):
+            if words[ii] @ words[jj]:
+                words[ii] = words[ii] + words[jj]
+                return ii
+    return None
+
+
+def _sum_of_two_squares(
+    field: type[galois.FieldArray], target: galois.FieldArray
+) -> tuple[galois.FieldArray, galois.FieldArray]:
+    """Field elements (a, b) with a**2 + b**2 == target; these exist in odd characteristic."""
+    for alpha in field.elements:  # roughly half of all "alpha" succeed, so this exits quickly
+        remainder = target - alpha * alpha
+        if field.is_square(remainder):
+            return alpha, _sqrt(remainder)
+    raise AssertionError("no solution to a**2 + b**2 = target")  # pragma: no cover
+
+
+def _sqrt(value: galois.FieldArray) -> galois.FieldArray:
+    """A square root of a (0-dimensional) field element that is known to be a square.
+
+    Delegates to np.sqrt, which computes finite-field square roots but must be given an array
+    rather than a 0-dimensional scalar (which it rejects over some extension fields).
+    """
+    return np.sqrt(np.atleast_1d(value))[0]
diff --git a/src/qldpc/math_test.py b/src/qldpc/math_test.py
index bbff624b1..2d54fc433 100644
--- a/src/qldpc/math_test.py
+++ b/src/qldpc/math_test.py
@@ -60,18 +60,53 @@ def test_nonzero_cols() -> None:
     assert np.array_equal(qldpc.math.first_nonzero_cols(tensor), [0])
 
 
-def test_dual_basis(pytestconfig: pytest.Config) -> None:
+def test_dual_basis() -> None:
     """Construct dual bases."""
-    np.random.seed(pytestconfig.getoption("randomly_seed"))
-
     field = galois.GF(2)
-    basis = field.Random((4, 5)).row_reduce()
-    basis = basis[qldpc.math.first_nonzero_cols(basis) < basis.shape[1]]
+    basis = field([[1, 1, 0, 0, 1], [0, 1, 1, 0, 0], [1, 0, 0, 1, 1]])
     dual_basis = qldpc.math.get_dual_basis(basis)
     assert np.array_equal(dual_basis @ basis.T, field.Identity(len(basis)))
 
     with pytest.raises(ValueError, match="wide matrices of full rank"):
-        qldpc.math.get_dual_basis(field.Random((2, 1)))
+        qldpc.math.get_dual_basis(field([[1], [0]]))
+
+
+def test_orthonormal_basis() -> None:
+    """Orthonormal bases over finite fields (arXiv:2503.19790, Algorithm 1 and Lemma 2)."""
+    # subspaces that admit an orthonormal basis: check that L @ L.T is the identity
+    with_basis = [
+        galois.GF(2).Zeros((0, 4)),  # the empty subspace
+        galois.GF(2)([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 0], [0, 0, 0, 0, 1, 1]]),  # Lemma 2
+        galois.GF(4)([[2, 0]]),  # characteristic 2, a rescaled unit vector
+        galois.GF(3)([[1, 0], [0, 1]]),  # odd characteristic, unit norms
+        galois.GF(7)([[1, 1]]),  # odd characteristic, a rescaled square norm
+        galois.GF(7)([[1, 2], [1, 3]]),  # odd characteristic, two non-square norms paired off
+    ]
+    # these inputs are full-rank bases crafted to exercise each branch of the orthonormalization
+    for vectors in with_basis:
+        field = type(vectors)
+        basis = qldpc.math.get_orthonormal_basis(vectors, promise_full_rank=True)
+        assert basis is not None and np.array_equal(basis @ basis.T, field.Identity(len(basis)))
+
+    # subspaces with no orthonormal basis: get_orthonormal_basis returns None
+    without_basis = [
+        # characteristic 2, an alternating form (every vector has even weight)
+        galois.GF(2)(
+            [[1, 1, 0, 0, 0, 0], [1, 0, 1, 0, 0, 0], [0, 1, 1, 1, 1, 0], [0, 0, 0, 1, 0, 1]]
+        ),
+        galois.GF(2)([[1, 1, 0, 0], [0, 0, 1, 1]]),  # a degenerate form
+        galois.GF(7)([[1, 2]]),  # odd characteristic, a non-square discriminant
+        galois.GF(7)([[1, 2, 3]]),  # odd characteristic, a degenerate (isotropic) form
+        galois.GF(7)([[1, 2, 3], [1, 2, 4]]),  # odd characteristic, isotropic then non-square
+    ]
+    for vectors in without_basis:
+        assert qldpc.math.get_orthonormal_basis(vectors, promise_full_rank=True) is None
+
+    # by default (promise_full_rank=False) the rows are first reduced to a basis of their row space
+    vectors = galois.GF(3)([[1, 0], [0, 1], [1, 1]])  # three dependent rows spanning all of GF(3)^2
+    basis = qldpc.math.get_orthonormal_basis(vectors)
+    assert basis is not None and basis.shape == (2, 2)
+    assert np.array_equal(basis @ basis.T, galois.GF(3).Identity(2))
 
 
 def test_block_matrix() -> None: