m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

View on GitHub

:heavy_check_mark: Matrix Bundle
(math/matrix/all.hpp)

Overview

math/matrix/all.hpp includes the complete dense and packed GF(2) matrix module.

Included Headers

Header Contents
math/matrix/adjugate.hpp Cubic-time adjugate matrix over a field, including singular matrices.
math/matrix/bit_matrix.hpp Packed GF(2) matrices, arithmetic, multiplication, elimination, rank, inverse, and linear systems.
math/matrix/characteristic_polynomial.hpp Characteristic polynomial of a square matrix over a field.
math/matrix/determinant_mod.hpp Determinant modulo an arbitrary positive, possibly composite modulus.
math/matrix/pfaffian.hpp Cubic-time Pfaffian of an alternating matrix.
math/matrix/hafnian.hpp Exact hafnian of a small symmetric matrix.
math/matrix/sparse_determinant.hpp Randomized black-box determinant of a sparse matrix over a finite field.
math/matrix/matrix.hpp Row-major dense matrices, arithmetic, multiplication, transposition, matrix-vector products, and powers.
math/matrix/linear_algebra.hpp Gaussian elimination, rank, determinant, inverse, and linear systems.

Depends on

Required by

Verified with

Code

#ifndef M1UNE_MATRIX_ALL_HPP
#define M1UNE_MATRIX_ALL_HPP 1

#include "adjugate.hpp"
#include "bit_matrix.hpp"
#include "characteristic_polynomial.hpp"
#include "determinant_mod.hpp"
#include "hafnian.hpp"
#include "linear_algebra.hpp"
#include "matrix.hpp"
#include "pfaffian.hpp"
#include "sparse_determinant.hpp"

#endif  // M1UNE_MATRIX_ALL_HPP
#line 1 "math/matrix/all.hpp"



#line 1 "math/matrix/adjugate.hpp"



#include <cassert>
#include <vector>

#line 1 "math/matrix/matrix.hpp"



#line 5 "math/matrix/matrix.hpp"
#include <cstddef>
#include <cstdint>
#include <utility>
#line 9 "math/matrix/matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
class Matrix {
   private:
    int _rows;
    int _cols;
    std::vector<T> _data;

    static std::size_t storage_size(int rows, int cols) {
        assert(rows >= 0);
        assert(cols >= 0);
        return std::size_t(rows) * std::size_t(cols);
    }

   public:
    using value_type = T;

    Matrix() : _rows(0), _cols(0) {}

    Matrix(int rows, int cols, const T& value = T())
        : _rows(rows), _cols(cols), _data(storage_size(rows, cols), value) {}

    Matrix(int rows, int cols, std::vector<T> values)
        : _rows(rows), _cols(cols), _data(std::move(values)) {
        assert(rows >= 0);
        assert(cols >= 0);
        assert(_data.size() == std::size_t(rows) * std::size_t(cols));
    }

    explicit Matrix(const std::vector<std::vector<T>>& values)
        : _rows(int(values.size())), _cols(values.empty() ? 0 : int(values[0].size())),
          _data(storage_size(_rows, _cols)) {
        for (int row = 0; row < _rows; row++) {
            assert(int(values[std::size_t(row)].size()) == _cols);
            for (int col = 0; col < _cols; col++) {
                (*this)[row][col] = values[std::size_t(row)][std::size_t(col)];
            }
        }
    }

    int rows() const {
        return _rows;
    }

    int cols() const {
        return _cols;
    }

    bool empty() const {
        return _rows == 0 || _cols == 0;
    }

    std::vector<T>& data() {
        return _data;
    }

    const std::vector<T>& data() const {
        return _data;
    }

    T* operator[](int row) {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    const T* operator[](int row) const {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    T& operator()(int row, int col) {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    const T& operator()(int row, int col) const {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    static Matrix identity(int size) {
        assert(size >= 0);
        Matrix result(size, size);
        for (int i = 0; i < size; i++) result[i][i] = T(1);
        return result;
    }

    Matrix transposed() const {
        Matrix result(_cols, _rows);
        for (int row = 0; row < _rows; row++) {
            for (int col = 0; col < _cols; col++) {
                result[col][row] = (*this)[row][col];
            }
        }
        return result;
    }

    void swap_rows(int first, int second) {
        assert(0 <= first && first < _rows);
        assert(0 <= second && second < _rows);
        if (first == second) return;
        for (int col = 0; col < _cols; col++) {
            std::swap((*this)[first][col], (*this)[second][col]);
        }
    }

    Matrix& operator+=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] += rhs._data[i];
        return *this;
    }

    Matrix& operator-=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] -= rhs._data[i];
        return *this;
    }

    Matrix& operator*=(const T& scalar) {
        for (T& value : _data) value *= scalar;
        return *this;
    }

    Matrix& operator/=(const T& scalar) {
        for (T& value : _data) value /= scalar;
        return *this;
    }

    Matrix& operator*=(const Matrix& rhs) {
        return *this = *this * rhs;
    }

    Matrix operator+() const {
        return *this;
    }

    Matrix operator-() const {
        Matrix result = *this;
        for (T& value : result._data) value = T() - value;
        return result;
    }

    friend Matrix operator+(Matrix lhs, const Matrix& rhs) {
        return lhs += rhs;
    }

    friend Matrix operator-(Matrix lhs, const Matrix& rhs) {
        return lhs -= rhs;
    }

    friend Matrix operator*(Matrix lhs, const T& rhs) {
        return lhs *= rhs;
    }

    friend Matrix operator*(const T& lhs, Matrix rhs) {
        return rhs *= lhs;
    }

