NicheLibrary
This documentation is automatically generated by NotLeonian/competitive-verifier (forked from competitive-verifier/competitive-verifier)
verify/standalone-dynamic-matrix-rank.test.cpp
- View this file on GitHub
- Last update: 2026-05-29 20:48:10+09:00
Depends on
Code
// competitive-verifier: STANDALONE
#include <cassert>
#include <vector>
#include "../math/matrix/dynamic-matrix-rank.hpp"
namespace {
struct F2 {
int value;
F2(int value = 0) : value(value & 1) {}
F2 &operator+=(const F2 &other) {
value ^= other.value;
return *this;
}
F2 &operator-=(const F2 &other) {
value ^= other.value;
return *this;
}
F2 &operator*=(const F2 &other) {
value &= other.value;
return *this;
}
F2 &operator/=(const F2 &other) {
assert(other.value == 1);
return *this;
}
[[maybe_unused]] friend F2 operator+(F2 lhs, const F2 &rhs) {
return lhs += rhs;
}
[[maybe_unused]] friend F2 operator-(F2 lhs, const F2 &rhs) {
return lhs -= rhs;
}
friend F2 operator*(F2 lhs, const F2 &rhs) { return lhs *= rhs; }
friend F2 operator/(F2 lhs, const F2 &rhs) { return lhs /= rhs; }
friend bool operator==(const F2 &lhs, const F2 &rhs) {
return lhs.value == rhs.value;
}
friend bool operator!=(const F2 &lhs, const F2 &rhs) {
return lhs.value != rhs.value;
}
};
template <class T> int brute_rank(std::vector<std::vector<T>> matrix) {
const int h = static_cast<int>(matrix.size());
if (h == 0) {
return 0;
}
const int w = static_cast<int>(matrix[0].size());
int rank = 0;
for (int column = 0; column < w; ++column) {
int pivot = -1;
for (int row = rank; row < h; ++row) {
if (matrix[row][column] != T()) {
pivot = row;
break;
}
}
if (pivot < 0) {
continue;
}
if (pivot != rank) {
matrix[pivot].swap(matrix[rank]);
}
const T inverse = T(1) / matrix[rank][column];
for (int j = column; j < w; ++j) {
matrix[rank][j] *= inverse;
}
for (int row = rank + 1; row < h; ++row) {
const T factor = matrix[row][column];
if (factor == T()) {
continue;
}
for (int j = column; j < w; ++j) {
matrix[row][j] -= matrix[rank][j] * factor;
}
}
++rank;
if (rank == h) {
break;
}
}
return rank;
}
std::vector<std::vector<F2>> matrix_from_mask(int h, int w, int mask) {
std::vector<std::vector<F2>> matrix(h, std::vector<F2>(w, F2()));
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
matrix[i][j] = F2((mask >> (i * w + j)) & 1);
}
}
return matrix;
}
std::vector<F2> vector_from_mask(int n, int mask) {
std::vector<F2> values(n, F2());
for (int i = 0; i < n; ++i) {
values[i] = F2((mask >> i) & 1);
}
return values;
}
std::vector<std::vector<F2>>
apply_rank_one(std::vector<std::vector<F2>> matrix,
const std::vector<F2> &column_vector,
const std::vector<F2> &row_vector) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
matrix[i][j] += column_vector[i] * row_vector[j];
}
}
return matrix;
}
void check_solver(const DynamicMatrixRank<F2> &solver,
const std::vector<std::vector<F2>> &matrix) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
assert(solver.rank() == brute_rank(matrix));
assert(solver.materialize_matrix() == matrix);
for (int row = 0; row < h; ++row) {
assert(solver.get_row(row) == matrix[row]);
}
for (int column = 0; column < w; ++column) {
std::vector<F2> expected(h, F2());
for (int row = 0; row < h; ++row) {
expected[row] = matrix[row][column];
}
assert(solver.get_column(column) == expected);
}
}
void check_updates(const std::vector<std::vector<F2>> &matrix) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
DynamicMatrixRank<F2> solver(matrix);
check_solver(solver, matrix);
for (int column_mask = 0; column_mask < (1 << h); ++column_mask) {
const std::vector<F2> column_vector = vector_from_mask(h, column_mask);
for (int row_mask = 0; row_mask < (1 << w); ++row_mask) {
const std::vector<F2> row_vector = vector_from_mask(w, row_mask);
const auto updated =
apply_rank_one(matrix, column_vector, row_vector);
const int expected = brute_rank(updated);
assert(solver.rank_after_rank_one_update(column_vector,
row_vector) == expected);
auto applied = solver;
assert(applied.apply_rank_one_update(column_vector, row_vector) ==
expected);
check_solver(applied, updated);
applied.build();
check_solver(applied, updated);
for (int row = 0; row < h; ++row) {
for (int new_row_mask = 0; new_row_mask < (1 << w);
++new_row_mask) {
const std::vector<F2> new_row =
vector_from_mask(w, new_row_mask);
auto row_updated = updated;
row_updated[row] = new_row;
const int row_expected = brute_rank(row_updated);
assert(applied.rank_after_row_replacement(row, new_row) ==
row_expected);
auto row_applied = applied;
assert(row_applied.apply_row_replacement(row, new_row) ==
row_expected);
check_solver(row_applied, row_updated);
}
}
for (int column = 0; column < w; ++column) {
for (int new_column_mask = 0; new_column_mask < (1 << h);
++new_column_mask) {
const std::vector<F2> new_column =
vector_from_mask(h, new_column_mask);
auto column_updated = updated;
for (int row = 0; row < h; ++row) {
column_updated[row][column] = new_column[row];
}
const int column_expected = brute_rank(column_updated);
assert(applied.rank_after_column_replacement(
column, new_column) == column_expected);
auto column_applied = applied;
assert(column_applied.apply_column_replacement(
column, new_column) == column_expected);
check_solver(column_applied, column_updated);
}
}
}
}
}
void self_test() {
{
DynamicMatrixRank<F2> solver;
solver.build(std::vector<std::vector<F2>>{});
assert(solver.rank() == 0);
assert(solver.materialize_matrix().empty());
}
for (int h = 1; h <= 3; ++h) {
for (int w = 1; w <= 3; ++w) {
const int states = 1 << (h * w);
for (int mask = 0; mask < states; ++mask) {
check_updates(matrix_from_mask(h, w, mask));
}
}
}
}
} // namespace
int main() {
self_test();
return 0;
}
#line 1 "verify/standalone-dynamic-matrix-rank.test.cpp"
// competitive-verifier: STANDALONE
#include <cassert>
#include <vector>
#line 1 "math/matrix/dynamic-matrix-rank.hpp"
// 行列の階数を求め、外積 1 項更新後の階数を判定する。
// さらに、内部状態を O(k(r + c)) で更新しつつ変更後の階数を返せる。
// 現在の行列は左右の階数分解と片側逆元で保持する。
// 1 行差し替え・1 列差し替えはそれぞれ e_i (b-a_i)^T, (b-a_j) e_j^T に帰着する。
// 前処理は O(rc min(r, c) + k^2(r + c))、各更新・各判定は O(k(r + c)) である。
// T は体を成し、零判定と四則演算ができることを仮定する。
#line 13 "math/matrix/dynamic-matrix-rank.hpp"
template <class T> struct DynamicMatrixRank {
int row_size = 0;
int column_size = 0;
int matrix_rank = 0;
std::vector<std::vector<T>> column_space_basis;
std::vector<std::vector<T>> row_space_basis;
std::vector<std::vector<T>> column_space_left_inverse;
std::vector<std::vector<T>> row_space_right_inverse;
DynamicMatrixRank() = default;
explicit DynamicMatrixRank(const std::vector<std::vector<T>> &matrix) {
build(matrix);
}
void build(const std::vector<std::vector<T>> &matrix) {
row_size = static_cast<int>(matrix.size());
column_size = row_size == 0 ? 0 : static_cast<int>(matrix[0].size());
for (int i = 1; i < row_size; ++i) {
assert(static_cast<int>(matrix[i].size()) == column_size);
}
const std::vector<int> basis_columns = find_independent_columns(matrix);
const std::vector<int> basis_rows =
find_independent_columns(transpose_matrix(matrix));
matrix_rank = static_cast<int>(basis_columns.size());
assert(static_cast<int>(basis_rows.size()) == matrix_rank);
std::vector<std::vector<T>> intersection_matrix(
matrix_rank, std::vector<T>(matrix_rank, T()));
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
intersection_matrix[i][j] =
matrix[basis_rows[i]][basis_columns[j]];
}
}
const std::vector<std::vector<T>> intersection_inverse =
inverse_matrix(intersection_matrix);
column_space_basis.assign(row_size, std::vector<T>(matrix_rank, T()));
for (int i = 0; i < row_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
column_space_basis[i][j] = matrix[i][basis_columns[j]];
}
}
row_space_basis.assign(matrix_rank, std::vector<T>(column_size, T()));
for (int i = 0; i < matrix_rank; ++i) {
for (int mid = 0; mid < matrix_rank; ++mid) {
const T value = intersection_inverse[i][mid];
if (value == T()) {
continue;
}
for (int j = 0; j < column_size; ++j) {
row_space_basis[i][j] += value * matrix[basis_rows[mid]][j];
}
}
}
column_space_left_inverse.assign(matrix_rank,
std::vector<T>(row_size, T()));
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
column_space_left_inverse[i][basis_rows[j]] =
intersection_inverse[i][j];
}
}
row_space_right_inverse.assign(column_size,
std::vector<T>(matrix_rank, T()));
for (int i = 0; i < matrix_rank; ++i) {
row_space_right_inverse[basis_columns[i]][i] = T(1);
}
}
void build() { build(materialize_matrix()); }
int rank() const { return matrix_rank; }
std::vector<T> get_row(int row_index) const {
assert(0 <= row_index && row_index < row_size);
std::vector<T> row(column_size, T());
for (int i = 0; i < matrix_rank; ++i) {
const T value = column_space_basis[row_index][i];
if (value == T()) {
continue;
}
for (int j = 0; j < column_size; ++j) {
row[j] += value * row_space_basis[i][j];
}
}
return row;
}
std::vector<T> get_column(int column_index) const {
assert(0 <= column_index && column_index < column_size);
std::vector<T> column(row_size, T());
for (int i = 0; i < row_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
column[i] +=
column_space_basis[i][j] * row_space_basis[j][column_index];
}
}
return column;
}
std::vector<std::vector<T>> materialize_matrix() const {
std::vector<std::vector<T>> matrix(row_size,
std::vector<T>(column_size, T()));
for (int i = 0; i < row_size; ++i) {
for (int mid = 0; mid < matrix_rank; ++mid) {
const T value = column_space_basis[i][mid];
if (value == T()) {
continue;
}
for (int j = 0; j < column_size; ++j) {
matrix[i][j] += value * row_space_basis[mid][j];
}
}
}
return matrix;
}
int rank_after_rank_one_update(const std::vector<T> &column_vector,
const std::vector<T> &row_vector) const {
return analyze_rank_one_update(column_vector, row_vector).next_rank;
}
int rank_after_row_replacement(int row_index,
const std::vector<T> &new_row) const {
assert(0 <= row_index && row_index < row_size);
assert(static_cast<int>(new_row.size()) == column_size);
std::vector<T> difference = new_row;
const std::vector<T> current_row = get_row(row_index);
for (int j = 0; j < column_size; ++j) {
difference[j] -= current_row[j];
}
std::vector<T> column_vector(row_size, T());
column_vector[row_index] = T(1);
return rank_after_rank_one_update(column_vector, difference);
}
int rank_after_column_replacement(int column_index,
const std::vector<T> &new_column) const {
assert(0 <= column_index && column_index < column_size);
assert(static_cast<int>(new_column.size()) == row_size);
std::vector<T> difference = new_column;
const std::vector<T> current_column = get_column(column_index);
for (int i = 0; i < row_size; ++i) {
difference[i] -= current_column[i];
}
std::vector<T> row_vector(column_size, T());
row_vector[column_index] = T(1);
return rank_after_rank_one_update(difference, row_vector);
}
int apply_rank_one_update(const std::vector<T> &column_vector,
const std::vector<T> &row_vector) {
RankOneUpdateInfo info =
analyze_rank_one_update(column_vector, row_vector);
if (info.next_rank == matrix_rank + 1) {
const int pivot_row = first_nonzero(info.column_residual);
const int pivot_column = first_nonzero(info.row_residual);
const std::vector<T> lambda = normalized_left_annihilator(
pivot_row, info.column_residual[pivot_row]);
const std::vector<T> rho = normalized_right_annihilator(
pivot_column, info.row_residual[pivot_column]);
const std::vector<std::vector<T>> old_row_space_right_inverse =
row_space_right_inverse;
append_column(column_space_basis, info.column_residual);
column_space_left_inverse.push_back(lambda);
row_space_basis.push_back(row_vector);
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < column_size; ++j) {
row_space_basis[i][j] += info.alpha[i] * row_vector[j];
}
}
std::vector<T> right_alpha(column_size, T());
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
right_alpha[i] +=
old_row_space_right_inverse[i][j] * info.alpha[j];
}
}
row_space_right_inverse.assign(
column_size, std::vector<T>(matrix_rank + 1, T()));
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
row_space_right_inverse[i][j] =
old_row_space_right_inverse[i][j] -
rho[i] * info.beta[j];
}
row_space_right_inverse[i][matrix_rank] =
T() - right_alpha[i] + rho[i] * info.schur;
}
++matrix_rank;
return matrix_rank;
}
if (!info.column_inside && info.row_inside) {
const int pivot_row = first_nonzero(info.column_residual);
const std::vector<T> lambda = normalized_left_annihilator(
pivot_row, info.column_residual[pivot_row]);
for (int i = 0; i < row_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
column_space_basis[i][j] += column_vector[i] * info.beta[j];
}
}
std::vector<std::vector<T>> new_left_inverse =
column_space_left_inverse;
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < row_size; ++j) {
new_left_inverse[i][j] -= info.alpha[i] * lambda[j];
}
}
column_space_left_inverse.swap(new_left_inverse);
return matrix_rank;
}
if (info.column_inside && !info.row_inside) {
const int pivot_column = first_nonzero(info.row_residual);
const std::vector<T> rho = normalized_right_annihilator(
pivot_column, info.row_residual[pivot_column]);
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < column_size; ++j) {
row_space_basis[i][j] += info.alpha[i] * row_vector[j];
}
}
std::vector<std::vector<T>> new_right_inverse =
row_space_right_inverse;
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
new_right_inverse[i][j] -= rho[i] * info.beta[j];
}
}
row_space_right_inverse.swap(new_right_inverse);
return matrix_rank;
}
if (info.schur != T()) {
std::vector<T> right_alpha(column_size, T());
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
right_alpha[i] +=
row_space_right_inverse[i][j] * info.alpha[j];
}
}
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < column_size; ++j) {
row_space_basis[i][j] += info.alpha[i] * row_vector[j];
}
}
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
row_space_right_inverse[i][j] -=
right_alpha[i] * (info.beta[j] / info.schur);
}
}
return matrix_rank;
}
const int removed = first_nonzero(info.alpha);
const T removed_inverse = T(1) / info.alpha[removed];
std::vector<T> right_alpha(column_size, T());
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
right_alpha[i] += row_space_right_inverse[i][j] * info.alpha[j];
}
}
std::vector<std::vector<T>> new_column_space_basis(
row_size, std::vector<T>(matrix_rank - 1, T()));
std::vector<std::vector<T>> new_row_space_basis(
matrix_rank - 1, std::vector<T>(column_size, T()));
std::vector<std::vector<T>> new_column_space_left_inverse(
matrix_rank - 1, std::vector<T>(row_size, T()));
std::vector<std::vector<T>> new_row_space_right_inverse(
column_size, std::vector<T>(matrix_rank - 1, T()));
int new_index = 0;
for (int old = 0; old < matrix_rank; ++old) {
if (old == removed) {
continue;
}
for (int i = 0; i < row_size; ++i) {
new_column_space_basis[i][new_index] =
column_space_basis[i][old] +
column_vector[i] * info.beta[old];
}
const T factor = info.alpha[old] * removed_inverse;
for (int j = 0; j < column_size; ++j) {
new_row_space_basis[new_index][j] =
row_space_basis[old][j] -
factor * row_space_basis[removed][j];
}
for (int j = 0; j < row_size; ++j) {
new_column_space_left_inverse[new_index][j] =
column_space_left_inverse[old][j] -
factor * column_space_left_inverse[removed][j];
}
for (int i = 0; i < column_size; ++i) {
new_row_space_right_inverse[i][new_index] =
row_space_right_inverse[i][old] +
right_alpha[i] * info.beta[old];
}
++new_index;
}
column_space_basis.swap(new_column_space_basis);
row_space_basis.swap(new_row_space_basis);
column_space_left_inverse.swap(new_column_space_left_inverse);
row_space_right_inverse.swap(new_row_space_right_inverse);
--matrix_rank;
return matrix_rank;
}
int apply_row_replacement(int row_index, const std::vector<T> &new_row) {
assert(0 <= row_index && row_index < row_size);
assert(static_cast<int>(new_row.size()) == column_size);
std::vector<T> difference = new_row;
const std::vector<T> current_row = get_row(row_index);
for (int j = 0; j < column_size; ++j) {
difference[j] -= current_row[j];
}
std::vector<T> column_vector(row_size, T());
column_vector[row_index] = T(1);
return apply_rank_one_update(column_vector, difference);
}
int apply_column_replacement(int column_index,
const std::vector<T> &new_column) {
assert(0 <= column_index && column_index < column_size);
assert(static_cast<int>(new_column.size()) == row_size);
std::vector<T> difference = new_column;
const std::vector<T> current_column = get_column(column_index);
for (int i = 0; i < row_size; ++i) {
difference[i] -= current_column[i];
}
std::vector<T> row_vector(column_size, T());
row_vector[column_index] = T(1);
return apply_rank_one_update(difference, row_vector);
}
private:
struct RankOneUpdateInfo {
std::vector<T> alpha;
std::vector<T> beta;
std::vector<T> column_residual;
std::vector<T> row_residual;
bool column_inside = false;
bool row_inside = false;
T schur = T();
int next_rank = 0;
};
RankOneUpdateInfo
analyze_rank_one_update(const std::vector<T> &column_vector,
const std::vector<T> &row_vector) const {
assert(static_cast<int>(column_vector.size()) == row_size);
assert(static_cast<int>(row_vector.size()) == column_size);
RankOneUpdateInfo info;
info.alpha = multiply_left_inverse(column_vector);
info.beta = multiply_right_inverse(row_vector);
const std::vector<T> projected_column =
reconstruct_column_vector(info.alpha);
const std::vector<T> projected_row = reconstruct_row_vector(info.beta);
info.column_residual.resize(row_size, T());
info.row_residual.resize(column_size, T());
info.column_inside = true;
info.row_inside = true;
for (int i = 0; i < row_size; ++i) {
info.column_residual[i] = column_vector[i] - projected_column[i];
if (info.column_residual[i] != T()) {
info.column_inside = false;
}
}
for (int j = 0; j < column_size; ++j) {
info.row_residual[j] = row_vector[j] - projected_row[j];
if (info.row_residual[j] != T()) {
info.row_inside = false;
}
}
info.schur = T(1);
for (int i = 0; i < matrix_rank; ++i) {
info.schur += info.beta[i] * info.alpha[i];
}
if (!info.column_inside && !info.row_inside) {
info.next_rank = matrix_rank + 1;
} else if (info.column_inside != info.row_inside) {
info.next_rank = matrix_rank;
} else {
info.next_rank = info.schur == T() ? matrix_rank - 1 : matrix_rank;
}
return info;
}
std::vector<T>
multiply_left_inverse(const std::vector<T> &column_vector) const {
std::vector<T> result(matrix_rank, T());
for (int i = 0; i < matrix_rank; ++i) {
for (int j = 0; j < row_size; ++j) {
result[i] += column_space_left_inverse[i][j] * column_vector[j];
}
}
return result;
}
std::vector<T>
multiply_right_inverse(const std::vector<T> &row_vector) const {
std::vector<T> result(matrix_rank, T());
for (int j = 0; j < column_size; ++j) {
const T value = row_vector[j];
if (value == T()) {
continue;
}
for (int i = 0; i < matrix_rank; ++i) {
result[i] += value * row_space_right_inverse[j][i];
}
}
return result;
}
std::vector<T>
reconstruct_column_vector(const std::vector<T> &coefficients) const {
std::vector<T> result(row_size, T());
for (int i = 0; i < row_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
result[i] += column_space_basis[i][j] * coefficients[j];
}
}
return result;
}
std::vector<T>
reconstruct_row_vector(const std::vector<T> &coefficients) const {
std::vector<T> result(column_size, T());
for (int i = 0; i < matrix_rank; ++i) {
const T value = coefficients[i];
if (value == T()) {
continue;
}
for (int j = 0; j < column_size; ++j) {
result[j] += value * row_space_basis[i][j];
}
}
return result;
}
std::vector<T> normalized_left_annihilator(int pivot_row,
const T &pivot_value) const {
std::vector<T> result(row_size, T());
result[pivot_row] = T(1);
for (int j = 0; j < matrix_rank; ++j) {
const T value = column_space_basis[pivot_row][j];
if (value == T()) {
continue;
}
for (int i = 0; i < row_size; ++i) {
result[i] -= value * column_space_left_inverse[j][i];
}
}
for (int i = 0; i < row_size; ++i) {
result[i] /= pivot_value;
}
return result;
}
std::vector<T> normalized_right_annihilator(int pivot_column,
const T &pivot_value) const {
std::vector<T> result(column_size, T());
result[pivot_column] = T(1);
for (int i = 0; i < column_size; ++i) {
for (int j = 0; j < matrix_rank; ++j) {
result[i] -= row_space_right_inverse[i][j] *
row_space_basis[j][pivot_column];
}
}
for (int i = 0; i < column_size; ++i) {
result[i] /= pivot_value;
}
return result;
}
static void append_column(std::vector<std::vector<T>> &matrix,
const std::vector<T> &column) {
assert(static_cast<int>(matrix.size()) ==
static_cast<int>(column.size()));
for (int i = 0; i < static_cast<int>(matrix.size()); ++i) {
matrix[i].push_back(column[i]);
}
}
static int first_nonzero(const std::vector<T> &vector) {
for (int i = 0; i < static_cast<int>(vector.size()); ++i) {
if (vector[i] != T()) {
return i;
}
}
assert(false);
return -1;
}
static std::vector<std::vector<T>>
transpose_matrix(const std::vector<std::vector<T>> &matrix) {
const int h = static_cast<int>(matrix.size());
if (h == 0) {
return {};
}
const int w = static_cast<int>(matrix[0].size());
std::vector<std::vector<T>> transposed(w, std::vector<T>(h, T()));
for (int i = 0; i < h; ++i) {
assert(static_cast<int>(matrix[i].size()) == w);
for (int j = 0; j < w; ++j) {
transposed[j][i] = matrix[i][j];
}
}
return transposed;
}
static std::vector<int>
find_independent_columns(const std::vector<std::vector<T>> &matrix) {
const int h = static_cast<int>(matrix.size());
if (h == 0) {
return {};
}
const int w = static_cast<int>(matrix[0].size());
std::vector<std::vector<T>> b = matrix;
int rank = 0;
std::vector<int> pivot_columns;
for (int column = 0; column < w; ++column) {
int pivot = -1;
for (int row = rank; row < h; ++row) {
if (b[row][column] != T()) {
pivot = row;
break;
}
}
if (pivot < 0) {
continue;
}
if (pivot != rank) {
b[pivot].swap(b[rank]);
}
const T inverse = T(1) / b[rank][column];
for (int j = column; j < w; ++j) {
b[rank][j] *= inverse;
}
for (int row = rank + 1; row < h; ++row) {
const T factor = b[row][column];
if (factor == T()) {
continue;
}
for (int j = column; j < w; ++j) {
b[row][j] -= b[rank][j] * factor;
}
}
pivot_columns.push_back(column);
++rank;
if (rank == h) {
break;
}
}
return pivot_columns;
}
static std::vector<std::vector<T>>
inverse_matrix(std::vector<std::vector<T>> matrix) {
const int n = static_cast<int>(matrix.size());
for (int i = 0; i < n; ++i) {
assert(static_cast<int>(matrix[i].size()) == n);
}
std::vector<std::vector<T>> inverse(n, std::vector<T>(n, T()));
for (int i = 0; i < n; ++i) {
inverse[i][i] = T(1);
}
for (int column = 0; column < n; ++column) {
int pivot = -1;
for (int row = column; row < n; ++row) {
if (matrix[row][column] != T()) {
pivot = row;
break;
}
}
assert(pivot >= 0);
if (pivot != column) {
matrix[pivot].swap(matrix[column]);
inverse[pivot].swap(inverse[column]);
}
const T diagonal_inverse = T(1) / matrix[column][column];
for (int j = 0; j < n; ++j) {
matrix[column][j] *= diagonal_inverse;
inverse[column][j] *= diagonal_inverse;
}
for (int row = 0; row < n; ++row) {
if (row == column) {
continue;
}
const T factor = matrix[row][column];
if (factor == T()) {
continue;
}
for (int j = 0; j < n; ++j) {
matrix[row][j] -= matrix[column][j] * factor;
inverse[row][j] -= inverse[column][j] * factor;
}
}
}
return inverse;
}
};
#line 7 "verify/standalone-dynamic-matrix-rank.test.cpp"
namespace {
struct F2 {
int value;
F2(int value = 0) : value(value & 1) {}
F2 &operator+=(const F2 &other) {
value ^= other.value;
return *this;
}
F2 &operator-=(const F2 &other) {
value ^= other.value;
return *this;
}
F2 &operator*=(const F2 &other) {
value &= other.value;
return *this;
}
F2 &operator/=(const F2 &other) {
assert(other.value == 1);
return *this;
}
[[maybe_unused]] friend F2 operator+(F2 lhs, const F2 &rhs) {
return lhs += rhs;
}
[[maybe_unused]] friend F2 operator-(F2 lhs, const F2 &rhs) {
return lhs -= rhs;
}
friend F2 operator*(F2 lhs, const F2 &rhs) { return lhs *= rhs; }
friend F2 operator/(F2 lhs, const F2 &rhs) { return lhs /= rhs; }
friend bool operator==(const F2 &lhs, const F2 &rhs) {
return lhs.value == rhs.value;
}
friend bool operator!=(const F2 &lhs, const F2 &rhs) {
return lhs.value != rhs.value;
}
};
template <class T> int brute_rank(std::vector<std::vector<T>> matrix) {
const int h = static_cast<int>(matrix.size());
if (h == 0) {
return 0;
}
const int w = static_cast<int>(matrix[0].size());
int rank = 0;
for (int column = 0; column < w; ++column) {
int pivot = -1;
for (int row = rank; row < h; ++row) {
if (matrix[row][column] != T()) {
pivot = row;
break;
}
}
if (pivot < 0) {
continue;
}
if (pivot != rank) {
matrix[pivot].swap(matrix[rank]);
}
const T inverse = T(1) / matrix[rank][column];
for (int j = column; j < w; ++j) {
matrix[rank][j] *= inverse;
}
for (int row = rank + 1; row < h; ++row) {
const T factor = matrix[row][column];
if (factor == T()) {
continue;
}
for (int j = column; j < w; ++j) {
matrix[row][j] -= matrix[rank][j] * factor;
}
}
++rank;
if (rank == h) {
break;
}
}
return rank;
}
std::vector<std::vector<F2>> matrix_from_mask(int h, int w, int mask) {
std::vector<std::vector<F2>> matrix(h, std::vector<F2>(w, F2()));
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
matrix[i][j] = F2((mask >> (i * w + j)) & 1);
}
}
return matrix;
}
std::vector<F2> vector_from_mask(int n, int mask) {
std::vector<F2> values(n, F2());
for (int i = 0; i < n; ++i) {
values[i] = F2((mask >> i) & 1);
}
return values;
}
std::vector<std::vector<F2>>
apply_rank_one(std::vector<std::vector<F2>> matrix,
const std::vector<F2> &column_vector,
const std::vector<F2> &row_vector) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
matrix[i][j] += column_vector[i] * row_vector[j];
}
}
return matrix;
}
void check_solver(const DynamicMatrixRank<F2> &solver,
const std::vector<std::vector<F2>> &matrix) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
assert(solver.rank() == brute_rank(matrix));
assert(solver.materialize_matrix() == matrix);
for (int row = 0; row < h; ++row) {
assert(solver.get_row(row) == matrix[row]);
}
for (int column = 0; column < w; ++column) {
std::vector<F2> expected(h, F2());
for (int row = 0; row < h; ++row) {
expected[row] = matrix[row][column];
}
assert(solver.get_column(column) == expected);
}
}
void check_updates(const std::vector<std::vector<F2>> &matrix) {
const int h = static_cast<int>(matrix.size());
const int w = h == 0 ? 0 : static_cast<int>(matrix[0].size());
DynamicMatrixRank<F2> solver(matrix);
check_solver(solver, matrix);
for (int column_mask = 0; column_mask < (1 << h); ++column_mask) {
const std::vector<F2> column_vector = vector_from_mask(h, column_mask);
for (int row_mask = 0; row_mask < (1 << w); ++row_mask) {
const std::vector<F2> row_vector = vector_from_mask(w, row_mask);
const auto updated =
apply_rank_one(matrix, column_vector, row_vector);
const int expected = brute_rank(updated);
assert(solver.rank_after_rank_one_update(column_vector,
row_vector) == expected);
auto applied = solver;
assert(applied.apply_rank_one_update(column_vector, row_vector) ==
expected);
check_solver(applied, updated);
applied.build();
check_solver(applied, updated);
for (int row = 0; row < h; ++row) {
for (int new_row_mask = 0; new_row_mask < (1 << w);
++new_row_mask) {
const std::vector<F2> new_row =
vector_from_mask(w, new_row_mask);
auto row_updated = updated;
row_updated[row] = new_row;
const int row_expected = brute_rank(row_updated);
assert(applied.rank_after_row_replacement(row, new_row) ==
row_expected);
auto row_applied = applied;
assert(row_applied.apply_row_replacement(row, new_row) ==
row_expected);
check_solver(row_applied, row_updated);
}
}
for (int column = 0; column < w; ++column) {
for (int new_column_mask = 0; new_column_mask < (1 << h);
++new_column_mask) {
const std::vector<F2> new_column =
vector_from_mask(h, new_column_mask);
auto column_updated = updated;
for (int row = 0; row < h; ++row) {
column_updated[row][column] = new_column[row];
}
const int column_expected = brute_rank(column_updated);
assert(applied.rank_after_column_replacement(
column, new_column) == column_expected);
auto column_applied = applied;
assert(column_applied.apply_column_replacement(
column, new_column) == column_expected);
check_solver(column_applied, column_updated);
}
}
}
}
}
void self_test() {
{
DynamicMatrixRank<F2> solver;
solver.build(std::vector<std::vector<F2>>{});
assert(solver.rank() == 0);
assert(solver.materialize_matrix().empty());
}
for (int h = 1; h <= 3; ++h) {
for (int w = 1; w <= 3; ++w) {
const int states = 1 << (h * w);
for (int mask = 0; mask < states; ++mask) {
check_updates(matrix_from_mask(h, w, mask));
}
}
}
}
} // namespace
int main() {
self_test();
return 0;
}