BitMatrix.h

Go to the documentation of this file.

#pragma once
// SPDX-FileCopyrightText: 2025 Nessan Fitzmaurice <nzznfitz+gh@icloud.com>
// SPDX-License-Identifier: MIT

 
#include <gf2/BitPolynomial.h>
#include <gf2/BitVector.h>
#include <gf2/RNG.h>
 
#include <string>
#include <optional>
#include <regex>
 
namespace gf2 {
 
// Forward declarations to avoid recursive inclusion issues
// clang-format off
template<Unsigned Word> class BitGauss;
template<Unsigned Word> class BitLU;
// clang-format on

template<Unsigned Word = usize>
class BitMatrix {
private:
    std::vector<BitVector<Word>> m_rows;
 
public:
    // The row type is a `BitVector` ...
    using row_type = BitVector<Word>;

    using word_type = Word;


    constexpr BitMatrix() : m_rows{} {}

    constexpr BitMatrix(usize n) : m_rows{} {
        if (n > 0) resize(n, n);
    }

    constexpr BitMatrix(usize m, usize n) : m_rows{} {
        if (m > 0 && n > 0) resize(m, n);
    }

    constexpr BitMatrix(std::vector<row_type> const& rows) : m_rows{rows} {
        gf2_assert(check_rows(m_rows), "Not all rows have the same size!");
    }

    constexpr BitMatrix(std::vector<BitVector<Word>>&& rows) : m_rows{std::move(rows)} {
        gf2_assert(check_rows(m_rows), "Not all rows have the same size!");
    }


    static constexpr BitMatrix zeros(usize m, usize n) { return BitMatrix{m, n}; }

    static constexpr BitMatrix zeros(usize m) { return BitMatrix{m, m}; }

    static constexpr BitMatrix ones(usize m, usize n) {
        auto rows = std::vector{m, BitVector<Word>::ones(n)};
        return BitMatrix{std::move(rows)};
    }

    static constexpr BitMatrix ones(usize m) { return BitMatrix::ones(m, m); }

    static constexpr BitMatrix alternating(usize m, usize n) {
        auto rows = std::vector{m, BitVector<Word>::alternating(n)};
        // Flip every other row.
        for (auto i = 1uz; i < m; i += 2) rows[i].flip_all();
        return BitMatrix{std::move(rows)};
    }

    static constexpr BitMatrix alternating(usize m) { return BitMatrix::alternating(m, m); }

    template<BitStore Lhs, BitStore Rhs>
        requires std::same_as<typename Lhs::word_type, Word> && std::same_as<typename Rhs::word_type, Word>
    static constexpr BitMatrix from_outer_product(Lhs const& u, Rhs const& v) {
        BitMatrix result;
        for (auto i = 0uz; i < u.size(); ++i) {
            row_type row = u.get(i) ? row_type(v) : row_type(v.size());
            result.m_rows.push_back(std::move(row));
        }
        return result;
    }

    template<BitStore Lhs, BitStore Rhs>
        requires std::same_as<typename Lhs::word_type, Word> && std::same_as<typename Rhs::word_type, Word>
    static constexpr BitMatrix from_outer_sum(Lhs const& u, Rhs const& v) {
        BitMatrix result;
        for (auto i = 0uz; i < u.size(); ++i) {
            result.m_rows.push_back(v);
            if (u.get(i)) result.m_rows.back().flip_all();
        }
        return result;
    }

    static constexpr BitMatrix from(usize m, usize n, std::invocable<usize, usize> auto f) {
        BitMatrix result{m, n};
        for (auto i = 0uz; i < m; ++i)
            for (auto j = 0uz; j < n; ++j)
                if (f(i, j)) result.set(i, j, true);
        return result;
    }


    static BitMatrix random(usize m, usize n, double p = 0.5, std::uint64_t seed = 0) {
        // Keep a single static RNG per thread for all calls to this method, seeded with entropy on the first call.
        thread_local RNG rng;
 
        // Edge case handling ...
        if (p < 0) return BitMatrix::zeros(m, n);
 
        // Scale p by 2^64 to remove floating point arithmetic from the main loop below.
        // If we determine p rounds to 1 then we can just set all elements to 1 and return early.
        p = p * 0x1p64 + 0.5;
        if (p >= 0x1p64) return BitMatrix::ones(m, n);
 
        // p does not round to 1 so we use a 64-bit URNG and check each draw against the 64-bit scaled p.
        auto scaled_p = static_cast<std::uint64_t>(p);
 
        // If a seed was provided, set the RNG's seed to it. Otherwise, we carry on from where we left off.
        std::uint64_t old_seed = rng.seed();
        if (seed != 0) rng.set_seed(seed);
 
        auto result = BitMatrix::zeros(m, n);
        for (auto i = 0uz; i < m; ++i) {
            for (auto j = 0uz; j < n; ++j) {
                if (rng.u64() < scaled_p) result.set(i, j, true);
            }
        }
 
        // Restore the old seed if necessary.
        if (seed != 0) rng.set_seed(old_seed);
 
        return result;
    }

    static BitMatrix seeded_random(usize m, usize n, std::uint64_t seed) { return random(m, n, 0.5, seed); }

    static BitMatrix seeded_random(usize m, std::uint64_t seed) { return random(m, m, 0.5, seed); }

    static BitMatrix biased_random(usize m, usize n, double p) { return random(m, n, p, 0); }

    static BitMatrix biased_random(usize m, double p) { return random(m, m, p, 0); }


    static constexpr BitMatrix zero(usize m) { return BitMatrix{m, m}; }

    static constexpr BitMatrix identity(usize m) {
        BitMatrix result{m, m};
        for (auto i = 0uz; i < m; ++i) result.set(i, i, true);
        return result;
    }

    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    static constexpr BitMatrix companion(Store const& top_row) {
        // Edge case:
        if (top_row.size() == 0) return BitMatrix{};
        auto result = BitMatrix::zero(top_row.size());
        result.m_rows[0].copy(top_row);
        result.set_sub_diagonal(1);
        return result;
    }

