BitLU.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/BitMatrix.h>
 
namespace gf2 {

template<Unsigned Word>
class BitLU {
private:
    BitMatrix<Word>    m_lu;    // The LU decomposition stored in a single bit-matrix.
    std::vector<usize> m_swaps; // The list of row swap instructions (a permutation in LAPACK format).
    usize              m_rank;  // The rank of the underlying matrix.
 
public:
    BitLU(BitMatrix<Word> const& A) : m_lu(A), m_swaps(A.rows(), 0uz), m_rank(A.rows()) {
        // Only handle square matrices
        gf2_assert(A.is_square(), "Matrix is {} x {} but it should be square!", A.rows(), A.cols());
 
        // Iterate through the bit-matrix.
        auto nr = m_lu.rows();
        auto nc = m_lu.cols();
        for (auto j = 0uz; j < nr; ++j) {
 
            // Initialize this element of the row swap instruction vector.
            m_swaps[j] = j;
 
            // Find a non-zero element on or below the diagonal in column -- a pivot.
            auto p = j;
            while (p < nr && !m_lu.get(p, j)) ++p;
 
            // Perhaps no such element exists in this column? If so the matrix is singular!
            if (p >= nr) {
                m_rank--;
                continue;
            }
 
            // If necessary, swap the pivot row  into place & record the rows swap instruction
            if (p != j) {
                m_lu.swap_rows(j, p);
                m_swaps[j] = p;
            }
 
            // At this point m_lu(j,j) == 1
            for (auto i = j + 1; i < nr; ++i) {
                if (m_lu.get(i, j)) {
                    for (auto k = j + 1; k < nc; ++k) {
                        auto tmp = m_lu.get(i, k) ^ m_lu.get(j, k);
                        m_lu.set(i, k, tmp);
                    }
                }
            }
        }
    }

    constexpr usize rank() const { return m_rank; }

    constexpr bool is_singular() const { return m_rank < m_lu.rows(); }

    constexpr bool determinant() const { return !is_singular(); }

    constexpr const BitMatrix<Word>& LU() const { return m_lu; }

    constexpr BitMatrix<Word> L() const { return m_lu.unit_lower(); }

    constexpr BitMatrix<Word> U() const { return m_lu.upper(); }

    constexpr const std::vector<usize>& swaps() const { return m_swaps; }

    std::vector<usize> permutation_vector() const {
        auto n = m_swaps.size();
 
        std::vector<usize> result(n);
        for (auto i = 0uz; i < n; ++i) result[i] = i;
        for (auto i = 0uz; i < n; ++i) std::swap(result[m_swaps[i]], result[i]);
        return result;
    }

    constexpr void permute(BitMatrix<Word>& B) const {
        auto n = m_swaps.size();
        gf2_assert(B.rows() == n, "Matrix has {} rows but the row-swap instruction vector has {}.", B.rows(), n);
        for (auto i = 0uz; i < n; ++i)
            if (m_swaps[i] != i) B.swap_rows(i, m_swaps[i]);
    }

    template<BitStore Store>
        requires(std::same_as<typename Store::word_type, Word>)
    constexpr void permute(Store& b) const {
        auto n = m_swaps.size();
        gf2_assert(b.size() == n, "Vector has {} elements but the row-swap instruction vector has {}.", b.size(), n);
        for (auto i = 0uz; i < n; ++i)
            if (m_swaps[i] != i) b.swap(i, m_swaps[i]);
    }

    template<BitStore Store>
        requires(std::same_as<typename Store::word_type, Word>)
    std::optional<BitVector<Word>> operator()(Store const& b) const {
        auto n = m_lu.rows();
        gf2_assert(b.size() == n, "RHS b has {} elements but the LHS matrix has {} rows.", b.size(), n);
 
        // Can we solve this at all?
        if (is_singular()) return std::nullopt;
 
        // Make a permuted copy of the the right hand side
        auto x = b;
        permute(x);
 
        // Forward substitution
        for (auto i = 0uz; i < n; ++i) {
            for (auto j = 0uz; j < i; ++j)
                if (m_lu(i, j)) x[i] ^= x[j];
        }
 
        // Backward substitution
        for (auto i = n; i--;) {
            for (auto j = i + 1; j < n; ++j)
                if (m_lu(i, j)) x[i] ^= x[j];
        }
 
        return x;
    }

    std::optional<BitMatrix<Word>> operator()(BitMatrix<Word> const& B) const {
        auto n = m_lu.rows();
        gf2_assert(B.rows() == n, "RHS B has {} rows but the LHS matrix has {} rows.", B.rows(), n);
 
        // Can we solve this at all?
        if (is_singular()) return std::nullopt;
 
        // Make a permuted copy of the the right hand side
        auto X = B;
        permute(X);
 
        // Solve for each column
        for (auto j = 0uz; j < n; ++j) {
 
            // Forward substitution
            for (auto i = 0uz; i < n; ++i) {
                for (auto k = 0uz; k < i; ++k)
                    if (m_lu(i, k)) X(i, j) ^= X(k, j);
            }
 
            // Backward substitution
            for (auto i = n; i--;) {
                for (auto k = i + 1; k < n; ++k)
                    if (m_lu(i, k)) X(i, j) ^= X(k, j);
            }
        }
 
        return X;
    }

    std::optional<BitMatrix<Word>> inverse() const {
        // We just solve the system A.A_inv = I for A_inv
        return operator()(BitMatrix<Word>::identity(m_lu.rows()));
    }
};
} // namespace gf2