    friend Matrix operator/(Matrix lhs, const T& rhs) {
        return lhs /= rhs;
    }

    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
        assert(lhs._cols == rhs._rows);
        Matrix result(lhs._rows, rhs._cols);
        for (int row = 0; row < lhs._rows; row++) {
            T* output = result[row];
            for (int middle = 0; middle < lhs._cols; middle++) {
                const T coefficient = lhs[row][middle];
                if (coefficient == T()) continue;
                const T* input = rhs[middle];
                for (int col = 0; col < rhs._cols; col++) {
                    output[col] += coefficient * input[col];
                }
            }
        }
        return result;
    }

    friend std::vector<T> operator*(const Matrix& lhs, const std::vector<T>& rhs) {
        assert(lhs._cols == int(rhs.size()));
        std::vector<T> result(std::size_t(lhs._rows));
        for (int row = 0; row < lhs._rows; row++) {
            T value = T();
            for (int col = 0; col < lhs._cols; col++) {
                value += lhs[row][col] * rhs[std::size_t(col)];
            }
            result[std::size_t(row)] = value;
        }
        return result;
    }

    friend std::vector<T> operator*(const std::vector<T>& lhs, const Matrix& rhs) {
        assert(int(lhs.size()) == rhs._rows);
        std::vector<T> result(std::size_t(rhs._cols));
        for (int row = 0; row < rhs._rows; row++) {
            if (lhs[std::size_t(row)] == T()) continue;
            for (int col = 0; col < rhs._cols; col++) {
                result[std::size_t(col)] += lhs[std::size_t(row)] * rhs[row][col];
            }
        }
        return result;
    }

    bool operator==(const Matrix& rhs) const {
        return _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
    }

    bool operator!=(const Matrix& rhs) const {
        return !(*this == rhs);
    }

    Matrix pow(std::uint64_t exponent) const {
        assert(_rows == _cols);
        Matrix result = identity(_rows);
        Matrix base = *this;
        while (exponent > 0) {
            if (exponent & 1) result *= base;
            exponent >>= 1;
            if (exponent > 0) base *= base;
        }
        return result;
    }
};

}  // namespace matrix
}  // namespace m1une


#line 8 "math/matrix/adjugate.hpp"

namespace m1une {
namespace matrix {

template <class T>
Matrix<T> adjugate(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    std::vector<int> pivot_columns;
    T pivot_product = T(1);
    bool negate = false;
    for (int col = 0; col < size && int(pivot_columns.size()) < size; col++) {
        const int pivot_row = int(pivot_columns.size());
        int pivot = pivot_row;
        while (pivot < size && augmented[pivot][col] == T()) pivot++;
        if (pivot == size) continue;
        if (pivot != pivot_row) {
            augmented.swap_rows(pivot, pivot_row);
            negate = !negate;
        }

        const T pivot_value = augmented[pivot_row][col];
        pivot_product *= pivot_value;
        const T inverse_pivot = T(1) / pivot_value;
        for (int index = col; index < size; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }
        for (int index = size; index < size * 2; index++) {
            augmented[pivot_row][index] *= inverse_pivot;
        }

        for (int row = 0; row < size; row++) {
            if (row == pivot_row || augmented[row][col] == T()) continue;
            const T factor = augmented[row][col];
            augmented[row][col] = T();
            for (int index = col + 1; index < size; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
            for (int index = size; index < size * 2; index++) {
                augmented[row][index] -= factor * augmented[pivot_row][index];
            }
        }
        pivot_columns.push_back(col);
    }

    const int rank = int(pivot_columns.size());
    Matrix<T> result(size, size);
    if (rank + 1 < size) return result;

    if (rank == size) {
        const T determinant = negate ? T() - pivot_product : pivot_product;
        for (int row = 0; row < size; row++) {
            for (int col = 0; col < size; col++) {
                result[row][col] = determinant * augmented[row][size + col];
            }
        }
        return result;
    }

    int free_column = 0;
    while (free_column < rank && pivot_columns[free_column] == free_column) {
        free_column++;
    }
    std::vector<T> right_null(size);
    right_null[free_column] = T(1);
    for (int row = 0; row < rank; row++) {
        right_null[pivot_columns[row]] = T() - augmented[row][free_column];
    }

    T scale = pivot_product;
    if (negate != bool((size - 1 + free_column) & 1)) scale = T() - scale;
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] =
                scale * right_null[row] * augmented[size - 1][size + col];
        }
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/bit_matrix.hpp"



#include <algorithm>
#include <bit>
#line 9 "math/matrix/bit_matrix.hpp"
#include <optional>
#include <string>
#include <string_view>
#line 14 "math/matrix/bit_matrix.hpp"

namespace m1une {
namespace matrix {

class BitMatrix {
   private:
    int _rows;
    int _cols;
    int _blocks;
    std::vector<std::uint64_t> _data;

    static int block_count(int cols) {
        assert(cols >= 0);
        return (cols + 63) / 64;
    }

    static std::size_t storage_size(int rows, int blocks) {
        assert(rows >= 0);
        return std::size_t(rows) * std::size_t(blocks);
    }

    std::size_t word_index(int row, int col) const {
        assert(0 <= row && row < _rows);
        assert(0 <= col && col < _cols);
        return std::size_t(row) * std::size_t(_blocks) +
               std::size_t(col / 64);
    }

    std::uint64_t trailing_mask() const {
        if ((_cols & 63) == 0) return ~std::uint64_t(0);
        return (std::uint64_t(1) << (_cols & 63)) - 1;
    }

   public:
    class BitReference {
       private:
        std::uint64_t* word;
        std::uint64_t mask;

       public:
        BitReference(std::uint64_t& word_value, std::uint64_t mask_value)
            : word(&word_value), mask(mask_value) {}

        operator bool() const {
            return (*word & mask) != 0;
        }

        BitReference& operator=(bool value) {
            if (value) {
                *word |= mask;
            } else {
                *word &= ~mask;
            }
            return *this;
        }

        BitReference& operator=(const BitReference& other) {
            return *this = bool(other);
        }