    static constexpr BitMatrix left_shift(usize n, usize p) {
        auto result = BitMatrix::zeros(n, n);
        result.set_super_diagonal(p);
        return result;
    }

    static BitMatrix right_shift(usize n, usize p) {
        auto result = BitMatrix::zeros(n, n);
        result.set_sub_diagonal(p);
        return result;
    }

    static BitMatrix left_rotation(usize n, usize p) {
        auto result = BitMatrix::zeros(n, n);
        for (auto i = 0uz; i < n; ++i) {
            auto j = (i + n - p) % n;
            result.set(i, j, true);
        }
        return result;
    }

    static BitMatrix right_rotation(usize n, usize p) {
        auto result = BitMatrix::zeros(n, n);
        for (auto i = 0uz; i < n; ++i) {
            auto j = (i + p) % n;
            result.set(i, j, true);
        }
        return result;
    }


    //
    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    static std::optional<BitMatrix> from_row_store(Store const& v, usize r) {
        // Edge case?
        if (v.size() == 0) return BitMatrix{};
 
        // Error handling ...
        if (r == 0 || v.size() % r != 0) return std::nullopt;
 
        // Number of columns (for sure this is an even division)
        usize c = v.size() / r;
 
        // Create the bit-matrix and copy the rows.
        BitMatrix result{r, c};
        for (auto i = 0uz; i < r; ++i) {
            auto begin = i * c;
            auto end = begin + c;
            result.m_rows[i].copy(v.span(begin, end));
        }
        return result;
    }

    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    static std::optional<BitMatrix> from_col_store(Store const& v, usize c) {
        // Edge case?
        if (v.size() == 0) return BitMatrix{};
 
        // Error handling ...
        if (c == 0 || v.size() % c != 0) return std::nullopt;
 
        // Number of columns (for sure this is an even division)
        usize r = v.size() / c;
 
        // Create the bit-matrix and copy the rows.
        BitMatrix result{r, c};
        usize     iv = 0;
        for (auto j = 0uz; j < c; ++j) {
            for (auto i = 0uz; i < r; ++i) {
                if (v.get(iv)) result.m_rows[i].set(j);
                iv++;
            }
        }
        return result;
    }


    static std::optional<BitMatrix> from_string(std::string_view s) {
        // Edge case
        if (s.empty()) return BitMatrix{};
 
        // We split the string into tokens using the standard regex library.
        std::string                src(s);
        auto                       delims = std::regex(R"([\s|,|;]+)");
        std::sregex_token_iterator iter{src.cbegin(), src.cend(), delims, -1};
        std::vector<std::string>   tokens{iter, {}};
 
        // Zap any empty tokens & check there is something to do
        tokens.erase(std::remove_if(tokens.begin(), tokens.end(), [](std::string_view x) { return x.empty(); }),
                     tokens.end());
        if (tokens.empty()) return BitMatrix{};
 
        // We hope to fill a BitMatrix.
        BitMatrix result;
 
        // Iterate through the possible rows
        usize n_rows = tokens.size();
        usize n_cols = 0;
        for (auto i = 0uz; i < n_rows; ++i) {
            // Attempt to parse the current token as a row for the bit-matrix.
            auto r = row_type::from_string(tokens[i]);
 
            // Parse failure?
            if (!r) return std::nullopt;
 
            // We've read a potentially valid BitMatrix row.
            if (i == 0) {
                // First row sets the number of columns
                n_cols = r->size();
                result.resize(n_rows, n_cols);
            } else {
                // Subsequent rows better have the same number of elements!
                if (r->size() != n_cols) return std::nullopt;
            }
            result.m_rows[i] = *r;
        }
        return result;
    }


    constexpr usize rows() const { return m_rows.size(); }

    constexpr usize cols() const { return m_rows.size() > 0 ? m_rows[0].size() : 0; }

    constexpr usize size() const { return rows() * cols(); }

    constexpr bool is_empty() const { return rows() == 0; }


    constexpr bool is_square() const { return rows() != 0 && rows() == cols(); }

    constexpr bool is_zero() const { return is_square() && none(); }

    constexpr bool is_identity() const {
        if (!is_square()) return false;
        for (auto i = 0uz; i < rows(); ++i) {
            row_type r = m_rows[i];
            r.flip(i);
            if (r.any()) return false;
        }
        return true;
    }

    constexpr bool is_symmetric() const {
        if (!is_square()) return false;
        for (auto i = 0uz; i < rows(); ++i) {
            for (auto j = 0uz; j < cols(); ++j) {
                if (get(i, j) != get(j, i)) return false;
            }
        }
        return true;
    }


    constexpr usize count_ones() const {
        usize count = 0;
        for (const auto& row : m_rows) count += row.count_ones();
        return count;
    }

    constexpr usize count_zeros() const { return size() - count_ones(); }

    constexpr usize count_ones_on_diagonal() const {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        usize result = 0;
        for (auto i = 0uz; i < rows(); ++i)
            if (get(i, i)) ++result;
        return result;
    }

    constexpr bool trace() const { return count_ones_on_diagonal() % 2 == 1; }


    constexpr bool any() const {
        for (const auto& row : m_rows)
            if (row.any()) return true;
        return false;
    }

    constexpr bool all() const {
        for (const auto& row : m_rows)
            if (!row.all()) return false;
        return true;
    }

    constexpr bool none() const {
        for (const auto& row : m_rows)
            if (!row.none()) return false;
        return true;
    }


    constexpr bool get(usize r, usize c) const {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        gf2_debug_assert(c < cols(), "Column index {} out of bounds [0,{})", c, cols());
        return m_rows[r].get(c);
    }

    constexpr bool operator()(usize r, usize c) const {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        gf2_debug_assert(c < cols(), "Column index {} out of bounds [0,{})", c, cols());
        return m_rows[r].get(c);
    }

