f2_matrix.hpp¶
Return the rank of a packed binary matrix. Gaussian elimination uses XOR as row addition; when the matrix is tall, transposition first minimizes the number of packed rows participating in elimination.
Verified by inverse_matrix_mod_2, matrix_det_mod_2, matrix_product_mod_2, matrix_rank_mod_2, system_of_linear_equations_mod_2.
\[
\displaystyle \mathrm{rank}_{\mathbf F_2}(A)
\]
Implementation¶
#ifndef NOYA_F2_MATRIX_HPP
#define NOYA_F2_MATRIX_HPP 1
/// @complexity Time: O(rows * columns * min(rows,columns) / 64) for
/// elimination, O(rows * inner * columns / 64) for multiplication.
/// Space: O(rows * columns / 64).
#include "noya/dynamic_bitset.hpp"
#include <algorithm>
#include <cassert>
#include <optional>
#include <utility>
#include <vector>
namespace noya {
using f2_matrix = std::vector<dynamic_bitset>;
namespace f2_matrix_detail {
inline std::size_t columns(const f2_matrix &matrix) {
if (matrix.empty()) {
return 0;
}
std::size_t result = matrix.front().size();
for (const auto &row : matrix) {
assert(row.size() == result);
}
return result;
}
inline f2_matrix transpose(const f2_matrix &matrix, std::size_t columns) {
f2_matrix result(columns, dynamic_bitset(matrix.size()));
for (std::size_t row = 0; row < matrix.size(); row++) {
assert(matrix[row].size() == columns);
for (std::size_t column = matrix[row].find_first(); column < columns;
column = matrix[row].find_next(column + 1)) {
result[column].set(row);
}
}
return result;
}
} // namespace f2_matrix_detail
/// @brief Return the rank of a packed binary matrix. Gaussian elimination
/// uses XOR as row addition; when the matrix is tall, transposition first
/// minimizes the number of packed rows participating in elimination.
inline int f2_matrix_rank(f2_matrix matrix, std::size_t columns) {
for (const auto &row : matrix) {
assert(row.size() == columns);
}
if (matrix.size() > columns) {
matrix = f2_matrix_detail::transpose(matrix, columns);
columns = matrix.empty() ? 0 : matrix.front().size();
}
int rank = 0;
for (std::size_t column = 0;
column < columns && rank < int(matrix.size()); column++) {
int pivot = rank;
while (pivot < int(matrix.size()) && !matrix[pivot][column]) {
pivot++;
}
if (pivot == int(matrix.size())) {
continue;
}
std::swap(matrix[pivot], matrix[rank]);
for (int row = rank + 1; row < int(matrix.size()); row++) {
if (matrix[row][column]) {
matrix[row] ^= matrix[rank];
}
}
rank++;
}
return rank;
}
inline int f2_matrix_rank(f2_matrix matrix) {
std::size_t columns = f2_matrix_detail::columns(matrix);
return f2_matrix_rank(std::move(matrix), columns);
}
/// @brief Return the determinant of a square binary matrix. Over F2 row
/// swaps have no sign cost, so the determinant is one exactly when every
/// column obtains a pivot.
inline bool f2_determinant(f2_matrix matrix) {
int size = int(matrix.size());
assert(f2_matrix_detail::columns(matrix) == std::size_t(size));
for (int column = 0; column < size; column++) {
int pivot = column;
while (pivot < size && !matrix[pivot][column]) {
pivot++;
}
if (pivot == size) {
return false;
}
std::swap(matrix[pivot], matrix[column]);
for (int row = column + 1; row < size; row++) {
if (matrix[row][column]) {
matrix[row] ^= matrix[column];
}
}
}
return true;
}
/// @brief Multiply two packed binary matrices. A set entry in a left row
/// selects the corresponding right row, and XOR of all selected rows is the
/// output row.
inline f2_matrix f2_matrix_product(const f2_matrix &left,
const f2_matrix &right,
std::size_t inner_columns,
std::size_t output_columns) {
for (const auto &row : left) {
assert(row.size() == inner_columns);
}
assert(right.size() == inner_columns);
for (const auto &row : right) {
assert(row.size() == output_columns);
}
f2_matrix result(left.size(), dynamic_bitset(output_columns));
for (std::size_t row = 0; row < left.size(); row++) {
for (std::size_t index = left[row].find_first(); index < inner_columns;
index = left[row].find_next(index + 1)) {
result[row] ^= right[index];
}
}
return result;
}
/// @brief Return the inverse of a square binary matrix, or nullopt when it is
/// singular. Gauss--Jordan elimination applies every row operation to an
/// identity matrix in parallel, leaving the inverse after the left side is
/// reduced to identity.
inline std::optional<f2_matrix> f2_matrix_inverse(f2_matrix matrix) {
int size = int(matrix.size());
assert(f2_matrix_detail::columns(matrix) == std::size_t(size));
f2_matrix inverse(size, dynamic_bitset(size));
for (int row = 0; row < size; row++) {
inverse[row].set(row);
}
for (int column = 0; column < size; column++) {
int pivot = column;
while (pivot < size && !matrix[pivot][column]) {
pivot++;
}
if (pivot == size) {
return std::nullopt;
}
std::swap(matrix[pivot], matrix[column]);
std::swap(inverse[pivot], inverse[column]);
for (int row = 0; row < size; row++) {
if (row != column && matrix[row][column]) {
matrix[row] ^= matrix[column];
inverse[row] ^= inverse[column];
}
}
}
return inverse;
}
struct f2_linear_system_solution {
bool consistent = false;
dynamic_bitset solution;
f2_matrix nullspace_basis;
std::vector<int> pivot_columns;
};
/// @brief Solve `matrix * x = right_hand_side` over F2. Reduced row echelon
/// form yields one solution by setting all free variables to zero; setting one
/// free variable at a time then gives a basis of the homogeneous nullspace.
inline f2_linear_system_solution
solve_f2_linear_system(f2_matrix matrix, std::size_t columns,
const dynamic_bitset &right_hand_side) {
int rows = int(matrix.size());
assert(right_hand_side.size() == std::size_t(rows));
for (const auto &row : matrix) {
assert(row.size() == columns);
}
std::vector<bool> right(rows);
for (int row = 0; row < rows; row++) {
right[row] = right_hand_side[row];
}
std::vector<int> pivots;
int rank = 0;
for (std::size_t column = 0; column < columns && rank < rows; column++) {
int pivot = rank;
while (pivot < rows && !matrix[pivot][column]) {
pivot++;
}
if (pivot == rows) {
continue;
}
std::swap(matrix[pivot], matrix[rank]);
bool temporary = right[pivot];
right[pivot] = right[rank];
right[rank] = temporary;
for (int row = 0; row < rows; row++) {
if (row != rank && matrix[row][column]) {
matrix[row] ^= matrix[rank];
right[row] = right[row] != right[rank];
}
}
pivots.push_back(int(column));
rank++;
}
for (int row = rank; row < rows; row++) {
if (right[row]) {
return {false, dynamic_bitset(columns), {}, std::move(pivots)};
}
}
f2_linear_system_solution result;
result.consistent = true;
result.solution = dynamic_bitset(columns);
result.pivot_columns = pivots;
std::vector<int> pivot_row(columns, -1);
for (int row = 0; row < rank; row++) {
int column = pivots[row];
pivot_row[column] = row;
result.solution.set(column, right[row]);
}
for (std::size_t free_column = 0; free_column < columns; free_column++) {
if (pivot_row[free_column] != -1) {
continue;
}
dynamic_bitset basis(columns);
basis.set(free_column);
for (int row = 0; row < rank; row++) {
if (matrix[row][free_column]) {
basis.set(pivots[row]);
}
}
result.nullspace_basis.push_back(std::move(basis));
}
return result;
}
} // namespace noya
#endif // NOYA_F2_MATRIX_HPP
#include <algorithm>
#include <bit>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <optional>
#include <utility>
#include <vector>
/// @complexity Time: O(rows * columns * min(rows,columns) / 64) for
/// elimination, O(rows * inner * columns / 64) for multiplication.
/// Space: O(rows * columns / 64).
/// @complexity Time: O(n / 64) for whole-bitset operations; O(1) bit access.
/// Space: O(n / 64).
namespace noya {
/// @brief Resizable packed bitset with bitwise operations, shifts, population
/// count, and efficient iteration over set bits.
struct dynamic_bitset {
using word_type = std::uint64_t;
static constexpr std::size_t bits_per_word = 64;
std::size_t bit_size = 0;
std::vector<word_type> words;
dynamic_bitset() = default;
explicit dynamic_bitset(std::size_t size, bool value = false)
: bit_size(size), words(word_count(size), value ? ~word_type{} : 0) {
trim();
}
std::size_t size() const { return bit_size; }
bool empty() const { return bit_size == 0; }
bool test(std::size_t position) const {
assert(position < bit_size);
return (words[position / bits_per_word] >> (position % bits_per_word)) & 1;
}
bool operator[](std::size_t position) const { return test(position); }
dynamic_bitset &set(std::size_t position, bool value = true) {
assert(position < bit_size);
word_type mask = word_type(1) << (position % bits_per_word);
if (value) {
words[position / bits_per_word] |= mask;
} else {
words[position / bits_per_word] &= ~mask;
}
return *this;
}
dynamic_bitset &reset(std::size_t position) { return set(position, false); }
dynamic_bitset &flip(std::size_t position) {
assert(position < bit_size);
words[position / bits_per_word] ^=
word_type(1) << (position % bits_per_word);
return *this;
}
dynamic_bitset &set() {
std::fill(words.begin(), words.end(), ~word_type{});
trim();
return *this;
}
dynamic_bitset &reset() {
std::fill(words.begin(), words.end(), word_type{});
return *this;
}
dynamic_bitset &flip() {
for (word_type &word : words) {
word = ~word;
}
trim();
return *this;
}
std::size_t count() const {
std::size_t result = 0;
for (word_type word : words) {
result += std::popcount(word);
}
return result;
}
bool any() const {
return std::any_of(words.begin(), words.end(),
[](word_type word) { return word != 0; });
}
bool none() const { return !any(); }
/// @brief Return the first set position at least position, or size() if no
/// such position exists.
std::size_t find_next(std::size_t position) const {
if (position >= bit_size) {
return bit_size;
}
std::size_t index = position / bits_per_word;
word_type word =
words[index] & (~word_type{} << (position % bits_per_word));
if (word != 0) {
return std::min(bit_size,
index * bits_per_word + std::size_t(std::countr_zero(word)));
}
for (index++; index < words.size(); index++) {
if (words[index] != 0) {
return std::min(
bit_size, index * bits_per_word +
std::size_t(std::countr_zero(words[index])));
}
}
return bit_size;
}
std::size_t find_first() const { return find_next(0); }
dynamic_bitset &operator&=(const dynamic_bitset &other) {
check_same_size(other);
for (std::size_t i = 0; i < words.size(); i++) {
words[i] &= other.words[i];
}
return *this;
}
dynamic_bitset &operator|=(const dynamic_bitset &other) {
check_same_size(other);
for (std::size_t i = 0; i < words.size(); i++) {
words[i] |= other.words[i];
}
return *this;
}
dynamic_bitset &operator^=(const dynamic_bitset &other) {
check_same_size(other);
for (std::size_t i = 0; i < words.size(); i++) {
words[i] ^= other.words[i];
}
return *this;
}
dynamic_bitset &operator<<=(std::size_t shift) {
if (shift >= bit_size) {
return reset();
}
std::size_t whole = shift / bits_per_word;
int part = int(shift % bits_per_word);
for (std::size_t i = words.size(); i-- > 0;) {
word_type value = 0;
if (i >= whole) {
value = words[i - whole] << part;
if (part != 0 && i > whole) {
value |= words[i - whole - 1] >> (bits_per_word - part);
}
}
words[i] = value;
}
trim();
return *this;
}
dynamic_bitset &operator>>=(std::size_t shift) {
if (shift >= bit_size) {
return reset();
}
std::size_t whole = shift / bits_per_word;
int part = int(shift % bits_per_word);
for (std::size_t i = 0; i < words.size(); i++) {
word_type value = 0;
if (i + whole < words.size()) {
value = words[i + whole] >> part;
if (part != 0 && i + whole + 1 < words.size()) {
value |= words[i + whole + 1] << (bits_per_word - part);
}
}
words[i] = value;
}
return *this;
}
friend dynamic_bitset operator&(dynamic_bitset first,
const dynamic_bitset &second) {
return first &= second;
}
friend dynamic_bitset operator|(dynamic_bitset first,
const dynamic_bitset &second) {
return first |= second;
}
friend dynamic_bitset operator^(dynamic_bitset first,
const dynamic_bitset &second) {
return first ^= second;
}
friend dynamic_bitset operator<<(dynamic_bitset value, std::size_t shift) {
return value <<= shift;
}
friend dynamic_bitset operator>>(dynamic_bitset value, std::size_t shift) {
return value >>= shift;
}
friend dynamic_bitset operator~(dynamic_bitset value) { return value.flip(); }
friend bool operator==(const dynamic_bitset &, const dynamic_bitset &) =
default;
private:
static std::size_t word_count(std::size_t size) {
return (size + bits_per_word - 1) / bits_per_word;
}
void trim() {
if (!words.empty() && bit_size % bits_per_word != 0) {
words.back() &=
(word_type(1) << (bit_size % bits_per_word)) - word_type(1);
}
}
void check_same_size(const dynamic_bitset &other) const {
assert(bit_size == other.bit_size);
}
};
} // namespace noya
namespace noya {
using f2_matrix = std::vector<dynamic_bitset>;
namespace f2_matrix_detail {
inline std::size_t columns(const f2_matrix &matrix) {
if (matrix.empty()) {
return 0;
}
std::size_t result = matrix.front().size();
for (const auto &row : matrix) {
assert(row.size() == result);
}
return result;
}
inline f2_matrix transpose(const f2_matrix &matrix, std::size_t columns) {
f2_matrix result(columns, dynamic_bitset(matrix.size()));
for (std::size_t row = 0; row < matrix.size(); row++) {
assert(matrix[row].size() == columns);
for (std::size_t column = matrix[row].find_first(); column < columns;
column = matrix[row].find_next(column + 1)) {
result[column].set(row);
}
}
return result;
}
} // namespace f2_matrix_detail
/// @brief Return the rank of a packed binary matrix. Gaussian elimination
/// uses XOR as row addition; when the matrix is tall, transposition first
/// minimizes the number of packed rows participating in elimination.
inline int f2_matrix_rank(f2_matrix matrix, std::size_t columns) {
for (const auto &row : matrix) {
assert(row.size() == columns);
}
if (matrix.size() > columns) {
matrix = f2_matrix_detail::transpose(matrix, columns);
columns = matrix.empty() ? 0 : matrix.front().size();
}
int rank = 0;
for (std::size_t column = 0;
column < columns && rank < int(matrix.size()); column++) {
int pivot = rank;
while (pivot < int(matrix.size()) && !matrix[pivot][column]) {
pivot++;
}
if (pivot == int(matrix.size())) {
continue;
}
std::swap(matrix[pivot], matrix[rank]);
for (int row = rank + 1; row < int(matrix.size()); row++) {
if (matrix[row][column]) {
matrix[row] ^= matrix[rank];
}
}
rank++;
}
return rank;
}
inline int f2_matrix_rank(f2_matrix matrix) {
std::size_t columns = f2_matrix_detail::columns(matrix);
return f2_matrix_rank(std::move(matrix), columns);
}
/// @brief Return the determinant of a square binary matrix. Over F2 row
/// swaps have no sign cost, so the determinant is one exactly when every
/// column obtains a pivot.
inline bool f2_determinant(f2_matrix matrix) {
int size = int(matrix.size());
assert(f2_matrix_detail::columns(matrix) == std::size_t(size));
for (int column = 0; column < size; column++) {
int pivot = column;
while (pivot < size && !matrix[pivot][column]) {
pivot++;
}
if (pivot == size) {
return false;
}
std::swap(matrix[pivot], matrix[column]);
for (int row = column + 1; row < size; row++) {
if (matrix[row][column]) {
matrix[row] ^= matrix[column];
}
}
}
return true;
}
/// @brief Multiply two packed binary matrices. A set entry in a left row
/// selects the corresponding right row, and XOR of all selected rows is the
/// output row.
inline f2_matrix f2_matrix_product(const f2_matrix &left,
const f2_matrix &right,
std::size_t inner_columns,
std::size_t output_columns) {
for (const auto &row : left) {
assert(row.size() == inner_columns);
}
assert(right.size() == inner_columns);
for (const auto &row : right) {
assert(row.size() == output_columns);
}
f2_matrix result(left.size(), dynamic_bitset(output_columns));
for (std::size_t row = 0; row < left.size(); row++) {
for (std::size_t index = left[row].find_first(); index < inner_columns;
index = left[row].find_next(index + 1)) {
result[row] ^= right[index];
}
}
return result;
}
/// @brief Return the inverse of a square binary matrix, or nullopt when it is
/// singular. Gauss--Jordan elimination applies every row operation to an
/// identity matrix in parallel, leaving the inverse after the left side is
/// reduced to identity.
inline std::optional<f2_matrix> f2_matrix_inverse(f2_matrix matrix) {
int size = int(matrix.size());
assert(f2_matrix_detail::columns(matrix) == std::size_t(size));
f2_matrix inverse(size, dynamic_bitset(size));
for (int row = 0; row < size; row++) {
inverse[row].set(row);
}
for (int column = 0; column < size; column++) {
int pivot = column;
while (pivot < size && !matrix[pivot][column]) {
pivot++;
}
if (pivot == size) {
return std::nullopt;
}
std::swap(matrix[pivot], matrix[column]);
std::swap(inverse[pivot], inverse[column]);
for (int row = 0; row < size; row++) {
if (row != column && matrix[row][column]) {
matrix[row] ^= matrix[column];
inverse[row] ^= inverse[column];
}
}
}
return inverse;
}
struct f2_linear_system_solution {
bool consistent = false;
dynamic_bitset solution;
f2_matrix nullspace_basis;
std::vector<int> pivot_columns;
};
/// @brief Solve `matrix * x = right_hand_side` over F2. Reduced row echelon
/// form yields one solution by setting all free variables to zero; setting one
/// free variable at a time then gives a basis of the homogeneous nullspace.
inline f2_linear_system_solution
solve_f2_linear_system(f2_matrix matrix, std::size_t columns,
const dynamic_bitset &right_hand_side) {
int rows = int(matrix.size());
assert(right_hand_side.size() == std::size_t(rows));
for (const auto &row : matrix) {
assert(row.size() == columns);
}
std::vector<bool> right(rows);
for (int row = 0; row < rows; row++) {
right[row] = right_hand_side[row];
}
std::vector<int> pivots;
int rank = 0;
for (std::size_t column = 0; column < columns && rank < rows; column++) {
int pivot = rank;
while (pivot < rows && !matrix[pivot][column]) {
pivot++;
}
if (pivot == rows) {
continue;
}
std::swap(matrix[pivot], matrix[rank]);
bool temporary = right[pivot];
right[pivot] = right[rank];
right[rank] = temporary;
for (int row = 0; row < rows; row++) {
if (row != rank && matrix[row][column]) {
matrix[row] ^= matrix[rank];
right[row] = right[row] != right[rank];
}
}
pivots.push_back(int(column));
rank++;
}
for (int row = rank; row < rows; row++) {
if (right[row]) {
return {false, dynamic_bitset(columns), {}, std::move(pivots)};
}
}
f2_linear_system_solution result;
result.consistent = true;
result.solution = dynamic_bitset(columns);
result.pivot_columns = pivots;
std::vector<int> pivot_row(columns, -1);
for (int row = 0; row < rank; row++) {
int column = pivots[row];
pivot_row[column] = row;
result.solution.set(column, right[row]);
}
for (std::size_t free_column = 0; free_column < columns; free_column++) {
if (pivot_row[free_column] != -1) {
continue;
}
dynamic_bitset basis(columns);
basis.set(free_column);
for (int row = 0; row < rank; row++) {
if (matrix[row][free_column]) {
basis.set(pivots[row]);
}
}
result.nullspace_basis.push_back(std::move(basis));
}
return result;
}
} // namespace noya