        void flip() {
            *word ^= mask;
        }
    };

    class RowReference {
       private:
        BitMatrix* matrix;
        int row;

       public:
        RowReference(BitMatrix& matrix_value, int row_value)
            : matrix(&matrix_value), row(row_value) {}

        BitReference operator[](int col) const {
            return (*matrix)(row, col);
        }
    };

    class ConstRowReference {
       private:
        const BitMatrix* matrix;
        int row;

       public:
        ConstRowReference(const BitMatrix& matrix_value, int row_value)
            : matrix(&matrix_value), row(row_value) {}

        bool operator[](int col) const {
            return (*matrix)(row, col);
        }
    };

    BitMatrix() : _rows(0), _cols(0), _blocks(0) {}

    BitMatrix(int rows, int cols, bool value = false)
        : _rows(rows),
          _cols(cols),
          _blocks(block_count(cols)),
          _data(
              storage_size(rows, _blocks),
              value ? ~std::uint64_t(0) : std::uint64_t(0)
          ) {
        assert(rows >= 0);
        if (value && _blocks > 0) {
            const std::uint64_t mask = trailing_mask();
            for (int row = 0; row < _rows; row++) {
                _data[
                    std::size_t(row + 1) * std::size_t(_blocks) - 1
                ] &= mask;
            }
        }
    }

    int rows() const {
        return _rows;
    }

    int cols() const {
        return _cols;
    }

    int blocks_per_row() const {
        return _blocks;
    }

    bool empty() const {
        return _rows == 0 || _cols == 0;
    }

    RowReference operator[](int row) {
        assert(0 <= row && row < _rows);
        return RowReference(*this, row);
    }

    ConstRowReference operator[](int row) const {
        assert(0 <= row && row < _rows);
        return ConstRowReference(*this, row);
    }

    BitReference operator()(int row, int col) {
        const std::size_t index = word_index(row, col);
        return BitReference(_data[index], std::uint64_t(1) << (col & 63));
    }

    bool operator()(int row, int col) const {
        const std::size_t index = word_index(row, col);
        return (_data[index] >> (col & 63)) & 1;
    }

    bool get(int row, int col) const {
        return (*this)(row, col);
    }

    void set(int row, int col, bool value = true) {
        (*this)(row, col) = value;
    }

    void reset(int row, int col) {
        set(row, col, false);
    }

    void flip(int row, int col) {
        (*this)(row, col).flip();
    }

    void clear() {
        std::fill(_data.begin(), _data.end(), std::uint64_t(0));
    }

    void set_row(int row, std::string_view bits) {
        assert(0 <= row && row < _rows);
        assert(int(bits.size()) == _cols);
        const std::size_t offset =
            std::size_t(row) * std::size_t(_blocks);
        std::fill(
            _data.begin() + std::ptrdiff_t(offset),
            _data.begin() + std::ptrdiff_t(offset + std::size_t(_blocks)),
            std::uint64_t(0)
        );
        for (int col = 0; col < _cols; col++) {
            assert(bits[std::size_t(col)] == '0' || bits[std::size_t(col)] == '1');
            if (bits[std::size_t(col)] == '1') set(row, col);
        }
    }

    std::string row_string(int row) const {
        assert(0 <= row && row < _rows);
        std::string result(std::size_t(_cols), '0');
        for (int col = 0; col < _cols; col++) {
            if (get(row, col)) result[std::size_t(col)] = '1';
        }
        return result;
    }

    static BitMatrix identity(int size) {
        assert(size >= 0);
        BitMatrix result(size, size);
        for (int index = 0; index < size; index++) result.set(index, index);
        return result;
    }

    BitMatrix transposed() const {
        BitMatrix result(_cols, _rows);
        for (int row = 0; row < _rows; row++) {
            for (int col = 0; col < _cols; col++) {
                if (get(row, col)) result.set(col, row);
            }
        }
        return result;
    }

    void swap_rows(int first, int second) {
        assert(0 <= first && first < _rows);
        assert(0 <= second && second < _rows);
        if (first == second) return;
        const std::size_t first_offset =
            std::size_t(first) * std::size_t(_blocks);
        const std::size_t second_offset =
            std::size_t(second) * std::size_t(_blocks);
        for (int block = 0; block < _blocks; block++) {
            std::swap(
                _data[first_offset + std::size_t(block)],
                _data[second_offset + std::size_t(block)]
            );
        }
    }

    void xor_rows(int target, int source, int first_col = 0) {
        assert(0 <= target && target < _rows);
        assert(0 <= source && source < _rows);
        assert(0 <= first_col && first_col <= _cols);
        if (first_col == _cols) return;
        const std::size_t target_offset =
            std::size_t(target) * std::size_t(_blocks);
        const std::size_t source_offset =
            std::size_t(source) * std::size_t(_blocks);
        const int first_block = first_col / 64;
        const int first_bit = first_col & 63;
        if (first_bit != 0) {
            const std::uint64_t mask = ~std::uint64_t(0) << first_bit;
            _data[target_offset + std::size_t(first_block)] ^=
                _data[source_offset + std::size_t(first_block)] & mask;
        } else {
            _data[target_offset + std::size_t(first_block)] ^=
                _data[source_offset + std::size_t(first_block)];
        }
        for (int block = first_block + 1; block < _blocks; block++) {
            _data[target_offset + std::size_t(block)] ^=
                _data[source_offset + std::size_t(block)];
        }
    }

    BitMatrix& operator^=(const BitMatrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t index = 0; index < _data.size(); index++) {
            _data[index] ^= rhs._data[index];
        }
        return *this;
    }

    BitMatrix& operator+=(const BitMatrix& rhs) {
        return *this ^= rhs;
    }

    BitMatrix& operator-=(const BitMatrix& rhs) {
        return *this ^= rhs;
    }