    constexpr void set(usize r, usize c, bool val = true) {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        gf2_debug_assert(c < cols(), "Column index {} out of bounds [0,{})", c, cols());
        m_rows[r].set(c, val);
    }

    constexpr BitRef<row_type> operator()(usize r, usize c) {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        gf2_debug_assert(c < cols(), "Column index {} out of bounds [0,{})", c, cols());
        return m_rows[r][c];
    }

    constexpr void flip(usize r, usize c) {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        gf2_debug_assert(c < cols(), "Column index {} out of bounds [0,{})", c, cols());
        m_rows[r].flip(c);
    }


    constexpr const row_type& row(usize r) const {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        return m_rows[r];
    }

    constexpr row_type& row(usize r) {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        return m_rows[r];
    }

    constexpr const row_type& operator[](usize r) const {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        return m_rows[r];
    }

    constexpr row_type& operator[](usize r) {
        gf2_debug_assert(r < rows(), "Row index {} out of bounds [0,{})", r, rows());
        return m_rows[r];
    }


    constexpr BitVector<Word> col(usize c) const {
        gf2_debug_assert(c < cols(), "Column {} out of bounds [0, {})", c, cols());
        auto result = BitVector<Word>::zeros(rows());
        for (auto r = 0uz; r < rows(); ++r) {
            if (get(r, c)) { result.set(r); }
        }
        return result;
    }


    constexpr void set_all(bool value = true) {
        for (auto& row : m_rows) row.set_all(value);
    }

    constexpr void flip_all() {
        for (auto& row : m_rows) row.flip_all();
    }


    constexpr void set_diagonal(bool val = true) {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows(); ++i) { set(i, i, val); }
    }

