diff --git a/CHANGELOG.md b/CHANGELOG.md
index ce0b0cb7f..6d362c6a0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,6 +5,10 @@
 
 # News
 
+## v0.9.16 - dev
+
+- Enhancements to `GF(2)` Linear Algebra: unexported, experimental `gf2_row_echelon_with_pivots!`, `gf2_nullspace`, `gf2_rowspace_basis`.
+
 ## v0.9.15 - 2024-12-22
 
 - `pftrajectories` now supports fast multiqubit measurements with `PauliMeasurement` in addition to the already supported single qubit measurements `sMX/Z/Y` and workarounds like `naive_syndrome_circuit`.
diff --git a/src/QuantumClifford.jl b/src/QuantumClifford.jl
index 6b90af7a8..25f84a218 100644
--- a/src/QuantumClifford.jl
+++ b/src/QuantumClifford.jl
@@ -1168,6 +1168,69 @@ function gf2_H_to_G(H)
     G[:,invperm(sindx)]
 end
 
+"""Performs in-place Gaussian elimination on a binary matrix and returns
+its *row echelon form*,*rank*, the *transformation matrix*, and the *pivot
+columns*. The transformation matrix that converts the original matrix into
+the row echelon form. The `full` parameter controls the extent of elimination:
+if `true`, only rows below the pivot are affected; if `false`, both above and
+below the pivot are eliminated."""
+function gf2_row_echelon_with_pivots!(M::AbstractMatrix{Int}; full=false)
+    r, c = size(M)
+    N = Matrix{Int}(LinearAlgebra.I, r, r)
+    p = 1
+    pivots = Int[]
+    for col in 1:c
+        @inbounds for row in p:r
+            if M[row, col] == 1
+                if row != p
+                    M[[row, p], :] .= M[[p, row], :]
+                    N[[row, p], :] .= N[[p, row], :]
+                end
+                break
+            end
+        end
+        if M[p, col] == 1
+            if !full
+                elim_range = p+1:r
+            else
+                elim_range = 1:r
+            end
+            @simd for j in elim_range
+                @inbounds if j != p && M[j, col] == 1
+                    M[j, :] .= (M[j, :] .+ M[p, :]) .% 2
+                    N[j, :] .= (N[j, :] .+ N[p, :]) .% 2
+                end
+            end
+            p += 1
+            push!(pivots, col)
+        end
+        if p > r
+            break
+        end
+    end
+    rank = p - 1
+    return M, rank, N, pivots
+end
+
+"""The nullspace of a binary matrix."""
+function gf2_nullspace(H::AbstractMatrix{Int})
+    m = size(H',1)
+    _, matrix_rank, transformation_matrix, _ = gf2_row_echelon_with_pivots!(copy(H)')
+    if m == matrix_rank
+        # By the rank-nullity theorem, if rank(M) = m, then nullity(M) = 0
+        return zeros(Bool, 1, m)
+    end
+    # Extract the nullspace from the transformation matrix
+    return transformation_matrix[matrix_rank+1:end, :]
+end
+
+"""The basis for the row space of the binary matrix."""
+function gf2_rowspace_basis(H::AbstractMatrix{Int})
+    pivots = gf2_row_echelon_with_pivots!(copy(H)')[4]
+    # Extract the rows corresponding to the pivot columns
+    H[pivots,:]
+end
+
 ##############################
 # Error classes
 ##############################
diff --git a/test/test_row_echelon_with_pivots.jl b/test/test_row_echelon_with_pivots.jl
new file mode 100644
index 000000000..98d72429e
--- /dev/null
+++ b/test/test_row_echelon_with_pivots.jl
@@ -0,0 +1,165 @@
+@testitem "Row echelon with pivots" begin
+    using Random
+    using Nemo
+    using Nemo: echelon_form, matrix, GF
+    using QuantumClifford
+    using QuantumClifford: gf2_row_echelon_with_pivots!, gf2_nullspace, gf2_rowspace_basis
+    test_sizes = [1,2,10,63,64,65,127,128,129]
+
+    @testset "GF(2) row echelon form with transformation matrix, pivots etc." begin
+        for n in test_sizes
+            for rep in 1:10
+                gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
+                for (i, mat) in enumerate(gf2_matrices)
+                    naive_echelon_form, _, transformation, _ = gf2_row_echelon_with_pivots!(Matrix{Int}(mat), full=true) # in-place
+                    # Check the correctness of the transformation matrix
+                    @test (transformation*mat) .%2 == naive_echelon_form
+                    # Check the correctness of Gaussian elimination
+                    @test naive_echelon_form == gf2_gausselim!(mat)
+                    # Consistency check with Nemo.jl's echelon_form
+                    nemo_mat =  matrix(GF(2), Matrix{Int}(mat))
+                    @test echelon_form(nemo_mat) == matrix(GF(2), naive_echelon_form)
+                end
+            end
+        end
+    end
+
+    function is_in_nullspace(A, x)
+        # Ensure x is the correct orientation
+        if size(x, 1) != size(A, 2)
+            x = transpose(x)
+        end
+        # Perform modulo 2 arithmetic: A * x must be zero mod 2
+        if size(x, 2) == 1  # x is a single column vector
+            result = A * x
+            return all(result .% 2 .== 0)  # Check if A * x = 0 mod 2
+        else  # x is a matrix, check each column vector
+            for i in 1:size(x, 2)
+                result = A * x[:, i] # Multiply A with the i-th column of x
+                if !all(result .% 2 .== 0) # Check if A * column = 0 mod 2
+                    return false
+                end
+            end
+            return true # All columns are in the null space mod 2
+        end
+    end
+
+    @testset "GF(2) nullspace of the binary matrix" begin
+        for n in test_sizes
+            for rep in 1:10
+                gf2_matrices = [rand(Bool, size, size) for size in test_sizes]
+                for (i, matrix) in enumerate(gf2_matrices)
+                    imat = Matrix{Int}(matrix)
+                    ns = gf2_nullspace(imat)
+                    @test is_in_nullspace(imat, ns)
+                end
+            end
+        end
+    end
+
+    @testset "Consistency check with ldpc" begin
+        # sanity checks for comparison to https://github.com/quantumgizmos/ldpc
+        # results compared with 'from ldpc.mod2 import nullspace, row_basis, row_echelon'
+        # Consistency check 1
+        H = [1 1 1; 1 1 1; 0 1 0]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 1 1; 0 1 0; 0 0 0]
+        @test rank == 2
+        @test transformation == [1 0 0; 0 0 1; 1 1 0]
+        @test pivots == [1, 2] # in python, it's [0, 1] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [1 0 1]
+        @test gf2_rowspace_basis(copy(H)) == [1 1 1; 0 1 0]
+        # Consistency check 2
+        H = [0 0 0 1 1 1 1;
+             0 1 1 0 0 1 1;
+             1 0 1 0 1 0 1]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 0 1 0 1 0 1;
+                               0 1 1 0 0 1 1;
+                               0 0 0 1 1 1 1]
+        @test rank == 3
+        @test transformation == [0 0 1;
+                                 0 1 0;
+                                 1 0 0]
+        @test pivots == [1, 2, 4] # in python, it's [0, 1, 3] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
+                                         0 1 1 1 1 0 0;
+                                         0 1 0 1 0 1 0;
+                                         0 0 1 1 0 0 1]
+        @test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
+                                              0 1 1 0 0 1 1;
+                                              1 0 1 0 1 0 1]
+        # Consistency check 3
+        H = [1 1 0; 0 1 1; 1 0 1]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 1 0;
+                               0 1 1;
+                               0 0 0]
+        @test rank == 2
+        @test transformation == [1 0 0;
+                                 0 1 0;
+                                 1 1 1]
+        @test pivots == [1,2 ] # in python, it's [0, 1] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [1 1 1]
+        @test gf2_rowspace_basis(copy(H)) == [1 1 0;
+                                              0 1 1]
+        # Consistency check 4
+        H = [1 1 0; 0 1 0; 0 0 1]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 1 0;
+                               0 1 0;
+                               0 0 1]
+        @test rank == 3
+        @test transformation == [1 0 0;
+                                 0 1 0;
+                                 0 0 1]
+        @test pivots == [1, 2, 3]  # in python, it's [0, 1, 2] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [0 0 0]
+        @test gf2_rowspace_basis(copy(H)) == [1 1 0;
+                                              0 1 0;
+                                              0 0 1]
+        # Consistency check 5
+        H = [1 1 0; 0 1 0; 0 0 1; 0 1 1]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 1 0;
+                               0 1 0;
+                               0 0 1;
+                               0 0 0]
+        @test rank == 3
+        @test transformation == [1 0 0 0;
+                                 0 1 0 0;
+                                 0 0 1 0;
+                                 0 1 1 1]
+        @test pivots == [1, 2, 3] # in python, it's [0, 1, 2] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [0 0 0]
+        @test gf2_rowspace_basis(copy(H)) == [1 1 0;
+                                              0 1 0;
+                                              0 0 1]
+        # Consistency check 6
+        H = [0 0 0 1 1 1 1;
+             0 1 1 0 0 1 1;
+             1 0 1 0 1 0 1]
+        echelon_form, rank, transformation, pivots = gf2_row_echelon_with_pivots!(copy(H)) # in-place
+        @test echelon_form == [1 0 1 0 1 0 1;
+                               0 1 1 0 0 1 1;
+                               0 0 0 1 1 1 1]
+        @test rank == 3
+        @test transformation == [0 0 1;
+                                 0 1 0;
+                                 1 0 0]
+        @test pivots == [1, 2, 4]  # in python, it's [0, 1, 3] due to zero-indexing
+        @test mod.((transformation*copy(H)), 2) == echelon_form
+        @test gf2_nullspace(copy(H)) == [1 1 1 0 0 0 0;
+                                         0 1 1 1 1 0 0;
+                                         0 1 0 1 0 1 0;
+                                         0 0 1 1 0 0 1]
+        @test gf2_rowspace_basis(copy(H)) == [0 0 0 1 1 1 1;
+                                              0 1 1 0 0 1 1;
+                                              1 0 1 0 1 0 1]
+    end
+end