f2_matrix.hpp¶
在 \(\mathrm{GF}(2)\) 上用位集完成矩阵消元、秩和线性方程；适合异或方程组。
\[
\displaystyle \mathrm{rank}_{\mathbf F_2}(A)
\]
Complexity: Time: O(rows * columns * min(rows,columns) / 64) for elimination, O(rows * inner * columns / 64) for multiplication. Space: O(rows * columns / 64).
AC 记录：inverse_matrix_mod_2, matrix_det_mod_2, matrix_product_mod_2, matrix_rank_mod_2, system_of_linear_equations_mod_2。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @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 &mat) {
if (mat.empty()) {
return 0;
}
std::size_t res = mat.front().size();
for (const auto &row : mat) {
assert(row.size() == res);
}
return res;
}
inline f2_matrix transpose(const f2_matrix &mat, std::size_t cs) {
f2_matrix res(cs, dynamic_bitset(mat.size()));
for (std::size_t row = 0; row < mat.size(); row++) {
assert(mat[row].size() == cs);
for (std::size_t col = mat[row].find_first(); col < cs;
col = mat[row].find_next(col + 1)) {
res[col].set(row);
}
}
return res;
}
} // 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 mat, std::size_t cs) {
for (const auto &row : mat) {
assert(row.size() == cs);
}
if (mat.size() > cs) {
mat = f2_matrix_detail::transpose(mat, cs);
cs = mat.empty() ? 0 : mat.front().size();
}
int ran = 0;
for (std::size_t col = 0; col < cs && ran < int(mat.size()); col++) {
int piv = ran;
while (piv < int(mat.size()) && !mat[piv][col]) {
piv++;
}
if (piv == int(mat.size())) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < int(mat.size()); row++) {
if (mat[row][col]) {
mat[row] ^= mat[ran];
}
}
ran++;
}
return ran;
}
inline int f2_matrix_rank(f2_matrix mat) {
std::size_t cs = f2_matrix_detail::columns(mat);
return f2_matrix_rank(std::move(mat), cs);
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return false;
}
std::swap(mat[piv], mat[col]);
for (int row = col + 1; row < sz; row++) {
if (mat[row][col]) {
mat[row] ^= mat[col];
}
}
}
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 &l, const f2_matrix &r,
std::size_t ni, std::size_t no) {
for (const auto &row : l) {
assert(row.size() == ni);
}
assert(r.size() == ni);
for (const auto &row : r) {
assert(row.size() == no);
}
f2_matrix res(l.size(), dynamic_bitset(no));
for (std::size_t row = 0; row < l.size(); row++) {
for (std::size_t idx = l[row].find_first(); idx < ni;
idx = l[row].find_next(idx + 1)) {
res[row] ^= r[idx];
}
}
return res;
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
f2_matrix inv(sz, dynamic_bitset(sz));
for (int row = 0; row < sz; row++) {
inv[row].set(row);
}
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return std::nullopt;
}
std::swap(mat[piv], mat[col]);
std::swap(inv[piv], inv[col]);
for (int row = 0; row < sz; row++) {
if (row != col && mat[row][col]) {
mat[row] ^= mat[col];
inv[row] ^= inv[col];
}
}
}
return inv;
}
struct f2_linear_system_solution {
bool ok = false;
dynamic_bitset sol;
f2_matrix ker;
std::vector<int> pc;
};
/// @brief Solve `mat * x = rhs` 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 mat, std::size_t cs,
const dynamic_bitset &rhs) {
int rs = int(mat.size());
assert(rhs.size() == std::size_t(rs));
for (const auto &row : mat) {
assert(row.size() == cs);
}
std::vector<bool> r(rs);
for (int row = 0; row < rs; row++) {
r[row] = rhs[row];
}
std::vector<int> pvt;
int ran = 0;
for (std::size_t col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && !mat[piv][col]) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
bool tmp = r[piv];
r[piv] = r[ran];
r[ran] = tmp;
for (int row = 0; row < rs; row++) {
if (row != ran && mat[row][col]) {
mat[row] ^= mat[ran];
r[row] = r[row] != r[ran];
}
}
pvt.push_back(int(col));
ran++;
}
for (int row = ran; row < rs; row++) {
if (r[row]) {
return {false, dynamic_bitset(cs), {}, std::move(pvt)};
}
}
f2_linear_system_solution res;
res.ok = true;
res.sol = dynamic_bitset(cs);
res.pc = pvt;
std::vector<int> pr(cs, -1);
for (int row = 0; row < ran; row++) {
int col = pvt[row];
pr[col] = row;
res.sol.set(col, r[row]);
}
for (std::size_t fc = 0; fc < cs; fc++) {
if (pr[fc] != -1) {
continue;
}
dynamic_bitset bas(cs);
bas.set(fc);
for (int row = 0; row < ran; row++) {
if (mat[row][fc]) {
bas.set(pvt[row]);
}
}
res.ker.push_back(std::move(bas));
}
return res;
}
} // namespace noya
#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 &mat) {
if (mat.empty()) {
return 0;
}
std::size_t res = mat.front().size();
for (const auto &row : mat) {
assert(row.size() == res);
}
return res;
}
inline f2_matrix transpose(const f2_matrix &mat, std::size_t cs) {
f2_matrix res(cs, dynamic_bitset(mat.size()));
for (std::size_t row = 0; row < mat.size(); row++) {
assert(mat[row].size() == cs);
for (std::size_t col = mat[row].find_first(); col < cs;
col = mat[row].find_next(col + 1)) {
res[col].set(row);
}
}
return res;
}
} // 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 mat, std::size_t cs) {
for (const auto &row : mat) {
assert(row.size() == cs);
}
if (mat.size() > cs) {
mat = f2_matrix_detail::transpose(mat, cs);
cs = mat.empty() ? 0 : mat.front().size();
}
int ran = 0;
for (std::size_t col = 0; col < cs && ran < int(mat.size()); col++) {
int piv = ran;
while (piv < int(mat.size()) && !mat[piv][col]) {
piv++;
}
if (piv == int(mat.size())) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < int(mat.size()); row++) {
if (mat[row][col]) {
mat[row] ^= mat[ran];
}
}
ran++;
}
return ran;
}
inline int f2_matrix_rank(f2_matrix mat) {
std::size_t cs = f2_matrix_detail::columns(mat);
return f2_matrix_rank(std::move(mat), cs);
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return false;
}
std::swap(mat[piv], mat[col]);
for (int row = col + 1; row < sz; row++) {
if (mat[row][col]) {
mat[row] ^= mat[col];
}
}
}
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 &l, const f2_matrix &r,
std::size_t ni, std::size_t no) {
for (const auto &row : l) {
assert(row.size() == ni);
}
assert(r.size() == ni);
for (const auto &row : r) {
assert(row.size() == no);
}
f2_matrix res(l.size(), dynamic_bitset(no));
for (std::size_t row = 0; row < l.size(); row++) {
for (std::size_t idx = l[row].find_first(); idx < ni;
idx = l[row].find_next(idx + 1)) {
res[row] ^= r[idx];
}
}
return res;
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
f2_matrix inv(sz, dynamic_bitset(sz));
for (int row = 0; row < sz; row++) {
inv[row].set(row);
}
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return std::nullopt;
}
std::swap(mat[piv], mat[col]);
std::swap(inv[piv], inv[col]);
for (int row = 0; row < sz; row++) {
if (row != col && mat[row][col]) {
mat[row] ^= mat[col];
inv[row] ^= inv[col];
}
}
}
return inv;
}
struct f2_linear_system_solution {
bool ok = false;
dynamic_bitset sol;
f2_matrix ker;
std::vector<int> pc;
};
/// @brief Solve `mat * x = rhs` 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 mat, std::size_t cs,
const dynamic_bitset &rhs) {
int rs = int(mat.size());
assert(rhs.size() == std::size_t(rs));
for (const auto &row : mat) {
assert(row.size() == cs);
}
std::vector<bool> r(rs);
for (int row = 0; row < rs; row++) {
r[row] = rhs[row];
}
std::vector<int> pvt;
int ran = 0;
for (std::size_t col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && !mat[piv][col]) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
bool tmp = r[piv];
r[piv] = r[ran];
r[ran] = tmp;
for (int row = 0; row < rs; row++) {
if (row != ran && mat[row][col]) {
mat[row] ^= mat[ran];
r[row] = r[row] != r[ran];
}
}
pvt.push_back(int(col));
ran++;
}
for (int row = ran; row < rs; row++) {
if (r[row]) {
return {false, dynamic_bitset(cs), {}, std::move(pvt)};
}
}
f2_linear_system_solution res;
res.ok = true;
res.sol = dynamic_bitset(cs);
res.pc = pvt;
std::vector<int> pr(cs, -1);
for (int row = 0; row < ran; row++) {
int col = pvt[row];
pr[col] = row;
res.sol.set(col, r[row]);
}
for (std::size_t fc = 0; fc < cs; fc++) {
if (pr[fc] != -1) {
continue;
}
dynamic_bitset bas(cs);
bas.set(fc);
for (int row = 0; row < ran; row++) {
if (mat[row][fc]) {
bas.set(pvt[row]);
}
}
res.ker.push_back(std::move(bas));
}
return res;
}
} // 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 W = 64;
std::size_t sz_ = 0;
std::vector<word_type> a;
dynamic_bitset() = default;
explicit dynamic_bitset(std::size_t siz, bool val = false)
: sz_(siz), a(word_count(siz), val ? ~word_type{} : 0) {
trim();
}
std::size_t size() const { return sz_; }
bool empty() const { return sz_ == 0; }
bool test(std::size_t pos) const {
assert(pos < sz_);
return (a[pos / W] >> (pos % W)) & 1;
}
bool operator[](std::size_t pos) const { return test(pos); }
dynamic_bitset &set(std::size_t pos, bool val = true) {
assert(pos < sz_);
word_type msk = word_type(1) << (pos % W);
if (val) {
a[pos / W] |= msk;
} else {
a[pos / W] &= ~msk;
}
return *this;
}
dynamic_bitset &reset(std::size_t pos) { return set(pos, false); }
dynamic_bitset &flip(std::size_t pos) {
assert(pos < sz_);
a[pos / W] ^= word_type(1) << (pos % W);
return *this;
}
dynamic_bitset &set() {
std::fill(a.begin(), a.end(), ~word_type{});
trim();
return *this;
}
dynamic_bitset &reset() {
std::fill(a.begin(), a.end(), word_type{});
return *this;
}
dynamic_bitset &flip() {
for (word_type &w : a) {
w = ~w;
}
trim();
return *this;
}
std::size_t count() const {
std::size_t res = 0;
for (word_type w : a) {
res += std::popcount(w);
}
return res;
}
bool any() const {
return std::any_of(a.begin(), a.end(), [](word_type w) { return w != 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 pos) const {
if (pos >= sz_) {
return sz_;
}
std::size_t idx = pos / W;
word_type w = a[idx] & (~word_type{} << (pos % W));
if (w != 0) {
return std::min(sz_, idx * W + std::size_t(std::countr_zero(w)));
}
for (idx++; idx < a.size(); idx++) {
if (a[idx] != 0) {
return std::min(sz_, idx * W + std::size_t(std::countr_zero(a[idx])));
}
}
return sz_;
}
std::size_t find_first() const { return find_next(0); }
dynamic_bitset &operator&=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] &= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator|=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] |= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator^=(const dynamic_bitset &rhs) {
check_same_size(rhs);
for (std::size_t i = 0; i < a.size(); i++) {
a[i] ^= rhs.a[i];
}
return *this;
}
dynamic_bitset &operator<<=(std::size_t k) {
if (k >= sz_) {
return reset();
}
std::size_t blk = k / W;
int rem = int(k % W);
for (std::size_t i = a.size(); i-- > 0;) {
word_type val = 0;
if (i >= blk) {
val = a[i - blk] << rem;
if (rem != 0 && i > blk) {
val |= a[i - blk - 1] >> (W - rem);
}
}
a[i] = val;
}
trim();
return *this;
}
dynamic_bitset &operator>>=(std::size_t k) {
if (k >= sz_) {
return reset();
}
std::size_t blk = k / W;
int rem = int(k % W);
for (std::size_t i = 0; i < a.size(); i++) {
word_type val = 0;
if (i + blk < a.size()) {
val = a[i + blk] >> rem;
if (rem != 0 && i + blk + 1 < a.size()) {
val |= a[i + blk + 1] << (W - rem);
}
}
a[i] = val;
}
return *this;
}
friend dynamic_bitset operator&(dynamic_bitset arr, const dynamic_bitset &b) {
return arr &= b;
}
friend dynamic_bitset operator|(dynamic_bitset arr, const dynamic_bitset &b) {
return arr |= b;
}
friend dynamic_bitset operator^(dynamic_bitset arr, const dynamic_bitset &b) {
return arr ^= b;
}
friend dynamic_bitset operator<<(dynamic_bitset val, std::size_t k) {
return val <<= k;
}
friend dynamic_bitset operator>>(dynamic_bitset val, std::size_t k) {
return val >>= k;
}
friend dynamic_bitset operator~(dynamic_bitset val) { return val.flip(); }
friend bool operator==(const dynamic_bitset &,
const dynamic_bitset &) = default;
private:
static std::size_t word_count(std::size_t siz) { return (siz + W - 1) / W; }
void trim() {
if (!a.empty() && sz_ % W != 0) {
a.back() &= (word_type(1) << (sz_ % W)) - word_type(1);
}
}
void check_same_size(const dynamic_bitset &rhs) const {
assert(sz_ == rhs.sz_);
}
};
} // namespace noya
namespace noya {
using f2_matrix = std::vector<dynamic_bitset>;
namespace f2_matrix_detail {
inline std::size_t columns(const f2_matrix &mat) {
if (mat.empty()) {
return 0;
}
std::size_t res = mat.front().size();
for (const auto &row : mat) {
assert(row.size() == res);
}
return res;
}
inline f2_matrix transpose(const f2_matrix &mat, std::size_t cs) {
f2_matrix res(cs, dynamic_bitset(mat.size()));
for (std::size_t row = 0; row < mat.size(); row++) {
assert(mat[row].size() == cs);
for (std::size_t col = mat[row].find_first(); col < cs;
col = mat[row].find_next(col + 1)) {
res[col].set(row);
}
}
return res;
}
} // 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 mat, std::size_t cs) {
for (const auto &row : mat) {
assert(row.size() == cs);
}
if (mat.size() > cs) {
mat = f2_matrix_detail::transpose(mat, cs);
cs = mat.empty() ? 0 : mat.front().size();
}
int ran = 0;
for (std::size_t col = 0; col < cs && ran < int(mat.size()); col++) {
int piv = ran;
while (piv < int(mat.size()) && !mat[piv][col]) {
piv++;
}
if (piv == int(mat.size())) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < int(mat.size()); row++) {
if (mat[row][col]) {
mat[row] ^= mat[ran];
}
}
ran++;
}
return ran;
}
inline int f2_matrix_rank(f2_matrix mat) {
std::size_t cs = f2_matrix_detail::columns(mat);
return f2_matrix_rank(std::move(mat), cs);
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return false;
}
std::swap(mat[piv], mat[col]);
for (int row = col + 1; row < sz; row++) {
if (mat[row][col]) {
mat[row] ^= mat[col];
}
}
}
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 &l, const f2_matrix &r,
std::size_t ni, std::size_t no) {
for (const auto &row : l) {
assert(row.size() == ni);
}
assert(r.size() == ni);
for (const auto &row : r) {
assert(row.size() == no);
}
f2_matrix res(l.size(), dynamic_bitset(no));
for (std::size_t row = 0; row < l.size(); row++) {
for (std::size_t idx = l[row].find_first(); idx < ni;
idx = l[row].find_next(idx + 1)) {
res[row] ^= r[idx];
}
}
return res;
}
/// @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 mat) {
int sz = int(mat.size());
assert(f2_matrix_detail::columns(mat) == std::size_t(sz));
f2_matrix inv(sz, dynamic_bitset(sz));
for (int row = 0; row < sz; row++) {
inv[row].set(row);
}
for (int col = 0; col < sz; col++) {
int piv = col;
while (piv < sz && !mat[piv][col]) {
piv++;
}
if (piv == sz) {
return std::nullopt;
}
std::swap(mat[piv], mat[col]);
std::swap(inv[piv], inv[col]);
for (int row = 0; row < sz; row++) {
if (row != col && mat[row][col]) {
mat[row] ^= mat[col];
inv[row] ^= inv[col];
}
}
}
return inv;
}
struct f2_linear_system_solution {
bool ok = false;
dynamic_bitset sol;
f2_matrix ker;
std::vector<int> pc;
};
/// @brief Solve `mat * x = rhs` 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 mat, std::size_t cs,
const dynamic_bitset &rhs) {
int rs = int(mat.size());
assert(rhs.size() == std::size_t(rs));
for (const auto &row : mat) {
assert(row.size() == cs);
}
std::vector<bool> r(rs);
for (int row = 0; row < rs; row++) {
r[row] = rhs[row];
}
std::vector<int> pvt;
int ran = 0;
for (std::size_t col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && !mat[piv][col]) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
bool tmp = r[piv];
r[piv] = r[ran];
r[ran] = tmp;
for (int row = 0; row < rs; row++) {
if (row != ran && mat[row][col]) {
mat[row] ^= mat[ran];
r[row] = r[row] != r[ran];
}
}
pvt.push_back(int(col));
ran++;
}
for (int row = ran; row < rs; row++) {
if (r[row]) {
return {false, dynamic_bitset(cs), {}, std::move(pvt)};
}
}
f2_linear_system_solution res;
res.ok = true;
res.sol = dynamic_bitset(cs);
res.pc = pvt;
std::vector<int> pr(cs, -1);
for (int row = 0; row < ran; row++) {
int col = pvt[row];
pr[col] = row;
res.sol.set(col, r[row]);
}
for (std::size_t fc = 0; fc < cs; fc++) {
if (pr[fc] != -1) {
continue;
}
dynamic_bitset bas(cs);
bas.set(fc);
for (int row = 0; row < ran; row++) {
if (mat[row][fc]) {
bas.set(pvt[row]);
}
}
res.ker.push_back(std::move(bas));
}
return res;
}
} // namespace noya