characteristic_polynomial.hpp

SECTIONMath INCLUDEnoya/math/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, \(O(n^{2})\) for triangular inputs. Space: \(O(n^{2})\).

AC 记录：characteristic_polynomial。

跳到代码 · GitHub ↗

Implementation

ReferenceSourceExpanded

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

/// @complexity Time: O(n^3) field operations, O(n^2) for triangular inputs.
/// 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);
  }
  // Triangular matrices have characteristic polynomial prod(x - a_ii).
  bool upper = true, lower = true;
  for (int row = 0; row < n && (upper || lower); row++) {
    for (int col = row + 1; col < n && (upper || lower); col++) {
      upper = upper && mat[col][row] == T{};
      lower = lower && mat[row][col] == T{};
    }
  }
  if (upper || lower) {
    std::vector<T> result(n + 1);
    result[0] = T(1);
    for (int row = 0; row < n; row++) {
      T negative{};
      negative -= mat[row][row];
      result[row + 1] = result[row];
      for (int deg = row; deg > 0; deg--) {
        result[deg] *= negative;
        result[deg] += result[deg - 1];
      }
      result[0] *= negative;
    }
    return result;
  }
  std::vector<std::pair<int, T>> ratios;
  ratios.reserve(n > 2 ? n : 0);
  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];
    // Group the commuting row eliminations before the matching column update.
    ratios.clear();
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) continue;
      ratios.emplace_back(row, r);
      mat[row][col] = T{};
      for (int nxt = col + 1; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
    }
    // Visit each row consecutively instead of repeatedly walking full columns.
    if (!ratios.empty()) {
      for (int prv = 0; prv < n; prv++) {
        T delta{};
        for (const auto &[row, r] : ratios) delta += r * mat[prv][row];
        mat[prv][col + 1] += delta;
      }
    }
  }

  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{});
    const auto &last = pre[sz - 1];
    T diagonal = mat[sz - 1][sz - 1];
    pre[sz][0] -= diagonal * last[0];
    for (int deg = 1; deg < sz; deg++) {
      pre[sz][deg] = last[deg - 1];
      pre[sz][deg] -= diagonal * last[deg];
    }
    pre[sz][sz] = T(1);
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      if (prd == T{}) break; // Every earlier term contains this same zero factor.
      T fct = mat[st - 1][sz - 1] * prd;
      if (fct == T{}) continue;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return std::move(pre[n]);
}

} // namespace noya

 // NOYA_CHARACTERISTIC_POLYNOMIAL_HPP

#ifndef NOYA_CHARACTERISTIC_POLYNOMIAL_HPP
#define NOYA_CHARACTERISTIC_POLYNOMIAL_HPP 1

/// @complexity Time: O(n^3) field operations, O(n^2) for triangular inputs.
/// 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);
  }
  // Triangular matrices have characteristic polynomial prod(x - a_ii).
  bool upper = true, lower = true;
  for (int row = 0; row < n && (upper || lower); row++) {
    for (int col = row + 1; col < n && (upper || lower); col++) {
      upper = upper && mat[col][row] == T{};
      lower = lower && mat[row][col] == T{};
    }
  }
  if (upper || lower) {
    std::vector<T> result(n + 1);
    result[0] = T(1);
    for (int row = 0; row < n; row++) {
      T negative{};
      negative -= mat[row][row];
      result[row + 1] = result[row];
      for (int deg = row; deg > 0; deg--) {
        result[deg] *= negative;
        result[deg] += result[deg - 1];
      }
      result[0] *= negative;
    }
    return result;
  }
  std::vector<std::pair<int, T>> ratios;
  ratios.reserve(n > 2 ? n : 0);
  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];
    // Group the commuting row eliminations before the matching column update.
    ratios.clear();
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) continue;
      ratios.emplace_back(row, r);
      mat[row][col] = T{};
      for (int nxt = col + 1; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
    }
    // Visit each row consecutively instead of repeatedly walking full columns.
    if (!ratios.empty()) {
      for (int prv = 0; prv < n; prv++) {
        T delta{};
        for (const auto &[row, r] : ratios) delta += r * mat[prv][row];
        mat[prv][col + 1] += delta;
      }
    }
  }

  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{});
    const auto &last = pre[sz - 1];
    T diagonal = mat[sz - 1][sz - 1];
    pre[sz][0] -= diagonal * last[0];
    for (int deg = 1; deg < sz; deg++) {
      pre[sz][deg] = last[deg - 1];
      pre[sz][deg] -= diagonal * last[deg];
    }
    pre[sz][sz] = T(1);
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      if (prd == T{}) break; // Every earlier term contains this same zero factor.
      T fct = mat[st - 1][sz - 1] * prd;
      if (fct == T{}) continue;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return std::move(pre[n]);
}

} // namespace noya

#endif // NOYA_CHARACTERISTIC_POLYNOMIAL_HPP

/// @complexity Time: O(n^3) field operations, O(n^2) for triangular inputs.
/// 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);
  }
  // Triangular matrices have characteristic polynomial prod(x - a_ii).
  bool upper = true, lower = true;
  for (int row = 0; row < n && (upper || lower); row++) {
    for (int col = row + 1; col < n && (upper || lower); col++) {
      upper = upper && mat[col][row] == T{};
      lower = lower && mat[row][col] == T{};
    }
  }
  if (upper || lower) {
    std::vector<T> result(n + 1);
    result[0] = T(1);
    for (int row = 0; row < n; row++) {
      T negative{};
      negative -= mat[row][row];
      result[row + 1] = result[row];
      for (int deg = row; deg > 0; deg--) {
        result[deg] *= negative;
        result[deg] += result[deg - 1];
      }
      result[0] *= negative;
    }
    return result;
  }
  std::vector<std::pair<int, T>> ratios;
  ratios.reserve(n > 2 ? n : 0);
  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];
    // Group the commuting row eliminations before the matching column update.
    ratios.clear();
    for (int row = col + 2; row < n; row++) {
      T r = mat[row][col] * inv;
      if (r == T{}) continue;
      ratios.emplace_back(row, r);
      mat[row][col] = T{};
      for (int nxt = col + 1; nxt < n; nxt++) {
        mat[row][nxt] -= r * mat[col + 1][nxt];
      }
    }
    // Visit each row consecutively instead of repeatedly walking full columns.
    if (!ratios.empty()) {
      for (int prv = 0; prv < n; prv++) {
        T delta{};
        for (const auto &[row, r] : ratios) delta += r * mat[prv][row];
        mat[prv][col + 1] += delta;
      }
    }
  }

  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{});
    const auto &last = pre[sz - 1];
    T diagonal = mat[sz - 1][sz - 1];
    pre[sz][0] -= diagonal * last[0];
    for (int deg = 1; deg < sz; deg++) {
      pre[sz][deg] = last[deg - 1];
      pre[sz][deg] -= diagonal * last[deg];
    }
    pre[sz][sz] = T(1);
    T prd = T(1);
    for (int st = sz - 1; st >= 1; st--) {
      prd *= mat[st][st - 1];
      if (prd == T{}) break; // Every earlier term contains this same zero factor.
      T fct = mat[st - 1][sz - 1] * prd;
      if (fct == T{}) continue;
      for (int deg = 0; deg < int(pre[st - 1].size()); deg++) {
        pre[sz][deg] -= fct * pre[st - 1][deg];
      }
    }
  }
  return std::move(pre[n]);
}

} // namespace noya

 // NOYA_CHARACTERISTIC_POLYNOMIAL_HPP