    constexpr void flip_diagonal() {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows(); ++i) { flip(i, i); }
    }

    constexpr void set_super_diagonal(usize d, bool val = true) {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows() - d; ++i) { set(i, i + d, val); }
    }

    constexpr void flip_super_diagonal(usize d) {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows() - d; ++i) { flip(i, i + d); }
    }

    constexpr void set_sub_diagonal(usize d, bool val = true) {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows() - d; ++i) { set(i + d, i, val); }
    }

    constexpr void flip_sub_diagonal(usize d) {
        gf2_debug_assert(is_square(), "BitMatrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows() - d; ++i) { flip(i + d, i); }
    }


    constexpr void resize(usize r, usize c) {
        // Edge case ...
        if (rows() == r && cols() == c) return;
 
        // Resizes to zero in either dimension is taken as a zap the lot instruction
        if (r == 0 || c == 0) r = c = 0;
 
        // Rows, then columns
        m_rows.resize(r);
        for (auto& row : m_rows) row.resize(c);
    }

    constexpr void clear() { resize(0, 0); }

    constexpr void make_square(usize n) { resize(n, n); }


    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    constexpr BitMatrix& append_row(Store const& row) {
        gf2_assert_eq(row.size(), cols(), "Row has {} elements but bit-matrix has {} columns!", row.size(), cols());
        m_rows.push_back(row);
        return *this;
    }

    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    constexpr BitMatrix& append_row(Store&& row) {
        gf2_assert_eq(row.size(), cols(), "Row has {} elements but bit-matrix has {} columns!", row.size(), cols());
        m_rows.push_back(std::move(row));
        return *this;
    }

    constexpr BitMatrix& append_rows(BitMatrix<Word> const& src) {
        gf2_assert_eq(src.cols(), cols(), "Source has {} columns but bit-matrix has {} columns!", src.cols(), cols());
        m_rows.insert(m_rows.end(), src.m_rows.begin(), src.m_rows.end());
        return *this;
    }

    constexpr BitMatrix& append_rows(BitMatrix<Word>&& src) {
        gf2_assert_eq(src.cols(), cols(), "Source has {} columns but bit-matrix has {} columns!", src.cols(), cols());
        m_rows.insert(m_rows.end(), std::make_move_iterator(src.m_rows.begin()),
                      std::make_move_iterator(src.m_rows.end()));
        return *this;
    }

    template<BitStore Store>
        requires std::same_as<typename Store::word_type, Word>
    constexpr BitMatrix& append_col(Store const& col) {
        gf2_assert_eq(col.size(), rows(), "Column has {} elements but bit-matrix has {} rows!", col.size(), rows());
        for (auto i = 0uz; i < rows(); ++i) m_rows[i].push(col.get(i));
        return *this;
    }

    constexpr BitMatrix& append_cols(BitMatrix<Word> const& src) {
        gf2_assert_eq(src.rows(), rows(), "Source has {} rows but bit-matrix has {} rows!", src.rows(), rows());
        for (auto i = 0uz; i < rows(); ++i) m_rows[i].append(src.m_rows[i]);
        return *this;
    }


    constexpr std::optional<BitVector<Word>> remove_row() {
        if (rows() == 0) return std::nullopt;
        auto row = std::move(m_rows.back());
        m_rows.pop_back();
        return row;
    }

    constexpr std::optional<BitMatrix<Word>> remove_rows(usize k) {
        if (rows() < k) return std::nullopt;
        auto begin = m_rows.end() - static_cast<std::ptrdiff_t>(k);
        auto end = m_rows.end();
        auto result = BitMatrix<Word>(std::vector<row_type>(begin, end));
        m_rows.erase(begin, end);
        return result;
    }

    constexpr std::optional<BitVector<Word>> remove_col() {
        if (cols() == 0) return std::nullopt;
        auto result = col(cols() - 1);
        for (auto i = 0uz; i < rows(); ++i) m_rows[i].pop();
        return result;
    }


    constexpr BitMatrix sub_matrix(usize r_start, usize r_end, usize c_start, usize c_end) const {
 
        // Check that the row range is valid.
        gf2_assert(r_start <= r_end, "Invalid row range");
        gf2_assert(r_end <= rows(), "Row range extends beyond the end of the bit-matrix");
 
        // Check that the column range is valid.
        gf2_assert(c_start <= c_end, "Invalid column range");
        gf2_assert(c_end <= cols(), "Column range extends beyond the right edge of the bit-matrix");
 
        // Get the number of rows and columns in the sub-matrix.
        auto r = r_end - r_start;
        auto c = c_end - c_start;
 
        // Create the sub-matrix.
        BitMatrix result{r, c};
        for (auto i = 0uz; i < r; ++i) result.m_rows[i].copy(m_rows[i + r_start].span(c_start, c_end));
 
        return result;
    }

    constexpr void replace_sub_matrix(usize top, usize left, BitMatrix<Word> const& src) {
        auto r = src.rows();
        auto c = src.cols();
        gf2_assert(top + r <= rows(), "Too many rows for the replacement sub-matrix to fit");
        gf2_assert(left + c <= cols(), "Too many columns for the replacement sub-matrix to fit");
        for (auto i = 0uz; i < r; ++i) m_rows[top + i].span(left, left + c).copy(src.m_rows[i]);
    }


    constexpr BitMatrix lower() const {
        // Edge case:
        if (is_empty()) return BitMatrix{};
 
        // Start with a copy of the bit-matrix.
        BitMatrix result = *this;
 
        // Set the upper triangular part to zero.
        auto nc = cols();
        for (auto i = 0uz; i < rows(); ++i) {
            auto first = i + 1;
            if (first < nc) result.m_rows[i].span(first, nc).set_all(false);
        }
        return result;
    }

    constexpr BitMatrix upper() const {
        // Edge case:
        if (is_empty()) return BitMatrix{};
 
        // Start with a copy of the bit-matrix.
        auto result = *this;
 
        // Set the upper triangular part to zero.
        auto c = cols();
        for (auto i = 0uz; i < rows(); ++i) {
            auto len = std::min(i, c);
            if (len > 0) result.m_rows[i].span(0, len).set_all(false);
        }
        return result;
    }

    constexpr BitMatrix strictly_lower() const {
        auto result = lower();
        result.set_diagonal(false);
        return result;
    }

    constexpr BitMatrix strictly_upper() const {
        auto result = upper();
        result.set_diagonal(false);
        return result;
    }

    constexpr BitMatrix unit_lower() const {
        auto result = lower();
        result.set_diagonal(true);
        return result;
    }

    constexpr BitMatrix unit_upper() const {
        auto result = upper();
        result.set_diagonal(true);
        return result;
    }


    constexpr void swap_rows(usize i, usize j) {
        gf2_debug_assert(i < rows(), "Row index {} out of bounds [0,{})", i, rows());
        gf2_debug_assert(j < rows(), "Row index {} out of bounds [0,{})", j, rows());
        std::swap(m_rows[i], m_rows[j]);
    }

    constexpr void swap_cols(usize i, usize j) {
        gf2_debug_assert(i < cols(), "Column index {} out of bounds [0,{})", i, cols());
        gf2_debug_assert(j < cols(), "Column index {} out of bounds [0,{})", j, cols());
        for (auto k = 0uz; k < rows(); ++k) m_rows[k].swap(i, j);
    }

    constexpr void add_identity() {
        gf2_debug_assert(is_square(), "This bit-matrix is {} x {} but it should be square!", rows(), cols());
        for (auto i = 0uz; i < rows(); ++i) flip(i, i);
    }


    constexpr BitMatrix transposed() const {
        auto      r = rows();
        auto      c = cols();
        BitMatrix result{c, r};
        for (auto i = 0uz; i < r; ++i) {
            for (auto j = 0uz; j < c; ++j) {
                if (get(i, j)) { result.set(j, i); }
            }
        }
        return result;
    }

    constexpr void transpose() {
        gf2_assert(is_square(), "`transpose()` requires a square matrix");
        for (auto i = 0uz; i < rows(); ++i) {
            for (auto j = 0uz; j < i; ++j) {
                if (get(i, j) != get(j, i)) {
                    flip(i, j);
                    flip(j, i);
                }
            }
        }
    }


    constexpr BitMatrix to_the(usize n, bool n_is_log2 = false) const {
        gf2_assert(is_square(), "Bit-matrix is {} x {} but it should be square!", rows(), cols());
 
        // Perhaps we just need lots of square steps?  Note that 2^0 = 1 so M^(2^0) = M.
        if (n_is_log2) {
            auto result = *this;
            for (auto i = 0uz; i < n; ++i) result = dot(result, result);
            return result;
        }
 
        // Otherwise we need square & multiply steps but we first handle the edge case:
        if (n == 0) return BitMatrix::identity(rows());
 
        // OK n != 0: Note that if e.g. n = 0b00010111 then std::bit_floor(n) = 0b00010000.
        usize n_bit = std::bit_floor(n);
 
        // Start with a copy of this bit-matrix which takes care of the most significant binary digit in n.
        auto result = *this;
        n_bit >>= 1;
 
        // More to go?
        while (n_bit) {
            // Always square and then multiply if necessary (i.e. if current bit in n is set).
            result = dot(result, result);
            if (n & n_bit) result = dot(result, *this);
 
            // On to the next bit position in n.
            n_bit >>= 1;
        }
        return result;
    }


    BitVector<Word> to_echelon_form() {
        gf2_assert(!is_empty(), "Bit-matrix must not be empty");
 
        // We return a bit-vector that shows which columns have a pivot -- start by assuming none.
        auto has_pivot = BitVector<Word>::zeros(cols());
 
        // The current row of the echelon form we are working on.
        auto r = 0uz;
        auto nr = rows();
 
        // Iterate over each column.
        for (auto j = 0uz; j < cols(); ++j) {
            // Find a non-zero entry in this column below the diagonal (a "pivot").
            auto p = r;
            while (p < nr && !get(p, j)) p++;
 
            // Did we find a pivot in this column?
            if (p < nr) {
                // Mark this column as having a pivot.
                has_pivot.set(j);
 
                // If necessary, swap the current row with the row that has the pivot.
                if (p != r) swap_rows(p, r);
 
                // Below the working row make sure column j is zero by elimination if necessary.
                auto row_r = m_rows[r];
                for (auto i = r + 1; i < nr; ++i)
                    if (get(i, j)) m_rows[i] ^= row_r;
 
                // Move to the next row. If we've reached the end of the matrix, we're done.
                r += 1;
                if (r == nr) break;
            }
        }
 
        // Return the bit-vector that shows which columns have a pivot.
        return has_pivot;
    }

    BitVector<Word> to_reduced_echelon_form() {
        // Start with the echelon form.
        auto has_pivot = to_echelon_form();
 
        // Iterate over each row from the bottom up - Gauss Jordan elimination.
        for (auto r = rows(); r--;) {
            // Find the first set bit in the current row if there is one.
            if (auto p = m_rows[r].first_set()) {
                // Clear out everything in column p *above* row r (already cleared out below the pivot).
                auto row_r = m_rows[r];
                for (auto i = 0uz; i < r; ++i)
                    if (get(i, *p)) m_rows[i] ^= row_r;
            }
        }
 
        // Return the bit-vector that shows which columns have a pivot.
        return has_pivot;
    }


    std::optional<BitMatrix> inverse() const {
        // The bit-matrix must be square & non-empty.
        if (is_empty() || !is_square()) return std::nullopt;
 
        // Create a copy of the bit-matrix & augment it with the identity matrix on the right.
        auto matrix = *this;
        auto nr = rows();
        auto nc = cols();
        matrix.append_cols(BitMatrix::identity(nr));
 
        // Transform the augmented matrix to reduced row-echelon form.
        matrix.to_reduced_echelon_form();
 
        // If all went well the left half is the identity matrix and the right half is the inverse.
        if (matrix.sub_matrix(0, nr, 0, nc).is_identity()) return matrix.sub_matrix(0, nr, nc, nc * 2);
 
        return std::nullopt;
    }

    static constexpr double probability_invertible(usize n) {
        // A 0x0 matrix is not well defined throw an error!
        if (n == 0) throw std::invalid_argument("Matrix should not be 0x0--likely an error somewhere upstream!");
 
        // Formula is p(n) = \prod_{k = 1}^{n} (1 - 2^{-k}) which runs out of juice once n hits any size at all!
        usize n_prod = std::min(n, static_cast<size_t>(std::numeric_limits<double>::digits));
 
        // Compute the product over the range that matters
        double product = 1;
        double pow2 = 1;
        while (n_prod--) {
            pow2 *= 0.5;
            product *= 1 - pow2;
        }
        return product;
    }

    static constexpr double probability_singular(usize n) { return 1 - probability_invertible(n); }


    BitLU<Word> LU() const { return BitLU<Word>{*this}; }


    template<BitStore Rhs>
        requires std::same_as<typename Rhs::word_type, Word>
    auto solver_for(Rhs const& b) const {
        return BitGauss<Word>{*this, b};
    }

    template<BitStore Rhs>
        requires std::same_as<typename Rhs::word_type, Word>
    auto x_for(Rhs const& b) const {
        return solver_for(b)();
    }


    BitPolynomial<Word> characteristic_polynomial() const {
        gf2_assert(is_square(), "Bit-matrix must be square not {} x {}", rows(), cols());
        return frobenius_matrix_characteristic_polynomial(frobenius_form());
    }

    static BitPolynomial<Word>
    frobenius_matrix_characteristic_polynomial(std::vector<BitVector<Word>> const& top_rows) {
        auto n_companions = top_rows.size();
        if (n_companions == 0) return BitPolynomial<Word>::zero();
 
        // Compute the product of the characteristic polynomials of the companion matrices.
        auto result = companion_matrix_characteristic_polynomial(top_rows[0]);
        for (auto i = 1uz; i < n_companions; ++i) { result *= companion_matrix_characteristic_polynomial(top_rows[i]); }
        return result;
    }

    static BitPolynomial<Word> companion_matrix_characteristic_polynomial(BitVector<Word> const& top_row) {
        auto n = top_row.size();
 
        // The characteristic polynomial is degree n with n + 1 coefficients (leading coefficient is 1).
        auto coeffs = BitVector<Word>::ones(n + 1);
 
        // The lower order coefficients are the top row of the companion matrix in reverse order.
        for (auto j = 0uz; j < n; ++j) { coeffs.set(n - j - 1, top_row.get(j)); }
        return BitPolynomial<Word>{std::move(coeffs)};
    }

    std::vector<BitVector<Word>> frobenius_form() const {
        // The bit-matrix must be square.
        gf2_assert(is_square(), "Bit-matrix must be square not {} x {}", rows(), cols());
 
        // Space for the top rows of the companion matrices which we will return.
        auto nr = rows();
        auto top_rows = std::vector<BitVector<Word>>{};
        top_rows.reserve(nr);
 
        // Make a working copy of the bit-matrix to work through using Danilevsky's algorithm.
        auto copy = *this;
        while (nr > 0) {
            auto companion = copy.danilevsky_step(nr);
            nr -= companion.size();
            top_rows.push_back(companion);
        }
        return top_rows;
    }


    constexpr void operator^=(BitMatrix<Word> const& rhs) {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        for (auto i = 0uz; i < rows(); ++i) row(i) ^= rhs.row(i);
    }

    constexpr void operator&=(BitMatrix<Word> const& rhs) {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        for (auto i = 0uz; i < rows(); ++i) row(i) &= rhs.row(i);
    }

    constexpr void operator|=(BitMatrix<Word> const& rhs) {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        for (auto i = 0uz; i < rows(); ++i) row(i) |= rhs.row(i);
    }

    constexpr auto operator^(BitMatrix<Word> const& rhs) const {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        auto result = *this;
        result ^= rhs;
        return result;
    }

    constexpr auto operator&(BitMatrix<Word> const& rhs) const {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        auto result = *this;
        result &= rhs;
        return result;
    }

    constexpr auto operator|(BitMatrix<Word> const& rhs) const {
        gf2_assert(rows() == rhs.rows(), "Row dimensions do not match: {} != {}.", rows(), rhs.rows());
        gf2_assert(cols() == rhs.cols(), "Column dimensions do not match: {} != {}.", cols(), rhs.cols());
        auto result = *this;
        result |= rhs;
        return result;
    }

    constexpr auto operator~() {
        auto result = *this;
        result.flip_all();
        return result;
    }


    constexpr void operator+=(BitMatrix<Word> const& rhs) { operator^=(rhs); }

    constexpr void operator-=(BitMatrix<Word> const& rhs) { operator^=(rhs); }

    constexpr auto operator+(BitMatrix<Word> const& rhs) const {
        auto result = *this;
        result += rhs;
        return result;
    }

    constexpr auto operator-(BitMatrix<Word> const& rhs) const {
        auto result = *this;
        result -= rhs;
        return result;
    }


    std::string to_binary_string(std::string_view row_sep = "\n", std::string_view bit_sep = "",
                                 std::string_view row_prefix = "", std::string_view row_suffix = "") const {
        usize nr = rows();
 
        // Edge case?
        if (nr == 0) return std::string{};
 
        usize       nc = cols();
        usize       per_row = nc + (nc - 1) * bit_sep.size() + row_prefix.size() + row_suffix.size();
        std::string result;
        result.reserve(nr * per_row);
 
        for (auto i = 0uz; i < nr; ++i) {
            result += m_rows[i].to_binary_string(bit_sep, row_prefix, row_suffix);
            if (i + 1 < nr) result += row_sep;
        }
        return result;
    }

    std::string to_compact_binary_string() const { return to_binary_string(" "); }

    std::string to_string() const { return to_binary_string("\n"); }

    std::string to_pretty_string() const { return to_binary_string("\n", " ", "\u2502", "\u2502"); }

    std::string to_hex_string(std::string_view row_sep = "\n") const {
        // Edge case.
        usize nr = rows();
        if (nr == 0) return std::string{};
 
        // Set up the return string with enough space for the entire matrix.
        usize       nc = cols();
        std::string result;
        result.reserve(nr * (nc / 4 + row_sep.size()));
 
        // Append each row to the result string.
        for (auto i = 0uz; i < nr; ++i) {
            result += m_rows[i].to_hex_string();
            if (i + 1 < nr) result += row_sep;
        }
        return result;
    }

    std::string to_compact_hex_string() const { return to_hex_string(" "); }
 
    // The usual output stream operator for a bit-store
    constexpr std::ostream& operator<<(std::ostream& s) const { return s << to_string(); }


    friend constexpr bool operator==(BitMatrix const& lhs, BitMatrix const& rhs) {
        // Edge case.
        if (&lhs == &rhs) return true;
 
        // Check the active words for equality.
        auto row_count = lhs.rows();
        if (rhs.rows() != row_count) return false;
        for (auto i = 0uz; i < row_count; ++i)
            if (lhs.m_rows[i] != rhs.m_rows[i]) return false;
 
        // Made it
        return true;
    }
 
private:
    // Runs a check to see that all the rows in a vector of rows have the same number of columns.
    constexpr bool check_rows(std::vector<row_type> const& rows) const {
        if (rows.empty()) return true;
        auto nc = rows[0].size();
        for (auto i = 1uz; i < rows.size(); ++i) {
            if (rows[i].size() != nc) return false;
        }
        return true;
    }
 
    // Performs a single step of Danilevsky's algorithm to reduce a bit-matrix to Frobenius form.
    //
    // Frobenius form is block diagonal with one or more companion matrices along the diagonal. All matrices can be
    // reduced to this form via a sequence of similarity transformations. This methods performs a single one of those
    // transformations.
    //
    // This `frobenius_form` function calls here with an N x N bit-matrix. In each call the method concentrates on just
    // the top-left n x n sub-matrix. On the first call, n should be set to N. The method returns the top row of the
    // companion matrix that is the transformation of the bottom-right (n-k) x (n-k) sub-matrix. The caller can store
    // that result, decrement n, and call again on the smaller top-left sub-matrix. It may be that the whole matrix
    // gets reduced to a single companion matrix in one step and then there will be no need to call again.
    //
    // The method tries to transform the n x n top-left sub-matrix to a companion matrix working from its bottom-right
    // corner up. It stops when it gets to a point where the bottom-right (n-k) x (n-k) sub-matrix is in companion form
    // and returns the top row of that sub-matrix. The caller can store that result, decrement n, and call again on
    // the smaller top-left sub-matrix.
    //
    // NOTE: This method panics if the bit-matrix is not square.
    BitVector<Word> danilevsky_step(usize n) {
        gf2_assert(n <= rows(), "No top-left {} x {} sub-matrix in a matrix with {} rows", n, n, rows());
 
        // Edge case: A 1 x 1 matrix is already in companion form.
        if (n == 1) return BitVector<Word>::constant(1, get(0, 0));
 
        // Step k of algorithm attempts to reduce row k to companion form.
        // By construction, rows k+1 or later are already in companion form.
        auto k = n - 1;
        while (k > 0) {
            // If row k's sub-diagonal is all zeros we look for an earlier column with a 1.
            // If found, we swap that column here & then swap the equivalent rows to preserve similarity.
            if (!get(k, k - 1)) {
                for (auto j = 0uz; j < k - 1; ++j) {
                    if (get(k, j)) {
                        swap_rows(j, k - 1);
                        swap_cols(j, k - 1);
                        break;
                    }
                }
            }
 
            // No joy? Perhaps we have a companion matrix in the lower left corner and can return its top row?
            if (!get(k, k - 1)) break;
 
            // Still no joy? The sub-diagonal is not all zeros so apply transform to make it so: self <- M^-1 * self *
            // M, where M is the identity matrix with the (k-1)'st row replaced by the k'th row of `self`. We can
            // sparsely represent M as just a clone of that k'th row of `self`.
            auto m = m_rows[k];
 
            // Note the M^-1 is the same as M and self <- M^-1 * self just alters a few of our elements.
            for (auto j = 0uz; j < n; ++j) set(k - 1, j, dot(m, col(j)));
 
            // We also use the sparsity of M when computing self <- self * M.
            for (auto i = 0uz; i < k; ++i) {
                for (auto j = 0uz; j < n; ++j) {
                    auto tmp = get(i, k - 1) && m.get(j);
                    if (j == k - 1) {
                        set(i, j, tmp);
                    } else {
                        set(i, j, get(i, j) ^ tmp);
                    }
                }
            }
 
            // Now put row k into companion form of all zeros with one on the sub-diagonal.
            // All the rows below k are already in companion form.
            m_rows[k].set_all(false);
            set(k, k - 1);
 
            // Done with row k
            k -= 1;
        }
 
        // At this point, k == 0 OR the bit-matrix has non-removable zero on the sub-diagonal of row k.
        // Either way, the bottom-right (n-k) x (n-k) sub-matrix, starting at self[k][k], is in companion form.
        // We return the top row of that companion sub-matrix.
        auto top_row = BitVector<Word>::zeros(n - k);
        for (auto j = 0uz; j < n - k; ++j) top_row.set(j, get(k, k + j));
 
        // Done
        return top_row;
    }
};
 
