characteristic_polynomial.hpp

SECTIONMath INCLUDEnoya/characteristic_polynomial.hpp

计算 \(\det(xI-A)\) 的全部系数；用于矩阵特征值关系、递推和矩阵函数。

\[ \displaystyle \chi_A(x)=\det(xI-A)=\sum_{i=0}^{n} c_i x^i \]

Complexity: Time: O(n^3) field operations. Space: O(n^2).

AC 记录：characteristic_polynomial。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

当前头文件，省略 include guard；依赖见 #include。

/// @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>> mat) {
  int n = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == n);
  }
  for (int col = 0; col + 2 < n; col++) {
    int piv = col + 1;
    while (piv < n && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == n) {
      continue;
    }
    if (piv != col + 1) {
      std::swap(mat[piv], mat[col + 1]);
      for (int row = 0; row < n; row++) {
        std::swap(mat[row][piv], mat[row][col + 1]);
      }
    }
    T inv = T(1) / mat[col + 1][col];
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) {
        continue;
      }
      for (int nxt = col; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
      for (int prv = 0; prv < n; prv++) {
        mat[prv][col + 1] += r * mat[prv][row];
      }
    }
  }

  std::vector<std::vector<T>> pre(n + 1);
  pre[0] = {T(1)};
  for (int sz = 1; sz <= n; sz++) {
    pre[sz].assign(sz + 1, T{});
    for (int deg = 0; deg < sz; deg++) {
      pre[sz][deg + 1] += pre[sz - 1][deg];
      pre[sz][deg] -= mat[sz - 1][sz - 1] * pre[sz - 1][deg];
    }
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      T fct = mat[st - 1][sz - 1] * prd;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return pre[n];
}

} // namespace noya

#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>> mat) {
  int n = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == n);
  }
  for (int col = 0; col + 2 < n; col++) {
    int piv = col + 1;
    while (piv < n && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == n) {
      continue;
    }
    if (piv != col + 1) {
      std::swap(mat[piv], mat[col + 1]);
      for (int row = 0; row < n; row++) {
        std::swap(mat[row][piv], mat[row][col + 1]);
      }
    }
    T inv = T(1) / mat[col + 1][col];
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) {
        continue;
      }
      for (int nxt = col; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
      for (int prv = 0; prv < n; prv++) {
        mat[prv][col + 1] += r * mat[prv][row];
      }
    }
  }

  std::vector<std::vector<T>> pre(n + 1);
  pre[0] = {T(1)};
  for (int sz = 1; sz <= n; sz++) {
    pre[sz].assign(sz + 1, T{});
    for (int deg = 0; deg < sz; deg++) {
      pre[sz][deg + 1] += pre[sz - 1][deg];
      pre[sz][deg] -= mat[sz - 1][sz - 1] * pre[sz - 1][deg];
    }
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      T fct = mat[st - 1][sz - 1] * prd;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return pre[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>> mat) {
  int n = int(mat.size());
  for (const auto &row : mat) {
    assert(int(row.size()) == n);
  }
  for (int col = 0; col + 2 < n; col++) {
    int piv = col + 1;
    while (piv < n && mat[piv][col] == T{}) {
      piv++;
    }
    if (piv == n) {
      continue;
    }
    if (piv != col + 1) {
      std::swap(mat[piv], mat[col + 1]);
      for (int row = 0; row < n; row++) {
        std::swap(mat[row][piv], mat[row][col + 1]);
      }
    }
    T inv = T(1) / mat[col + 1][col];
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) {
        continue;
      }
      for (int nxt = col; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
      for (int prv = 0; prv < n; prv++) {
        mat[prv][col + 1] += r * mat[prv][row];
      }
    }
  }

  std::vector<std::vector<T>> pre(n + 1);
  pre[0] = {T(1)};
  for (int sz = 1; sz <= n; sz++) {
    pre[sz].assign(sz + 1, T{});
    for (int deg = 0; deg < sz; deg++) {
      pre[sz][deg + 1] += pre[sz - 1][deg];
      pre[sz][deg] -= mat[sz - 1][sz - 1] * pre[sz - 1][deg];
    }
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      T fct = mat[st - 1][sz - 1] * prd;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return pre[n];
}

} // namespace noya