    BitMatrix& operator*=(const BitMatrix& rhs) {
        return *this = *this * rhs;
    }

    friend BitMatrix operator^(BitMatrix lhs, const BitMatrix& rhs) {
        return lhs ^= rhs;
    }

    friend BitMatrix operator+(BitMatrix lhs, const BitMatrix& rhs) {
        return lhs += rhs;
    }

    friend BitMatrix operator-(BitMatrix lhs, const BitMatrix& rhs) {
        return lhs -= rhs;
    }

    friend BitMatrix operator*(const BitMatrix& lhs, const BitMatrix& rhs) {
        assert(lhs._cols == rhs._rows);
        BitMatrix result(lhs._rows, rhs._cols);
        for (int row = 0; row < lhs._rows; row++) {
            const std::size_t lhs_offset =
                std::size_t(row) * std::size_t(lhs._blocks);
            const std::size_t result_offset =
                std::size_t(row) * std::size_t(result._blocks);
            for (int lhs_block = 0; lhs_block < lhs._blocks; lhs_block++) {
                std::uint64_t word =
                    lhs._data[lhs_offset + std::size_t(lhs_block)];
                while (word != 0) {
                    const int bit = std::countr_zero(word);
                    const int middle = lhs_block * 64 + bit;
                    const std::size_t rhs_offset =
                        std::size_t(middle) * std::size_t(rhs._blocks);
                    for (int block = 0; block < rhs._blocks; block++) {
                        result._data[result_offset + std::size_t(block)] ^=
                            rhs._data[rhs_offset + std::size_t(block)];
                    }
                    word &= word - 1;
                }
            }
        }
        return result;
    }

    bool operator==(const BitMatrix& rhs) const {
        return
            _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
    }

    bool operator!=(const BitMatrix& rhs) const {
        return !(*this == rhs);
    }

    BitMatrix pow(std::uint64_t exponent) const {
        assert(_rows == _cols);
        BitMatrix result = identity(_rows);
        BitMatrix base = *this;
        while (exponent > 0) {
            if (exponent & 1) result *= base;
            exponent >>= 1;
            if (exponent > 0) base *= base;
        }
        return result;
    }
};

namespace bit_matrix_detail {

inline std::vector<int> row_reduce(
    BitMatrix& matrix,
    int pivot_col_limit,
    bool reduced
) {
    assert(0 <= pivot_col_limit && pivot_col_limit <= matrix.cols());
    std::vector<int> pivot_columns;
    int pivot_row = 0;
    for (
        int col = 0;
        col < pivot_col_limit && pivot_row < matrix.rows();
        col++
    ) {
        int pivot = -1;
        for (int row = pivot_row; row < matrix.rows(); row++) {
            if (matrix.get(row, col)) {
                pivot = row;
                break;
            }
        }
        if (pivot == -1) continue;
        matrix.swap_rows(pivot_row, pivot);

        const int first_row = reduced ? 0 : pivot_row + 1;
        for (int row = first_row; row < matrix.rows(); row++) {
            if (row != pivot_row && matrix.get(row, col)) {
                matrix.xor_rows(row, pivot_row, col);
            }
        }
        pivot_columns.push_back(col);
        pivot_row++;
    }
    return pivot_columns;
}

}  // namespace bit_matrix_detail

struct BitRowReduction {
    BitMatrix matrix;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }
};

inline BitRowReduction reduced_row_echelon_form(BitMatrix matrix) {
    BitRowReduction result;
    result.pivot_columns = bit_matrix_detail::row_reduce(
        matrix,
        matrix.cols(),
        true
    );
    result.matrix = std::move(matrix);
    return result;
}

inline int matrix_rank(BitMatrix matrix) {
    if (matrix.rows() > matrix.cols()) matrix = matrix.transposed();
    return int(bit_matrix_detail::row_reduce(
        matrix,
        matrix.cols(),
        false
    ).size());
}

inline bool determinant(const BitMatrix& matrix) {
    assert(matrix.rows() == matrix.cols());
    return matrix_rank(matrix) == matrix.rows();
}

inline std::optional<BitMatrix> inverse(const BitMatrix& matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    BitMatrix augmented(size, 2 * size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            if (matrix.get(row, col)) augmented.set(row, col);
        }
        augmented.set(row, size + row);
    }

    const std::vector<int> pivots = bit_matrix_detail::row_reduce(
        augmented,
        size,
        true
    );
    if (int(pivots.size()) != size) return std::nullopt;

    BitMatrix result(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            if (augmented.get(row, size + col)) result.set(row, col);
        }
    }
    return result;
}

struct BitLinearSystemResult {
    bool consistent = false;
    std::vector<bool> particular_solution;
    std::vector<std::vector<bool>> nullspace_basis;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }

    int nullity() const {
        return consistent ? int(nullspace_basis.size()) : 0;
    }

    bool has_unique_solution() const {
        return consistent && nullspace_basis.empty();
    }
};

