characteristic_polynomial.hpp¶
Return coefficients c_i of det(xI-A)=sum c_i x^i over a field in O(n^3), using a similarity reduction to upper Hessenberg form.
Verified by characteristic_polynomial.
\[
\displaystyle \chi_A(x)=\det(xI-A)=\sum_{i=0}^{n} c_i x^i
\]
Implementation¶
#ifndef NOYA_CHARACTERISTIC_POLYNOMIAL_HPP
#define NOYA_CHARACTERISTIC_POLYNOMIAL_HPP 1
/// @complexity Time: O(n^3) field operations.
/// Space: O(n^2).
#include <cassert>
#include <utility>
#include <vector>
namespace noya {
/// @brief Return coefficients c_i of det(xI-A)=sum c_i x^i over a field in
/// O(n^3), using a similarity reduction to upper Hessenberg form.
template <class T>
std::vector<T>
characteristic_polynomial(std::vector<std::vector<T>> matrix) {
int n = int(matrix.size());
for (const auto &row : matrix) {
assert(int(row.size()) == n);
}
for (int column = 0; column + 2 < n; column++) {
int pivot = column + 1;
while (pivot < n && matrix[pivot][column] == T{}) {
pivot++;
}
if (pivot == n) {
continue;
}
if (pivot != column + 1) {
std::swap(matrix[pivot], matrix[column + 1]);
for (int row = 0; row < n; row++) {
std::swap(matrix[row][pivot], matrix[row][column + 1]);
}
}
T inverse = T(1) / matrix[column + 1][column];
for (int row = column + 2; row < n; row++) {
T ratio = matrix[row][column] * inverse;
if (ratio == T{}) {
continue;
}
for (int next = column; next < n; next++) {
matrix[row][next] -= ratio * matrix[column + 1][next];
}
for (int previous = 0; previous < n; previous++) {
matrix[previous][column + 1] += ratio * matrix[previous][row];
}
}
}
std::vector<std::vector<T>> prefix(n + 1);
prefix[0] = {T(1)};
for (int size = 1; size <= n; size++) {
prefix[size].assign(size + 1, T{});
for (int degree = 0; degree < size; degree++) {
prefix[size][degree + 1] += prefix[size - 1][degree];
prefix[size][degree] -=
matrix[size - 1][size - 1] * prefix[size - 1][degree];
}
T product = T(1);
for (int start = size - 1; start >= 1; start--) {
product *= matrix[start][start - 1];
T factor = matrix[start - 1][size - 1] * product;
for (int degree = 0; degree < int(prefix[start - 1].size()); degree++) {
prefix[size][degree] -= factor * prefix[start - 1][degree];
}
}
}
return prefix[n];
}
} // namespace noya
#endif // NOYA_CHARACTERISTIC_POLYNOMIAL_HPP
#include <cassert>
#include <utility>
#include <vector>
/// @complexity Time: O(n^3) field operations.
/// Space: O(n^2).
namespace noya {
/// @brief Return coefficients c_i of det(xI-A)=sum c_i x^i over a field in
/// O(n^3), using a similarity reduction to upper Hessenberg form.
template <class T>
std::vector<T>
characteristic_polynomial(std::vector<std::vector<T>> matrix) {
int n = int(matrix.size());
for (const auto &row : matrix) {
assert(int(row.size()) == n);
}
for (int column = 0; column + 2 < n; column++) {
int pivot = column + 1;
while (pivot < n && matrix[pivot][column] == T{}) {
pivot++;
}
if (pivot == n) {
continue;
}
if (pivot != column + 1) {
std::swap(matrix[pivot], matrix[column + 1]);
for (int row = 0; row < n; row++) {
std::swap(matrix[row][pivot], matrix[row][column + 1]);
}
}
T inverse = T(1) / matrix[column + 1][column];
for (int row = column + 2; row < n; row++) {
T ratio = matrix[row][column] * inverse;
if (ratio == T{}) {
continue;
}
for (int next = column; next < n; next++) {
matrix[row][next] -= ratio * matrix[column + 1][next];
}
for (int previous = 0; previous < n; previous++) {
matrix[previous][column + 1] += ratio * matrix[previous][row];
}
}
}
std::vector<std::vector<T>> prefix(n + 1);
prefix[0] = {T(1)};
for (int size = 1; size <= n; size++) {
prefix[size].assign(size + 1, T{});
for (int degree = 0; degree < size; degree++) {
prefix[size][degree + 1] += prefix[size - 1][degree];
prefix[size][degree] -=
matrix[size - 1][size - 1] * prefix[size - 1][degree];
}
T product = T(1);
for (int start = size - 1; start >= 1; start--) {
product *= matrix[start][start - 1];
T factor = matrix[start - 1][size - 1] * product;
for (int degree = 0; degree < int(prefix[start - 1].size()); degree++) {
prefix[size][degree] -= factor * prefix[start - 1][degree];
}
}
}
return prefix[n];
}
} // namespace noya