// --------------------------------------------------------------------------------------------------------------------
// Matrix multiplication ...
// -------------------------------------------------------------------------------------------------------------------

template<Unsigned Word, BitStore Rhs>
    requires std::same_as<typename Rhs::word_type, Word>
constexpr auto
dot(BitMatrix<Word> const& lhs, Rhs const& rhs) {
    gf2_assert_eq(lhs.cols(), rhs.size(), "Incompatible dimensions: {} != {}", lhs.cols(), rhs.size());
    auto n_rows = lhs.rows();
    auto result = BitVector<Word>::zeros(n_rows);
    for (auto i = 0uz; i < n_rows; ++i) {
        if (dot(lhs.row(i), rhs)) result.set(i, true);
    }
    return result;
}

template<Unsigned Word, BitStore Rhs>
    requires std::same_as<typename Rhs::word_type, Word>
constexpr auto
operator*(BitMatrix<Word> const& lhs, Rhs const& rhs) {
    return dot(lhs, rhs);
}

template<Unsigned Word, BitStore Lhs>
    requires std::same_as<typename Lhs::word_type, Word>
constexpr auto
dot(Lhs const& lhs, BitMatrix<Word> const& rhs) {
    gf2_assert_eq(lhs.size(), rhs.rows(), "Incompatible dimensions: {} != {}", lhs.size(), rhs.rows());
    auto n_cols = rhs.cols();
    auto result = BitVector<Word>::zeros(n_cols);
    for (auto j = 0uz; j < n_cols; ++j) {
        if (dot(lhs, rhs.col(j))) result.set(j, true);
    }
    return result;
}