inline BitLinearSystemResult solve_linear_system(
    const BitMatrix& coefficients,
    const std::vector<bool>& constants
) {
    assert(coefficients.rows() == int(constants.size()));
    const int equation_count = coefficients.rows();
    const int variable_count = coefficients.cols();
    BitMatrix augmented(equation_count, variable_count + 1);
    for (int row = 0; row < equation_count; row++) {
        for (int col = 0; col < variable_count; col++) {
            if (coefficients.get(row, col)) augmented.set(row, col);
        }
        if (constants[std::size_t(row)]) augmented.set(row, variable_count);
    }

    BitLinearSystemResult result;
    result.pivot_columns = bit_matrix_detail::row_reduce(
        augmented,
        variable_count,
        true
    );
    for (int row = result.rank(); row < equation_count; row++) {
        if (augmented.get(row, variable_count)) return result;
    }

    result.consistent = true;
    result.particular_solution.assign(std::size_t(variable_count), false);
    std::vector<bool> is_pivot(std::size_t(variable_count), false);
    for (int row = 0; row < result.rank(); row++) {
        const int col = result.pivot_columns[std::size_t(row)];
        is_pivot[std::size_t(col)] = true;
        result.particular_solution[std::size_t(col)] =
            augmented.get(row, variable_count);
    }

    for (int free_col = 0; free_col < variable_count; free_col++) {
        if (is_pivot[std::size_t(free_col)]) continue;
        std::vector<bool> direction(std::size_t(variable_count), false);
        direction[std::size_t(free_col)] = true;
        for (int row = 0; row < result.rank(); row++) {
            const int pivot_col = result.pivot_columns[std::size_t(row)];
            direction[std::size_t(pivot_col)] = augmented.get(row, free_col);
        }
        result.nullspace_basis.push_back(std::move(direction));
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/characteristic_polynomial.hpp"



#line 8 "math/matrix/characteristic_polynomial.hpp"

#line 10 "math/matrix/characteristic_polynomial.hpp"

namespace m1une {
namespace matrix {

template <class T>
std::vector<T> characteristic_polynomial(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();

    for (int col = 0; col + 2 < size; col++) {
        int pivot = col + 1;
        while (pivot < size && matrix[pivot][col] == T()) pivot++;
        if (pivot == size) continue;

        if (pivot != col + 1) {
            matrix.swap_rows(pivot, col + 1);
            for (int row = 0; row < size; row++) {
                std::swap(matrix[row][pivot], matrix[row][col + 1]);
            }
        }

        const T inverse_pivot = T(1) / matrix[col + 1][col];
        for (int row = col + 2; row < size; row++) {
            if (matrix[row][col] == T()) continue;
            const T factor = matrix[row][col] * inverse_pivot;
            for (int j = col; j < size; j++) {
                matrix[row][j] -= factor * matrix[col + 1][j];
            }
            for (int i = 0; i < size; i++) {
                matrix[i][col + 1] += factor * matrix[i][row];
            }
        }
    }

    std::vector<std::vector<T>> polynomial(std::size_t(size + 1));
    polynomial[0].assign(1, T(1));
    for (int leading_size = 1; leading_size <= size; leading_size++) {
        const int last = leading_size - 1;
        polynomial[std::size_t(leading_size)].assign(
            std::size_t(leading_size + 1),
            T()
        );
        const std::vector<T>& previous =
            polynomial[std::size_t(leading_size - 1)];
        std::vector<T>& current = polynomial[std::size_t(leading_size)];

        for (int degree = 0; degree < leading_size; degree++) {
            current[std::size_t(degree)] -=
                previous[std::size_t(degree)] * matrix[last][last];
            current[std::size_t(degree + 1)] +=
                previous[std::size_t(degree)];
        }

        T subdiagonal_product = T(1);
        for (int row = last - 1; row >= 0; row--) {
            subdiagonal_product *= matrix[row + 1][row];
            const T factor = subdiagonal_product * matrix[row][last];
            if (factor == T()) continue;
            for (int degree = 0; degree <= row; degree++) {
                current[std::size_t(degree)] -=
                    factor * polynomial[std::size_t(row)][std::size_t(degree)];
            }
        }
    }
    return polynomial[std::size_t(size)];
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/determinant_mod.hpp"



#line 6 "math/matrix/determinant_mod.hpp"
#include <type_traits>

#line 9 "math/matrix/determinant_mod.hpp"

namespace m1une {
namespace matrix {

namespace detail {

inline std::uint64_t determinant_multiply_mod(std::uint64_t lhs,
                                              std::uint64_t rhs,
                                              std::uint64_t modulus) {
    return std::uint64_t(static_cast<unsigned __int128>(lhs) * rhs % modulus);
}

inline std::uint64_t determinant_subtract_product_mod(
    std::uint64_t value, std::uint64_t lhs, std::uint64_t rhs,
    std::uint64_t modulus) {
    const std::uint64_t product = determinant_multiply_mod(lhs, rhs, modulus);
    return std::uint64_t((static_cast<unsigned __int128>(value) + modulus - product) %
                         modulus);
}

inline std::uint64_t determinant_add_products_mod(
    std::uint64_t first_lhs, std::uint64_t first_rhs,
    std::uint64_t second_lhs, std::uint64_t second_rhs,
    std::uint64_t modulus) {
    const std::uint64_t first =
        determinant_multiply_mod(first_lhs, first_rhs, modulus);
    const std::uint64_t second =
        determinant_multiply_mod(second_lhs, second_rhs, modulus);
    return std::uint64_t((static_cast<unsigned __int128>(first) + second) % modulus);
}

template <class Integer>
std::uint64_t determinant_normalize(Integer value, std::uint64_t modulus) {
    static_assert(std::is_integral_v<Integer>);
    static_assert(sizeof(Integer) <= sizeof(std::uint64_t));
    if constexpr (std::is_signed_v<Integer>) {
        __int128 residue = static_cast<__int128>(value) % static_cast<__int128>(modulus);
        if (residue < 0) residue += modulus;
        return std::uint64_t(residue);
    } else {
        return std::uint64_t(static_cast<unsigned __int128>(value) % modulus);
    }
}

}  // namespace detail

template <class Integer>
std::uint64_t determinant_mod(const Matrix<Integer>& matrix,
                              std::uint64_t modulus) {
    static_assert(std::is_integral_v<Integer>);
    assert(matrix.rows() == matrix.cols());
    assert(modulus > 0);
    const int size = matrix.rows();
    if (size == 0) return std::uint64_t(1) % modulus;

    Matrix<std::uint64_t> reduced(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            reduced[row][col] =
                detail::determinant_normalize(matrix[row][col], modulus);
        }
    }

    std::uint64_t result = std::uint64_t(1) % modulus;
    bool negate = false;
    for (int col = 0; col < size; col++) {
        int pivot = col;
        while (pivot < size && reduced[pivot][col] == 0) pivot++;
        if (pivot == size) return 0;
        if (pivot != col) {
            reduced.swap_rows(pivot, col);
            negate = !negate;
        }

        for (int row = col + 1; row < size; row++) {
            std::uint64_t upper = reduced[col][col];
            std::uint64_t lower = reduced[row][col];
            if (lower == 0) continue;

            std::uint64_t upper_upper = 1 % modulus;
            std::uint64_t upper_lower = 0;
            std::uint64_t lower_upper = 0;
            std::uint64_t lower_lower = 1 % modulus;
            while (upper != 0 && lower != 0) {
                if (upper < lower) {
                    const std::uint64_t quotient = lower / upper;
                    lower -= quotient * upper;
                    lower_upper = detail::determinant_subtract_product_mod(
                        lower_upper, quotient, upper_upper, modulus);
                    lower_lower = detail::determinant_subtract_product_mod(
                        lower_lower, quotient, upper_lower, modulus);
                } else {
                    const std::uint64_t quotient = upper / lower;
                    upper -= quotient * lower;
                    upper_upper = detail::determinant_subtract_product_mod(
                        upper_upper, quotient, lower_upper, modulus);
                    upper_lower = detail::determinant_subtract_product_mod(
                        upper_lower, quotient, lower_lower, modulus);
                }
            }

            for (int index = col; index < size; index++) {
                const std::uint64_t old_upper = reduced[col][index];
                const std::uint64_t old_lower = reduced[row][index];
                reduced[col][index] = detail::determinant_add_products_mod(
                    upper_upper, old_upper, upper_lower, old_lower, modulus);
                reduced[row][index] = detail::determinant_add_products_mod(
                    lower_upper, old_upper, lower_lower, old_lower, modulus);
            }
            if (upper == 0) {
                reduced.swap_rows(col, row);
                negate = !negate;
            }
        }

        result = detail::determinant_multiply_mod(
            result, reduced[col][col], modulus);
        if (result == 0) return 0;
    }
    return negate ? modulus - result : result;
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/hafnian.hpp"



#line 7 "math/matrix/hafnian.hpp"

#line 9 "math/matrix/hafnian.hpp"

namespace m1une {
namespace matrix {
namespace internal {

template <class T>
class HafnianSolver {
    using Polynomial = std::vector<T>;
    using PolynomialMatrix = std::vector<std::vector<Polynomial>>;

    int _degree;

    void add_shifted_product(Polynomial& result, const Polynomial& first,
                             const Polynomial& second) const {
        for (int first_degree = 0; first_degree < _degree; first_degree++) {
            for (int second_degree = 0;
                 first_degree + second_degree + 1 < _degree;
                 second_degree++) {
                result[first_degree + second_degree + 1] +=
                    first[first_degree] * second[second_degree];
            }
        }
    }

    Polynomial solve(PolynomialMatrix matrix) const {
        if (matrix.empty()) {
            Polynomial result(_degree);
            result[0] = T(1);
            return result;
        }

        std::vector<Polynomial> first = std::move(matrix.back());
        matrix.pop_back();
        std::vector<Polynomial> second = std::move(matrix.back());
        matrix.pop_back();
        const int remaining = int(matrix.size());
        Polynomial first_to_pair = std::move(first[remaining]);

        Polynomial result = solve(matrix);
        for (T& coefficient : result) coefficient = T() - coefficient;

        for (int row = 0; row < remaining; row++) {
            for (int col = 0; col < row; col++) {
                add_shifted_product(matrix[row][col], first[row], second[col]);
                add_shifted_product(matrix[row][col], second[row], first[col]);
            }
        }

        Polynomial with_connections = solve(std::move(matrix));
        add_shifted_product(result, first_to_pair, with_connections);
        for (int degree = 0; degree < _degree; degree++) {
            result[degree] += with_connections[degree];
        }
        return result;
    }

   public:
    explicit HafnianSolver(int size) : _degree(size / 2 + 1) {}

    T operator()(const Matrix<T>& matrix) const {
        const int size = matrix.rows();
        PolynomialMatrix polynomial_matrix(size);
        for (int row = 0; row < size; row++) {
            polynomial_matrix[row].assign(row, Polynomial(_degree));
            for (int col = 0; col < row; col++) {
                polynomial_matrix[row][col][0] = matrix[row][col];
            }
        }
        return solve(std::move(polynomial_matrix)).back();
    }
};

}  // namespace internal

// Returns the hafnian of an even-dimensional symmetric zero-diagonal matrix.
template <class T>
T hafnian(const Matrix<T>& matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    assert(size % 2 == 0);

#ifndef NDEBUG
    for (int row = 0; row < size; row++) {
        assert(matrix[row][row] == T());
        for (int col = row + 1; col < size; col++) {
            assert(matrix[row][col] == matrix[col][row]);
        }
    }
#endif

    return internal::HafnianSolver<T>(size)(matrix);
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/linear_algebra.hpp"



#line 7 "math/matrix/linear_algebra.hpp"

#line 9 "math/matrix/linear_algebra.hpp"

namespace m1une {
namespace matrix {

template <class T>
constexpr T default_epsilon() {
    if constexpr (std::is_floating_point_v<T>) {
        return T(1e-10);
    } else {
        return T();
    }
}

namespace detail {

template <class T>
T matrix_abs(T value) {
    return value < T() ? T() - value : value;
}

template <class T>
bool is_zero(const T& value, const T& eps) {
    if constexpr (std::is_floating_point_v<T>) {
        return matrix_abs(value) <= eps;
    } else {
        (void)eps;
        return value == T();
    }
}

template <class T>
int choose_pivot(const Matrix<T>& matrix, int first_row, int col, const T& eps) {
    int pivot = -1;
    if constexpr (std::is_floating_point_v<T>) {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (is_zero(matrix[row][col], eps)) continue;
            if (pivot == -1 || matrix_abs(matrix[pivot][col]) < matrix_abs(matrix[row][col])) {
                pivot = row;
            }
        }
    } else {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (!is_zero(matrix[row][col], eps)) {
                pivot = row;
                break;
            }
        }
    }
    return pivot;
}

template <class T>
std::vector<int> row_reduce(Matrix<T>& matrix, int pivot_col_limit, const T& eps,
                            bool reduced) {
    std::vector<int> pivot_columns;
    int pivot_row = 0;
    for (int col = 0; col < pivot_col_limit && pivot_row < matrix.rows(); col++) {
        int pivot = choose_pivot(matrix, pivot_row, col, eps);
        if (pivot == -1) continue;
        matrix.swap_rows(pivot_row, pivot);

        const T pivot_value = matrix[pivot_row][col];
        if (reduced) {
            for (int j = col; j < matrix.cols(); j++) matrix[pivot_row][j] /= pivot_value;
        }

        const int first_row = reduced ? 0 : pivot_row + 1;
        for (int row = first_row; row < matrix.rows(); row++) {
            if (row == pivot_row || is_zero(matrix[row][col], eps)) continue;
            T factor = matrix[row][col];
            if (!reduced) factor /= pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < matrix.cols(); j++) {
                matrix[row][j] -= factor * matrix[pivot_row][j];
            }
        }

        pivot_columns.push_back(col);
        pivot_row++;
    }

    if constexpr (std::is_floating_point_v<T>) {
        for (T& value : matrix.data()) {
            if (is_zero(value, eps)) value = T();
        }
    }
    return pivot_columns;
}

}  // namespace detail

template <class T>
struct RowReduction {
    Matrix<T> matrix;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }
};

template <class T>
RowReduction<T> reduced_row_echelon_form(Matrix<T> matrix,
                                         T eps = default_epsilon<T>()) {
    RowReduction<T> result;
    result.pivot_columns = detail::row_reduce(matrix, matrix.cols(), eps, true);
    result.matrix = std::move(matrix);
    return result;
}

template <class T>
int matrix_rank(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    return int(detail::row_reduce(matrix, matrix.cols(), eps, false).size());
}

template <class T>
T determinant(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    T result = T(1);
    bool negate = false;

    for (int col = 0; col < size; col++) {
        int pivot = detail::choose_pivot(matrix, col, col, eps);
        if (pivot == -1) return T();
        if (pivot != col) {
            matrix.swap_rows(pivot, col);
            negate = !negate;
        }

        const T pivot_value = matrix[col][col];
        result *= pivot_value;
        for (int row = col + 1; row < size; row++) {
            if (detail::is_zero(matrix[row][col], eps)) continue;
            const T factor = matrix[row][col] / pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < size; j++) {
                matrix[row][j] -= factor * matrix[col][j];
            }
        }
    }
    return negate ? T() - result : result;
}

template <class T>
std::optional<Matrix<T>> inverse(const Matrix<T>& matrix,
                                 T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    const std::vector<int> pivots = detail::row_reduce(augmented, size, eps, true);
    if (int(pivots.size()) != size) return std::nullopt;

    Matrix<T> result(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] = augmented[row][size + col];
        }
    }
    return result;
}

template <class T>
struct LinearSystemResult {
    bool consistent = false;
    std::vector<T> particular_solution;
    std::vector<std::vector<T>> nullspace_basis;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }

    int nullity() const {
        return consistent ? int(nullspace_basis.size()) : 0;
    }

    bool has_unique_solution() const {
        return consistent && nullspace_basis.empty();
    }
};

template <class T>
LinearSystemResult<T> solve_linear_system(const Matrix<T>& coefficients,
                                          const std::vector<T>& constants,
                                          T eps = default_epsilon<T>()) {
    assert(coefficients.rows() == int(constants.size()));
    const int equation_count = coefficients.rows();
    const int variable_count = coefficients.cols();
    Matrix<T> augmented(equation_count, variable_count + 1);
    for (int row = 0; row < equation_count; row++) {
        for (int col = 0; col < variable_count; col++) {
            augmented[row][col] = coefficients[row][col];
        }
        augmented[row][variable_count] = constants[std::size_t(row)];
    }

    LinearSystemResult<T> result;
    result.pivot_columns =
        detail::row_reduce(augmented, variable_count, eps, true);

    for (int row = result.rank(); row < equation_count; row++) {
        bool zero_left = true;
        for (int col = 0; col < variable_count; col++) {
            if (!detail::is_zero(augmented[row][col], eps)) {
                zero_left = false;
                break;
            }
        }
        if (zero_left && !detail::is_zero(augmented[row][variable_count], eps)) {
            return result;
        }
    }

    result.consistent = true;
    result.particular_solution.assign(std::size_t(variable_count), T());
    std::vector<bool> is_pivot(std::size_t(variable_count), false);
    for (int row = 0; row < result.rank(); row++) {
        const int col = result.pivot_columns[std::size_t(row)];
        is_pivot[std::size_t(col)] = true;
        result.particular_solution[std::size_t(col)] = augmented[row][variable_count];
    }

    for (int free_col = 0; free_col < variable_count; free_col++) {
        if (is_pivot[std::size_t(free_col)]) continue;
        std::vector<T> direction(static_cast<std::size_t>(variable_count));
        direction[std::size_t(free_col)] = T(1);
        for (int row = 0; row < result.rank(); row++) {
            const int pivot_col = result.pivot_columns[std::size_t(row)];
            direction[std::size_t(pivot_col)] = T() - augmented[row][free_col];
        }
        result.nullspace_basis.push_back(std::move(direction));
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/pfaffian.hpp"



#line 6 "math/matrix/pfaffian.hpp"

#line 8 "math/matrix/pfaffian.hpp"

namespace m1une {
namespace matrix {

// Returns the Pfaffian of an even-dimensional alternating matrix over a field.
template <class T>
T pfaffian(Matrix<T> matrix) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    assert(size % 2 == 0);

#ifndef NDEBUG
    for (int row = 0; row < size; row++) {
        assert(matrix[row][row] == T());
        for (int col = row + 1; col < size; col++) {
            assert(matrix[row][col] == T() - matrix[col][row]);
        }
    }
#endif

    T result = T(1);
    for (int first = 0; first < size; first += 2) {
        int pivot = first + 1;
        while (pivot < size && matrix[first][pivot] == T()) pivot++;
        if (pivot == size) return T();

        if (pivot != first + 1) {
            matrix.swap_rows(pivot, first + 1);
            for (int row = 0; row < size; row++) {
                std::swap(matrix[row][pivot], matrix[row][first + 1]);
            }
            result = T() - result;
        }

        const int second = first + 1;
        const T pivot_value = matrix[first][second];
        result *= pivot_value;
        const T inverse_pivot = T(1) / pivot_value;

        for (int row = second + 1; row < size; row++) {
            for (int col = row + 1; col < size; col++) {
                matrix[row][col] +=
                    (matrix[second][row] * matrix[first][col] -
                     matrix[first][row] * matrix[second][col]) *
                    inverse_pivot;
                matrix[col][row] = T() - matrix[row][col];
            }
        }
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une


#line 1 "math/matrix/sparse_determinant.hpp"



#line 8 "math/matrix/sparse_determinant.hpp"

namespace m1une {
namespace matrix {

template <class T>
struct SparseMatrixEntry {
    int row;
    int col;
    T value;
};

namespace internal {

struct SparseDeterminantRandom {
    std::uint64_t state;

    explicit SparseDeterminantRandom(std::uint64_t seed) : state(seed) {}

    std::uint64_t operator()() {
        std::uint64_t value = (state += 0x9e3779b97f4a7c15ULL);
        value = (value ^ (value >> 30)) * 0xbf58476d1ce4e5b9ULL;
        value = (value ^ (value >> 27)) * 0x94d049bb133111ebULL;
        return value ^ (value >> 31);
    }
};

template <class T>
std::vector<T> berlekamp_massey(const std::vector<T>& sequence) {
    std::vector<T> recurrence(1, T(1));
    std::vector<T> previous(1, T(1));
    int degree = 0;
    int shift = 1;
    T previous_discrepancy = T(1);

    for (int index = 0; index < int(sequence.size()); index++) {
        T discrepancy = sequence[index];
        for (int i = 1; i <= degree; i++) {
            discrepancy += recurrence[i] * sequence[index - i];
        }
        if (discrepancy == T()) {
            shift++;
            continue;
        }

        const T factor = discrepancy / previous_discrepancy;
        std::vector<T> old_recurrence = recurrence;
        if (int(recurrence.size()) < int(previous.size()) + shift) {
            recurrence.resize(previous.size() + std::size_t(shift), T());
        }
        for (int i = 0; i < int(previous.size()); i++) {
            recurrence[i + shift] -= factor * previous[i];
        }

        if (2 * degree <= index) {
            degree = index + 1 - degree;
            previous = std::move(old_recurrence);
            previous_discrepancy = discrepancy;
            shift = 1;
        } else {
            shift++;
        }
    }
    recurrence.resize(std::size_t(degree + 1));
    return recurrence;
}

}  // namespace internal

// Randomized black-box determinant over a finite field. random_nonzero must
// return independent nonzero field elements.
template <class T, class RandomValue>
T sparse_determinant_with_randomizer(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    RandomValue random_nonzero
) {
    assert(size >= 0);
    for (const SparseMatrixEntry<T>& entry : entries) {
        assert(0 <= entry.row && entry.row < size);
        assert(0 <= entry.col && entry.col < size);
    }
    if (size == 0) return T(1);

    auto random_vector = [&]() {
        std::vector<T> result(size);
        for (T& value : result) {
            value = random_nonzero();
            assert(value != T());
        }
        return result;
    };

    while (true) {
        std::vector<T> diagonal = random_vector();
        std::vector<T> left = random_vector();
        std::vector<T> state = random_vector();
        std::vector<T> sequence(std::size_t(2 * size));

        for (int step = 0; step < 2 * size; step++) {
            for (int i = 0; i < size; i++) sequence[step] += left[i] * state[i];
            for (int i = 0; i < size; i++) state[i] *= diagonal[i];

            std::vector<T> next(size);
            for (const SparseMatrixEntry<T>& entry : entries) {
                next[entry.row] += entry.value * state[entry.col];
            }
            state = std::move(next);
        }

        std::vector<T> recurrence = internal::berlekamp_massey(sequence);
        if (recurrence.back() == T()) return T();
        if (int(recurrence.size()) != size + 1) continue;

        T determinant = recurrence.back();
        if (size % 2 == 1) determinant = T() - determinant;
        for (const T& value : diagonal) determinant /= value;
        return determinant;
    }
}

template <class T>
T sparse_determinant(
    int size, const std::vector<SparseMatrixEntry<T>>& entries,
    std::uint64_t seed = 0x243f6a8885a308d3ULL
) {
    const std::uint64_t modulus = T::mod();
    assert(modulus > 1);
    internal::SparseDeterminantRandom random(seed);
    auto random_nonzero = [&]() {
        return T(1 + random() % (modulus - 1));
    };
    return sparse_determinant_with_randomizer<T>(size, entries, random_nonzero);
}

}  // namespace matrix
}  // namespace m1une


#line 13 "math/matrix/all.hpp"
Back to top page