linear_algebra.hpp¶
在域上求矩阵秩、行列式、线性方程组解与零空间基；适合一般模线性代数题。
\[
\displaystyle A\mathbf{x}=\mathbf{b},\qquad \{\mathbf{x}\}=\mathbf{x}_0+\ker A
\]
Complexity: Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse. Space: O(rows * columns).
AC 记录：matrix_det, matrix_rank, system_of_linear_equations。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse.
/// Space: O(rows * columns).
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Exact zero predicate used by elimination routines by default.
template <class T> struct exact_zero {
bool operator()(const T &val) const { return val == T{}; }
};
/// @brief Consistency flag, one solution, nullspace basis, and pivot columns.
template <class T> struct linear_system_solution {
bool ok = false;
std::vector<T> sol;
std::vector<std::vector<T>> ker;
std::vector<int> pc;
};
namespace linear_algebra_internal {
template <class T> int column_count(const std::vector<std::vector<T>> &mat) {
if (mat.empty()) {
return 0;
}
int cs = int(mat[0].size());
for (const auto &row : mat) {
assert(int(row.size()) == cs);
}
return cs;
}
} // namespace linear_algebra_internal
/// @brief Compute matrix rank over a field.
template <class T, class IsZero = exact_zero<T>>
int matrix_rank(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int rs = int(mat.size());
int cs = linear_algebra_internal::column_count(mat);
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
}
ran++;
}
return ran;
}
/// @brief Compute the determinant of a square matrix over a field.
template <class T, class IsZero = exact_zero<T>>
T determinant(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int n = int(mat.size());
assert(linear_algebra_internal::column_count(mat) == n);
T res = T(1);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(mat[piv][col])) {
piv++;
}
if (piv == n) {
return T{};
}
if (piv != col) {
std::swap(mat[piv], mat[col]);
res = -res;
}
T pv = mat[col][col];
res *= pv;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / pv;
for (int j = col; j < n; j++) {
mat[row][j] -= r * mat[col][j];
}
}
}
return res;
}
/// @brief Solve A*x=b and return one solution plus a basis of the nullspace.
template <class T, class IsZero = exact_zero<T>>
linear_system_solution<T> solve_linear(std::vector<std::vector<T>> mat,
std::vector<T> rhs, IsZero iz = {}) {
int rs = int(mat.size());
assert(int(rhs.size()) == rs);
int cs = linear_algebra_internal::column_count(mat);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
std::swap(rhs[piv], rhs[ran]);
T inv = T(1) / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[ran][j] *= inv;
}
rhs[ran] *= inv;
for (int row = 0; row < rs; row++) {
if (row == ran || iz(mat[row][col])) {
continue;
}
T r = mat[row][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
rhs[row] -= r * rhs[ran];
}
pc.push_back(col);
ran++;
}
for (int row = ran; row < rs; row++) {
bool az = true;
for (int col = 0; col < cs; col++) {
az &= iz(mat[row][col]);
}
if (az && !iz(rhs[row])) {
return {};
}
}
linear_system_solution<T> res;
res.ok = true;
res.sol.assign(cs, T{});
res.pc = pc;
std::vector<bool> ip(cs);
for (int row = 0; row < ran; row++) {
int col = pc[row];
ip[col] = true;
res.sol[col] = rhs[row];
}
for (int fc = 0; fc < cs; fc++) {
if (ip[fc]) {
continue;
}
std::vector<T> bv(cs, T{});
bv[fc] = T(1);
for (int row = 0; row < ran; row++) {
bv[pc[row]] = -mat[row][fc];
}
res.ker.push_back(std::move(bv));
}
return res;
}
} // namespace noya
#ifndef NOYA_LINEAR_ALGEBRA_HPP
#define NOYA_LINEAR_ALGEBRA_HPP 1
/// @complexity Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse.
/// Space: O(rows * columns).
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Exact zero predicate used by elimination routines by default.
template <class T> struct exact_zero {
bool operator()(const T &val) const { return val == T{}; }
};
/// @brief Consistency flag, one solution, nullspace basis, and pivot columns.
template <class T> struct linear_system_solution {
bool ok = false;
std::vector<T> sol;
std::vector<std::vector<T>> ker;
std::vector<int> pc;
};
namespace linear_algebra_internal {
template <class T> int column_count(const std::vector<std::vector<T>> &mat) {
if (mat.empty()) {
return 0;
}
int cs = int(mat[0].size());
for (const auto &row : mat) {
assert(int(row.size()) == cs);
}
return cs;
}
} // namespace linear_algebra_internal
/// @brief Compute matrix rank over a field.
template <class T, class IsZero = exact_zero<T>>
int matrix_rank(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int rs = int(mat.size());
int cs = linear_algebra_internal::column_count(mat);
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
}
ran++;
}
return ran;
}
/// @brief Compute the determinant of a square matrix over a field.
template <class T, class IsZero = exact_zero<T>>
T determinant(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int n = int(mat.size());
assert(linear_algebra_internal::column_count(mat) == n);
T res = T(1);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(mat[piv][col])) {
piv++;
}
if (piv == n) {
return T{};
}
if (piv != col) {
std::swap(mat[piv], mat[col]);
res = -res;
}
T pv = mat[col][col];
res *= pv;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / pv;
for (int j = col; j < n; j++) {
mat[row][j] -= r * mat[col][j];
}
}
}
return res;
}
/// @brief Solve A*x=b and return one solution plus a basis of the nullspace.
template <class T, class IsZero = exact_zero<T>>
linear_system_solution<T> solve_linear(std::vector<std::vector<T>> mat,
std::vector<T> rhs, IsZero iz = {}) {
int rs = int(mat.size());
assert(int(rhs.size()) == rs);
int cs = linear_algebra_internal::column_count(mat);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
std::swap(rhs[piv], rhs[ran]);
T inv = T(1) / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[ran][j] *= inv;
}
rhs[ran] *= inv;
for (int row = 0; row < rs; row++) {
if (row == ran || iz(mat[row][col])) {
continue;
}
T r = mat[row][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
rhs[row] -= r * rhs[ran];
}
pc.push_back(col);
ran++;
}
for (int row = ran; row < rs; row++) {
bool az = true;
for (int col = 0; col < cs; col++) {
az &= iz(mat[row][col]);
}
if (az && !iz(rhs[row])) {
return {};
}
}
linear_system_solution<T> res;
res.ok = true;
res.sol.assign(cs, T{});
res.pc = pc;
std::vector<bool> ip(cs);
for (int row = 0; row < ran; row++) {
int col = pc[row];
ip[col] = true;
res.sol[col] = rhs[row];
}
for (int fc = 0; fc < cs; fc++) {
if (ip[fc]) {
continue;
}
std::vector<T> bv(cs, T{});
bv[fc] = T(1);
for (int row = 0; row < ran; row++) {
bv[pc[row]] = -mat[row][fc];
}
res.ker.push_back(std::move(bv));
}
return res;
}
} // namespace noya
#endif // NOYA_LINEAR_ALGEBRA_HPP
#include <algorithm>
#include <cassert>
#include <utility>
#include <vector>
/// @complexity Time: O(rows * columns * min(rows,columns)) elimination; O(n^3) square determinant/inverse.
/// Space: O(rows * columns).
namespace noya {
/// @brief Exact zero predicate used by elimination routines by default.
template <class T> struct exact_zero {
bool operator()(const T &val) const { return val == T{}; }
};
/// @brief Consistency flag, one solution, nullspace basis, and pivot columns.
template <class T> struct linear_system_solution {
bool ok = false;
std::vector<T> sol;
std::vector<std::vector<T>> ker;
std::vector<int> pc;
};
namespace linear_algebra_internal {
template <class T> int column_count(const std::vector<std::vector<T>> &mat) {
if (mat.empty()) {
return 0;
}
int cs = int(mat[0].size());
for (const auto &row : mat) {
assert(int(row.size()) == cs);
}
return cs;
}
} // namespace linear_algebra_internal
/// @brief Compute matrix rank over a field.
template <class T, class IsZero = exact_zero<T>>
int matrix_rank(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int rs = int(mat.size());
int cs = linear_algebra_internal::column_count(mat);
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
for (int row = ran + 1; row < rs; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
}
ran++;
}
return ran;
}
/// @brief Compute the determinant of a square matrix over a field.
template <class T, class IsZero = exact_zero<T>>
T determinant(std::vector<std::vector<T>> mat, IsZero iz = {}) {
int n = int(mat.size());
assert(linear_algebra_internal::column_count(mat) == n);
T res = T(1);
for (int col = 0; col < n; col++) {
int piv = col;
while (piv < n && iz(mat[piv][col])) {
piv++;
}
if (piv == n) {
return T{};
}
if (piv != col) {
std::swap(mat[piv], mat[col]);
res = -res;
}
T pv = mat[col][col];
res *= pv;
for (int row = col + 1; row < n; row++) {
if (iz(mat[row][col])) {
continue;
}
T r = mat[row][col] / pv;
for (int j = col; j < n; j++) {
mat[row][j] -= r * mat[col][j];
}
}
}
return res;
}
/// @brief Solve A*x=b and return one solution plus a basis of the nullspace.
template <class T, class IsZero = exact_zero<T>>
linear_system_solution<T> solve_linear(std::vector<std::vector<T>> mat,
std::vector<T> rhs, IsZero iz = {}) {
int rs = int(mat.size());
assert(int(rhs.size()) == rs);
int cs = linear_algebra_internal::column_count(mat);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < cs && ran < rs; col++) {
int piv = ran;
while (piv < rs && iz(mat[piv][col])) {
piv++;
}
if (piv == rs) {
continue;
}
std::swap(mat[piv], mat[ran]);
std::swap(rhs[piv], rhs[ran]);
T inv = T(1) / mat[ran][col];
for (int j = col; j < cs; j++) {
mat[ran][j] *= inv;
}
rhs[ran] *= inv;
for (int row = 0; row < rs; row++) {
if (row == ran || iz(mat[row][col])) {
continue;
}
T r = mat[row][col];
for (int j = col; j < cs; j++) {
mat[row][j] -= r * mat[ran][j];
}
rhs[row] -= r * rhs[ran];
}
pc.push_back(col);
ran++;
}
for (int row = ran; row < rs; row++) {
bool az = true;
for (int col = 0; col < cs; col++) {
az &= iz(mat[row][col]);
}
if (az && !iz(rhs[row])) {
return {};
}
}
linear_system_solution<T> res;
res.ok = true;
res.sol.assign(cs, T{});
res.pc = pc;
std::vector<bool> ip(cs);
for (int row = 0; row < ran; row++) {
int col = pc[row];
ip[col] = true;
res.sol[col] = rhs[row];
}
for (int fc = 0; fc < cs; fc++) {
if (ip[fc]) {
continue;
}
std::vector<T> bv(cs, T{});
bv[fc] = T(1);
for (int row = 0; row < ran; row++) {
bv[pc[row]] = -mat[row][fc];
}
res.ker.push_back(std::move(bv));
}
return res;
}
} // namespace noya