template<Unsigned Word, BitStore Lhs>
    requires std::same_as<typename Lhs::word_type, Word>
constexpr auto
operator*(Lhs const& lhs, BitMatrix<Word> const& rhs) {
    return dot(lhs, rhs);
}

template<Unsigned Word>
constexpr auto
dot(BitMatrix<Word> const& lhs, BitMatrix<Word> const& rhs) {
    gf2_assert_eq(lhs.cols(), rhs.rows(), "Incompatible dimensions: {} != {}", lhs.cols(), rhs.rows());
 
    auto n_rows = lhs.rows();
    auto n_cols = rhs.cols();
    auto result = BitMatrix<Word>::zeros(n_rows, n_cols);
 
    // Row access is cheap, columns expensive, so arrange things to pull out columns as few times as possible.
    for (auto j = 0uz; j < n_cols; ++j) {
        auto rhs_col = rhs.col(j);
        for (auto i = 0uz; i < n_rows; ++i) {
            if (dot(lhs.row(i), rhs_col)) result.set(i, j, true);
        }
    }
    return result;
}

template<Unsigned Word>
constexpr auto
operator*(BitMatrix<Word> const& lhs, BitMatrix<Word> const& rhs) {
    return dot(lhs, rhs);
}
 
// --------------------------------------------------------------------------------------------------------------------
// Some utility methods to print multiple matrices & vectors side-by-side ...
// -------------------------------------------------------------------------------------------------------------------

