adjugate_matrix.hpp¶
Compute the adjugate of a square matrix over a field. Gauss--Jordan elimination tracks the row-operation matrix E with E A = R. For full rank, E = A^{-1}. For rank n-1, adj(R) is an outer product of the unique right-null vector and the final zero-row selector; multiplying by E recovers adj(A). Rank at most n-2 makes every (n-1)-minor zero.
Verified by adjugate_matrix.
\[
\displaystyle \operatorname{adj}(A)A=\det(A)I
\]
Implementation¶
#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. Gauss--Jordan
/// elimination tracks the row-operation matrix E with `E A = R`. For full
/// rank, `E = A^{-1}`. For rank n-1, `adj(R)` is an outer product of the
/// unique right-null vector and the final zero-row selector; multiplying by E
/// recovers `adj(A)`. 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>> matrix) {
int size = int(matrix.size());
for (const auto &row : matrix) {
assert(int(row.size()) == size);
}
std::vector<std::vector<T>> operations(size,
std::vector<T>(size, T{}));
for (int row = 0; row < size; row++) {
operations[row][row] = T(1);
}
T determinant = T(1);
std::vector<int> pivot_columns;
int rank = 0;
for (int column = 0; column < size && rank < size; column++) {
int pivot = rank;
while (pivot < size && matrix[pivot][column] == T{}) {
pivot++;
}
if (pivot == size) {
continue;
}
if (pivot != rank) {
std::swap(matrix[pivot], matrix[rank]);
std::swap(operations[pivot], operations[rank]);
determinant = -determinant;
}
T pivot_value = matrix[rank][column];
determinant *= pivot_value;
T inverse = T(1) / pivot_value;
for (int next = column; next < size; next++) {
matrix[rank][next] *= inverse;
}
for (int next = 0; next < size; next++) {
operations[rank][next] *= inverse;
}
for (int row = 0; row < size; row++) {
if (row == rank || matrix[row][column] == T{}) {
continue;
}
T ratio = matrix[row][column];
for (int next = column; next < size; next++) {
matrix[row][next] -= ratio * matrix[rank][next];
}
for (int next = 0; next < size; next++) {
operations[row][next] -= ratio * operations[rank][next];
}
}
pivot_columns.push_back(column);
rank++;
}
std::vector<std::vector<T>> result(size,
std::vector<T>(size, T{}));
if (rank == size) {
for (int row = 0; row < size; row++) {
for (int column = 0; column < size; column++) {
result[row][column] = determinant * operations[row][column];
}
}
return result;
}
if (rank + 1 != size) {
return result;
}
std::vector<bool> is_pivot(size);
for (int column : pivot_columns) {
is_pivot[column] = true;
}
int free_column = 0;
while (is_pivot[free_column]) {
free_column++;
}
std::vector<T> right_null(size, T{});
right_null[free_column] = T(1);
for (int row = 0; row < rank; row++) {
right_null[pivot_columns[row]] = -matrix[row][free_column];
}
T scale = determinant;
if ((rank + free_column) % 2 != 0) {
scale = -scale;
}
for (int row = 0; row < size; row++) {
for (int column = 0; column < size; column++) {
result[row][column] =
scale * right_null[row] * operations[rank][column];
}
}
return result;
}
} // namespace noya
#endif // NOYA_ADJUGATE_MATRIX_HPP
#include <algorithm>
#include <cassert>
#include <vector>
/// @complexity Time: O(n^3).
/// Space: O(n^2).
namespace noya {
/// @brief Compute the adjugate of a square matrix over a field. Gauss--Jordan
/// elimination tracks the row-operation matrix E with `E A = R`. For full
/// rank, `E = A^{-1}`. For rank n-1, `adj(R)` is an outer product of the
/// unique right-null vector and the final zero-row selector; multiplying by E
/// recovers `adj(A)`. 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>> matrix) {
int size = int(matrix.size());
for (const auto &row : matrix) {
assert(int(row.size()) == size);
}
std::vector<std::vector<T>> operations(size,
std::vector<T>(size, T{}));
for (int row = 0; row < size; row++) {
operations[row][row] = T(1);
}
T determinant = T(1);
std::vector<int> pivot_columns;
int rank = 0;
for (int column = 0; column < size && rank < size; column++) {
int pivot = rank;
while (pivot < size && matrix[pivot][column] == T{}) {
pivot++;
}
if (pivot == size) {
continue;
}
if (pivot != rank) {
std::swap(matrix[pivot], matrix[rank]);
std::swap(operations[pivot], operations[rank]);
determinant = -determinant;
}
T pivot_value = matrix[rank][column];
determinant *= pivot_value;
T inverse = T(1) / pivot_value;
for (int next = column; next < size; next++) {
matrix[rank][next] *= inverse;
}
for (int next = 0; next < size; next++) {
operations[rank][next] *= inverse;
}
for (int row = 0; row < size; row++) {
if (row == rank || matrix[row][column] == T{}) {
continue;
}
T ratio = matrix[row][column];
for (int next = column; next < size; next++) {
matrix[row][next] -= ratio * matrix[rank][next];
}
for (int next = 0; next < size; next++) {
operations[row][next] -= ratio * operations[rank][next];
}
}
pivot_columns.push_back(column);
rank++;
}
std::vector<std::vector<T>> result(size,
std::vector<T>(size, T{}));
if (rank == size) {
for (int row = 0; row < size; row++) {
for (int column = 0; column < size; column++) {
result[row][column] = determinant * operations[row][column];
}
}
return result;
}
if (rank + 1 != size) {
return result;
}
std::vector<bool> is_pivot(size);
for (int column : pivot_columns) {
is_pivot[column] = true;
}
int free_column = 0;
while (is_pivot[free_column]) {
free_column++;
}
std::vector<T> right_null(size, T{});
right_null[free_column] = T(1);
for (int row = 0; row < rank; row++) {
right_null[pivot_columns[row]] = -matrix[row][free_column];
}
T scale = determinant;
if ((rank + free_column) % 2 != 0) {
scale = -scale;
}
for (int row = 0; row < size; row++) {
for (int column = 0; column < size; column++) {
result[row][column] =
scale * right_null[row] * operations[rank][column];
}
}
return result;
}
} // namespace noya