BitGauss.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/BitVec.h>
#include <gf2/BitMat.h>
 
namespace gf2 {

template<unsigned_word Word>
class BitGauss {
private:
    BitMat<Word>       m_A;         // The reduced row echelon form of the matrix `A`.
    BitVec<Word>       m_b;         // The equivalent reduced row echelon form of the vector `b`.
    usize              m_rank;      // The rank of the matrix `A`.
    usize              m_solutions; // The number of solutions to the system.
    std::vector<usize> m_free;      // The indices of the free variables in the system.
 
public:
    template<typename Store>
    BitGauss(const BitMat<Word>& A, const BitStore<Store>& b)
        requires store_word_is<Store, Word>
    {
        gf2_assert(A.is_square(), "Matrix is {} x {} but it should be square!", A.rows(), A.cols());
        gf2_assert(A.rows() == b.size(), "Matrix has {} rows but the RHS vector has {} elements.", A.rows(), b.size());
 
        // Copy A & augment it with b on the right.
        m_A = A;
        m_A.append_col(b);
 
        // Reduce the augmented matrix to reduced row echelon form.
        auto has_pivot = m_A.to_reduced_echelon_form();
 
        // Grab the last column of the reduced augmented matrix as a standalone vector, removing it from the matrix.
        m_b = m_A.remove_col().value();
 
        // Also pop the last element of the `has_pivot` bit-vector as it corresponds to the column that we just removed.
        has_pivot.pop();
 
        // The rank is the number of columns with a pivot (it is also the number if non-zero rows in the reduced
        // matrix).
        m_rank = has_pivot.count_ones();
 
        // Any columns without a pivot are free variables.
        for (auto i = 0uz; i < has_pivot.size(); ++i) {
            if (!has_pivot[i]) m_free.push_back(i);
        }
 
        // Check that all zero rows in the reduced matrix are matched by zero elements in the equivalent reduced right
        // hand side vector (this is consistency). Zero rows are at the bottom from `rank` to the end.
        auto consistent = true;
        for (auto i = m_rank; i < m_A.rows(); ++i) {
            if (m_b[i]) {
                consistent = false;
                break;
            }
        }
 
        // The number of solutions we can index into is either 0 or 2^f where f is the number of free variables.
        // However, for practical reasons we limit the number of solutions to `min(2^f, 2^63)`.
        m_solutions = 0;
        if (consistent) {
            auto f = m_free.size();
            auto b_max = static_cast<usize>(BITS<Word> - 1);
            auto f_max = std::min(f, b_max);
            m_solutions = 1uz << f_max;
        }
    }

    constexpr usize rank() const { return m_rank; }

    constexpr usize free_count() const { return m_free.size(); }

    constexpr bool is_underdetermined() const { return !m_free.empty(); }

    constexpr bool is_consistent() const { return m_solutions > 0; }

    constexpr usize solution_count() const { return m_solutions; }

    std::optional<BitVec<Word>> operator()() const {
        if (!is_consistent()) return std::nullopt;
 
        // Create a random starting point.
        auto x = BitVec<Word>::random(m_b.size());
 
        // All non-free variables will be overwritten by back substitution.
        back_substitute_into(x);
        return x;
    }

    std::optional<BitVec<Word>> operator()(usize i_solution) const {
        if (!is_consistent()) { return std::nullopt; }
        if (i_solution > solution_count()) { return std::nullopt; }
 
        // We start with a zero vector and then set the free variable slots to the fixed bit pattern for `i`.
        auto x = BitVec<Word>::zeros(m_b.size());
        for (auto f = 0uz; f < m_free.size(); ++f) {
            x[m_free[f]] = (i_solution & 1);
            i_solution >>= 1;
        }
 
        // Back substitution will now overwrite any non-free variables in `x` with their correct values.
        back_substitute_into(x);
        return x;
    }
 
private:
    // Helper function that performs back substitution to solve for the non-free variables in `x`.
    void back_substitute_into(BitVec<Word>& x) const {
        // Iterate from the bottom up, starting at the first non-zero row, solving for the non-free variables in `x`.
        for (auto i = m_rank; i--;) {
            auto j = m_A[i].first_set().value();
            x.set(j, m_b[i]);
            for (auto k = j + 1; k < x.size(); ++k)
                if (m_A(i, k)) { x.set(j, x[j] ^ x[k]); }
        }
    }
};
 
} // namespace gf2