template<Unsigned Word, BitStore Rhs>
    requires std::same_as<typename Rhs::word_type, Word>
std::string
string_for(BitMatrix<Word> const& A, Rhs const& b) {
 
    // If either is empty there was likely a bug somewhere!
    gf2_assert(!A.is_empty(), "Matrix A is empty which is likely an error!");
    gf2_assert(!b.is_empty(), "Vector b is empty which is likely an error!");
 
    // Most rows we ever print
    auto nr = std::max(A.rows(), b.size());
 
    // Space filler might be needed if the matrix has fewer rows than the vector
    std::string A_fill(A.cols(), ' ');
 
    // Little lambda that turns a bool into a string
    auto b2s = [](bool x) { return x ? "1" : "0"; };
 
    // Build the result string
    std::string result;
    for (auto r = 0uz; r < nr; ++r) {
        auto A_str = r < A.rows() ? A[r].to_string() : A_fill;
        auto b_str = r < b.size() ? b2s(b[r]) : " ";
        result += A_str + "\t" + b_str;
        if (r + 1 < nr) result += "\n";
    }
 
    return result;
}

template<Unsigned Word, BitStore Rhs>
    requires std::same_as<typename Rhs::word_type, Word>
std::string
string_for(BitMatrix<Word> const& A, Rhs const& b, Rhs const& c) {
 
    // If either is empty there was likely a bug somewhere!
    gf2_assert(!A.is_empty(), "Matrix A is empty which is likely an error!");
    gf2_assert(!b.is_empty(), "Vector b is empty which is likely an error!");
    gf2_assert(!c.is_empty(), "Vector c is empty which is likely an error!");
 
    // Most rows we ever print
    auto nr = std::max({A.rows(), b.size(), c.size()});
 
    // Space filler might be needed if the matrix has fewer rows than the vector
    std::string A_fill(A.cols(), ' ');
 
    // Little lambda that turns a bool into a string
    auto b2s = [](bool x) { return x ? "1" : "0"; };
 
    // Build the result string
    std::string result;
    for (auto r = 0uz; r < nr; ++r) {
        auto A_str = r < A.rows() ? A[r].to_string() : A_fill;
        auto b_str = r < b.size() ? b2s(b[r]) : " ";
        auto c_str = r < c.size() ? b2s(c[r]) : " ";
        result += A_str + "\t" + b_str + "\t" + c_str;
        if (r + 1 < nr) result += "\n";
    }
 
    return result;
}

