adjugate_matrix.hpp¶
计算矩阵的伴随矩阵，即使矩阵奇异也能得到所有代数余子式。
\[
\displaystyle \operatorname{adj}(A)A=\det(A)I
\]
Complexity: Time: \(O(n^{3})\). Space: \(O(n^{2})\).
AC 记录：adjugate_matrix。
Implementation¶
当前头文件，省略 include guard；依赖见 #include。
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <vector>
namespace noya {
/// @brief Compute the adjugate of a square matrix over a field.
/// Forward elimination tracks E A = U. Full rank uses back substitution for
/// A^{-1}; rank n-1 uses a right-null vector and E's last (left-null) row.
/// Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
for (int row = 0; row < sz; row++) {
ops[row][row] = T(1);
}
T det = T(1);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < sz && ran < sz; col++) {
int piv = ran;
while (piv < sz && mat[piv][col] == T{}) {
piv++;
}
if (piv == sz) {
continue;
}
if (piv != ran) {
std::swap(mat[piv], mat[ran]);
std::swap(ops[piv], ops[ran]);
det = -det;
}
T pv = mat[ran][col];
det *= pv;
T inv = T(1) / pv;
for (int nxt = col; nxt < sz; nxt++) {
mat[ran][nxt] *= inv;
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[ran][nxt] *= inv;
}
for (int row = ran + 1; row < sz; row++) {
if (mat[row][col] == T{}) {
continue;
}
T r = mat[row][col];
for (int nxt = col; nxt < sz; nxt++) {
mat[row][nxt] -= r * mat[ran][nxt];
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[row][nxt] -= r * ops[ran][nxt];
}
}
pc.push_back(col);
ran++;
}
if (ran == sz) {
// Solve U * A^{-1} = E; U itself need not be reduced above its pivots.
for (int col = sz - 1; col >= 0; --col) {
for (int row = 0; row < col; ++row) {
T r = mat[row][col];
if (r == T{}) continue;
for (int nxt = 0; nxt < sz; ++nxt) ops[row][nxt] -= r * ops[col][nxt];
}
}
for (auto &row : ops) for (auto &value : row) value *= det;
return ops;
}
std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
if (ran + 1 != sz) {
return res;
}
std::vector<bool> ip(sz);
for (int col : pc) {
ip[col] = true;
}
int fc = 0;
while (ip[fc]) {
fc++;
}
std::vector<T> rv(sz, T{});
rv[fc] = T(1);
for (int row = ran - 1; row >= 0; --row) {
for (int col = pc[row] + 1; col < sz; ++col) {
rv[pc[row]] -= mat[row][col] * rv[col];
}
}
T scl = det;
if ((ran + fc) % 2 != 0) {
scl = -scl;
}
for (int row = 0; row < sz; row++) {
for (int col = 0; col < sz; col++) {
res[row][col] = scl * rv[row] * ops[ran][col];
}
}
return res;
}
} // namespace noya
// NOYA_ADJUGATE_MATRIX_HPP
#ifndef NOYA_ADJUGATE_MATRIX_HPP
#define NOYA_ADJUGATE_MATRIX_HPP 1
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <vector>
namespace noya {
/// @brief Compute the adjugate of a square matrix over a field.
/// Forward elimination tracks E A = U. Full rank uses back substitution for
/// A^{-1}; rank n-1 uses a right-null vector and E's last (left-null) row.
/// Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
for (int row = 0; row < sz; row++) {
ops[row][row] = T(1);
}
T det = T(1);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < sz && ran < sz; col++) {
int piv = ran;
while (piv < sz && mat[piv][col] == T{}) {
piv++;
}
if (piv == sz) {
continue;
}
if (piv != ran) {
std::swap(mat[piv], mat[ran]);
std::swap(ops[piv], ops[ran]);
det = -det;
}
T pv = mat[ran][col];
det *= pv;
T inv = T(1) / pv;
for (int nxt = col; nxt < sz; nxt++) {
mat[ran][nxt] *= inv;
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[ran][nxt] *= inv;
}
for (int row = ran + 1; row < sz; row++) {
if (mat[row][col] == T{}) {
continue;
}
T r = mat[row][col];
for (int nxt = col; nxt < sz; nxt++) {
mat[row][nxt] -= r * mat[ran][nxt];
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[row][nxt] -= r * ops[ran][nxt];
}
}
pc.push_back(col);
ran++;
}
if (ran == sz) {
// Solve U * A^{-1} = E; U itself need not be reduced above its pivots.
for (int col = sz - 1; col >= 0; --col) {
for (int row = 0; row < col; ++row) {
T r = mat[row][col];
if (r == T{}) continue;
for (int nxt = 0; nxt < sz; ++nxt) ops[row][nxt] -= r * ops[col][nxt];
}
}
for (auto &row : ops) for (auto &value : row) value *= det;
return ops;
}
std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
if (ran + 1 != sz) {
return res;
}
std::vector<bool> ip(sz);
for (int col : pc) {
ip[col] = true;
}
int fc = 0;
while (ip[fc]) {
fc++;
}
std::vector<T> rv(sz, T{});
rv[fc] = T(1);
for (int row = ran - 1; row >= 0; --row) {
for (int col = pc[row] + 1; col < sz; ++col) {
rv[pc[row]] -= mat[row][col] * rv[col];
}
}
T scl = det;
if ((ran + fc) % 2 != 0) {
scl = -scl;
}
for (int row = 0; row < sz; row++) {
for (int col = 0; col < sz; col++) {
res[row][col] = scl * rv[row] * ops[ran][col];
}
}
return res;
}
} // namespace noya
#endif // NOYA_ADJUGATE_MATRIX_HPP
/// @complexity Time: O(n^3).
/// Space: O(n^2).
#include <algorithm>
#include <cassert>
#include <vector>
namespace noya {
/// @brief Compute the adjugate of a square matrix over a field.
/// Forward elimination tracks E A = U. Full rank uses back substitution for
/// A^{-1}; rank n-1 uses a right-null vector and E's last (left-null) row.
/// Rank at most n-2 makes every (n-1)-minor zero.
template <class T>
std::vector<std::vector<T>> adjugate_matrix(std::vector<std::vector<T>> mat) {
int sz = int(mat.size());
for (const auto &row : mat) {
assert(int(row.size()) == sz);
}
std::vector<std::vector<T>> ops(sz, std::vector<T>(sz, T{}));
for (int row = 0; row < sz; row++) {
ops[row][row] = T(1);
}
T det = T(1);
std::vector<int> pc;
int ran = 0;
for (int col = 0; col < sz && ran < sz; col++) {
int piv = ran;
while (piv < sz && mat[piv][col] == T{}) {
piv++;
}
if (piv == sz) {
continue;
}
if (piv != ran) {
std::swap(mat[piv], mat[ran]);
std::swap(ops[piv], ops[ran]);
det = -det;
}
T pv = mat[ran][col];
det *= pv;
T inv = T(1) / pv;
for (int nxt = col; nxt < sz; nxt++) {
mat[ran][nxt] *= inv;
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[ran][nxt] *= inv;
}
for (int row = ran + 1; row < sz; row++) {
if (mat[row][col] == T{}) {
continue;
}
T r = mat[row][col];
for (int nxt = col; nxt < sz; nxt++) {
mat[row][nxt] -= r * mat[ran][nxt];
}
for (int nxt = 0; nxt < sz; nxt++) {
ops[row][nxt] -= r * ops[ran][nxt];
}
}
pc.push_back(col);
ran++;
}
if (ran == sz) {
// Solve U * A^{-1} = E; U itself need not be reduced above its pivots.
for (int col = sz - 1; col >= 0; --col) {
for (int row = 0; row < col; ++row) {
T r = mat[row][col];
if (r == T{}) continue;
for (int nxt = 0; nxt < sz; ++nxt) ops[row][nxt] -= r * ops[col][nxt];
}
}
for (auto &row : ops) for (auto &value : row) value *= det;
return ops;
}
std::vector<std::vector<T>> res(sz, std::vector<T>(sz, T{}));
if (ran + 1 != sz) {
return res;
}
std::vector<bool> ip(sz);
for (int col : pc) {
ip[col] = true;
}
int fc = 0;
while (ip[fc]) {
fc++;
}
std::vector<T> rv(sz, T{});
rv[fc] = T(1);
for (int row = ran - 1; row >= 0; --row) {
for (int col = pc[row] + 1; col < sz; ++col) {
rv[pc[row]] -= mat[row][col] * rv[col];
}
}
T scl = det;
if ((ran + fc) % 2 != 0) {
scl = -scl;
}
for (int row = 0; row < sz; row++) {
for (int col = 0; col < sz; col++) {
res[row][col] = scl * rv[row] * ops[ran][col];
}
}
return res;
}
} // namespace noya
// NOYA_ADJUGATE_MATRIX_HPP