template<Unsigned Word, BitStore Rhs>
    requires std::same_as<typename Rhs::word_type, Word>
std::string
string_for(BitMatrix<Word> const& A, Rhs const& b, Rhs const& c, Rhs const& d) {
 
    // If either is empty there was likely a bug somewhere!
    gf2_assert(!A.is_empty(), "Matrix A is empty which is likely an error!");
    gf2_assert(!b.is_empty(), "Vector b is empty which is likely an error!");
    gf2_assert(!c.is_empty(), "Vector c is empty which is likely an error!");
    gf2_assert(!d.is_empty(), "Vector d is empty which is likely an error!");
 
    // Most rows we ever print
    auto nr = std::max({A.rows(), b.size(), c.size(), d.size()});
 
    // Space filler might be needed if the matrix has fewer rows than the vector
    std::string A_fill(A.cols(), ' ');
 
    // Little lambda that turns a bool into a string
    auto b2s = [](bool x) { return x ? "1" : "0"; };
 
    // Build the result string
    std::string result;
    for (auto r = 0uz; r < nr; ++r) {
        auto A_str = r < A.rows() ? A[r].to_string() : A_fill;
        auto b_str = r < b.size() ? b2s(b[r]) : " ";
        auto c_str = r < c.size() ? b2s(c[r]) : " ";
        auto d_str = r < d.size() ? b2s(d[r]) : " ";
        result += A_str + "\t" + b_str + "\t" + c_str + "\t" + d_str;
        if (r + 1 < nr) result += "\n";
    }
 
    return result;
}

template<Unsigned Word>
std::string
string_for(BitMatrix<Word> const& A, BitMatrix<Word> const& B) {
 
    // If either is empty there was likely a bug somewhere!
    gf2_assert(!A.is_empty(), "Matrix A is empty which is likely an error!");
    gf2_assert(!B.is_empty(), "Matrix B is empty which is likely an error!");
 
    // Most rows we ever print
    auto nr = std::max(A.rows(), B.rows());
 
    // Space filler might be needed if the matrix has fewer rows than the vector
    std::string A_fill(A.cols(), ' ');
    std::string B_fill(B.cols(), ' ');
 
    // Build the result string
    std::string result;
    for (auto r = 0uz; r < nr; ++r) {
        auto A_str = r < A.rows() ? A[r].to_string() : A_fill;
        auto B_str = r < B.rows() ? B[r].to_string() : B_fill;
        result += A_str + "\t" + B_str;
        if (r + 1 < nr) result += "\n";
    }
 
    return result;
}

template<Unsigned Word>
std::string
string_for(BitMatrix<Word> const& A, BitMatrix<Word> const& B, BitMatrix<Word> const& C) {
 
    // If either is empty there was likely a bug somewhere!
    gf2_assert(!A.is_empty(), "Matrix A is empty which is likely an error!");
    gf2_assert(!B.is_empty(), "Matrix B is empty which is likely an error!");
    gf2_assert(!C.is_empty(), "Matrix C is empty which is likely an error!");
 
    // Most rows we ever print
    auto nr = std::max({A.rows(), B.rows(), C.rows()});
 
    // Space filler might be needed if the matrix has fewer rows than the vector
    std::string A_fill(A.cols(), ' ');
    std::string B_fill(B.cols(), ' ');
    std::string C_fill(C.cols(), ' ');
 
    // Build the result string
    std::string result;
    for (auto r = 0uz; r < nr; ++r) {
        auto A_str = r < A.rows() ? A[r].to_string() : A_fill;
        auto B_str = r < B.rows() ? B[r].to_string() : B_fill;
        auto C_str = r < C.rows() ? C[r].to_string() : C_fill;
        result += A_str + "\t" + B_str + "\t" + C_str;
        if (r + 1 < nr) result += "\n";
    }
 
    return result;
}
 
} // namespace gf2
 
// --------------------------------------------------------------------------------------------------------------------
// Specialises `std::formatter` to handle bit-matrices ...
// -------------------------------------------------------------------------------------------------------------------

template<gf2::Unsigned Word>
struct std::formatter<gf2::BitMatrix<Word>> {

    constexpr auto parse(std::format_parse_context const& ctx) {
        auto it = ctx.begin();
        while (it != ctx.end() && *it != '}') {
            switch (*it) {
                case 'p': m_pretty = true; break;
                case 'x': m_hex = true; break;
                default: m_error = true;
            }
            ++it;
        }
        return it;
    }

    template<class FormatContext>
    auto format(gf2::BitMatrix<Word> const& rhs, FormatContext& ctx) const {
        // Was there a format specification error?
        if (m_error) return std::format_to(ctx.out(), "'UNRECOGNIZED FORMAT SPECIFIER FOR BIT-MATRIX'");
 
        // Special handling requested?
        if (m_hex) return std::format_to(ctx.out(), "{}", rhs.to_hex_string());
        if (m_pretty) return std::format_to(ctx.out(), "{}", rhs.to_pretty_string());
 
        // Default
        return std::format_to(ctx.out(), "{}", rhs.to_string());
    }
 
    bool m_hex = false;
    bool m_pretty = false;
    bool m_error = false;
};