特征多项式的 n^3 求法（注解）

https://oi-wiki.org/math/linear-algebra/char-poly/

OI-wiki 写的很好。这里只是一些注解。

相似矩阵

定理：\(A\) 的特征多项式与 \(PAP^{-1}\) 的相同。不证了。

使用高斯消元进行相似变换

要想将 \(A\) 相似变换成上 Hessenberg 矩阵（上海森堡矩阵），首先需要知道初等行变换对应的矩阵 \(P\) 的逆长什么样，以及它右乘 \(A\) 会使 \(A\) 变成什么，这样才能求出 \(PAP^{-1}\)。

交换两行

交换 \(p\) 行和 \(q\) 行，想象我们是在单位矩阵上做的，那么这个 \(P\) 应该就是：

• \(P[i][i]=1(i\neq p, q)\)
• \(P[p][q] = P[q][p]=1\)
• 其余位置为 \(0\)。

它的逆矩阵，根据定义，\(P^{-1}P=1\)，用求矩阵逆的做法（将矩阵消为单位阵所用的变换打在单位阵上，新的矩阵就是逆），那么将 \(P\) 还原为单位矩阵，只需要“交换 \(p\) 行和 \(q\) 行”即可，也就是说 \(P\) 的逆和 \(P\) 一样（或者交换列也可以得到一样的 \(P^{-1}\)）。如果 \(P^{-1}\) 要右乘 \(A\)，你可以自己模拟一下，拿 \(A\) 的一行的转置去扫 \(P\)，发现实际上就是交换了 \(p\) 列和 \(q\) 列。

结论：相似变换中，交换 \(p\) 行和 \(q\) 行后需要交换 \(p\) 列和 \(q\) 列。

一行的倍数加到另一行

\(p\) 行的 \(k\) 倍加到 \(q\) 行去，想象我们是在单位矩阵上做的，那么这个 \(P\) 应该就是：

• \(P[i][i]=1\)。
• \(P[q][p]=k\)。
• 其余位置为 \(0\)。

与上面的过程一样，其逆为：将 \(p\) 行的 \(-k\) 倍加到 \(q\) 行。右乘的效果为：将 \(q\) 列的 \(-k\) 倍加到 \(p\) 列。

结论：相似变换中，将 \(p\) 行的 \(-k\) 倍加到 \(q\) 行后需要将 \(q\) 列的 \(-k\) 倍加到 \(p\) 列。

高斯消元到底怎么消

枚举 \(i\)，消去第 \(i\) 列中次对角线下的所有元素，使用第 \(i\) 列次对角线上的元素与下面的元素进行对决。需要交换的是 \(\geq i+1\) 的行和列，不会影响到正在消元的第 \(i\) 列和已经完成的前 \(i-1\) 列。

Hessenberg 算法（海森堡算法）

目标是 \(O(n^3)\) 求出上海森堡矩阵的特征多项式。考虑记只考虑 \(A[1..i][1..i]\) 的元素的特征多项式为 \(p_i(x)\)。注意 \(p_0(x)=1\)。

从小到大枚举 \(i\)，去算特征多项式的定义 \(|xI-A|\)，考虑使用代数余子式技术展开最后一行：这里抄一个 \(n=3\) 的例子

\[\begin{aligned} p_3(x)&= \det(xI_3-H_3)\\ &=\begin{vmatrix} x-\alpha_1&-h_{12}&-h_{13}\\ -\beta_2&x-\alpha_2&-h_{23}\\ &-\beta_3&x-\alpha_3 \end{vmatrix}\\ &=(x-\alpha_3)\cdot (-1)^{3+3}p_2(x)-\beta_3\cdot (-1)^{3+2} \begin{vmatrix} x-\alpha_1&-h_{13}\\ -\beta_2&-h_{23} \end{vmatrix}\\ \end{aligned} \]

最下面就是两个代数余子式，就是说代数余子式一定是去选零多的行，例如每次都展开最后一行那么也就展开出两个代数余子式。展开完了以后，第一项可以计算，第二项需要递归下去，你想象一下 \(n=4\) 时你需要递归的是

\[\cdots+\beta_4\begin{vmatrix} x-\alpha_1&-h_{12}&-h_{14}\\ -\beta_2&x-\alpha_2&-h_{24}\\ &-\beta_3&-h_{34} \end{vmatrix}= -\beta_4h_{34}p_2(x)+\beta_4\beta_3\begin{vmatrix} x-\alpha_1&-h_{14}\\ -\beta_2&-h_{24} \end{vmatrix} \]

注意那个符号，主对角线上 \((-1)^{i+j}=1\)，次对角线自然就是 \(-1\)。这里已经初见端倪了，我们想象一下：

a  h  h  h  h  h  h
b  a  h  h  h  h  h
   b  a  h  h  h  h
      b  a  h  h  h
         b  a  h  h
            b  a  h
               b  a

我们这个递归就是从下往上走一段 \(\beta\)，然后突然往最右边一列 \(h\) 上撞，递归到已经算过的特征多项式，举例：

a  h  h  -  -  -  -         
b  a  h  -  -  -  -          
   b  a  -  -  -  -
      -  -  -  - [h]
        [b] -  -  -
           [b] -  -
              [b] -

中括号括起来的是我们到的需要取乘积的点，减号就是被划去的元素，剩下的就是 \(p_3(x)\)。

我们需要对次对角线的每个非空后缀都这么算一次，并求和（当然因为 \(h\) 是负的所以求和以后要减掉），加上划去右下角的那一项。

这时候再看这个式子就豁然开朗了：

\[p_i(x)=(x-\alpha_i)p_{i-1}(x)- \sum_{m=1}^{i-1}h_{i-m,i} \left( \prod_{j=i-m+1}^{i}\beta_j \right) p_{i-m-1}(x) \]

特征多项式优化矩阵快速幂

Cayley–Hamilton 定理（哈密顿-凯莱定理）：方阵 \(A\) 的特征多项式为 \(p(x)\)，则 \(p(A)=0\)。

想要求 \(A^k\)。先获取特征多项式 \(p(x)\)。令 \(f(x)=x^k\)，则我们要求的就是 \(f(A)\)，如果将 \(f(x)\bmod p(x)\) 这个新的多项式记作 \(h(x)\)，则我们说 \(f(A)=h(A)\)，这是在说什么呢？什么叫 \(\bmod 0\)，实际上是 \(f(x)=g(x)p(x)+h(x)\) 其中 \(g(x)\) 不管是什么反正是多项式，当代入 \(x=A\) 时 \(f(A)=g(A)p(A)+h(A)=h(A)\)。\(h(A)\) 的次数是 \(O(n)\) 的，这样就能得到优化了。

应用到常系数齐次线性递推就是：\(p(x)\) 可以提前手算，然后我们要的答案就是

\[A^kv=\sum_{i=0}^{n-1}h_iA^iv \]

\[(A^kv)_0=\sum_{i=0}^{n-1}h_i(A^iv)_0=\sum_{i=0}^{n-1}h_iv_i \]

\(O(n)\) 搞完，这下的工作就是求 \(x^k\bmod p(x)\)，这个的话，右边的东西快速幂，一边快速幂一遍对 \(p(x)\) 取模，根据取模算法的复杂度，有 \(O(n^2\log k)\) 或 \(O(n\log n\log k)\) 两种复杂度。

总结

你已经会了基于相似变换和海森堡算法的 \(O(n^3)\) 求矩阵特征多项式的做法，可以配合哈密顿-凯莱定理优化常系数齐次线性递推。求出特征多项式的所有根得到矩阵的所有特征值，然后就是另外一段关于特征值的故事了。

代码

【模板】特征多项式

不知道为什么跑的这么慢。

#include <bits/stdc++.h>
using namespace std;
#ifdef LOCAL
#define debug(...) fprintf(stderr, ##__VA_ARGS__)
#else
#define endl "\n"
#define debug(...) void(0)
#endif
using LL = long long;
template <class T>
using must_int = enable_if_t<is_integral<T>::value, int>;
template <unsigned umod>
struct modint {/*{{{*/
  static constexpr int mod = umod;
  unsigned v;
  modint() : v(0) {}
  template <class T, must_int<T> = 0>
  modint(T x) {
    x %= mod;
    v = x < 0 ? x + mod : x;
  }
  modint operator+() const { return *this; }
  modint operator-() const { return modint() - *this; }
  friend int raw(const modint &self) { return self.v; }
  friend ostream& operator<<(ostream& os, const modint &self) { 
    return os << raw(self); 
  }
  modint &operator+=(const modint &rhs) {
    v += rhs.v;
    if (v >= umod) v -= umod;
    return *this;
  }
  modint &operator-=(const modint &rhs) {
    v -= rhs.v;
    if (v >= umod) v += umod;
    return *this;
  }
  modint &operator*=(const modint &rhs) {
    v = 1ull * v * rhs.v % umod;
    return *this;
  }
  modint &operator/=(const modint &rhs) {
    assert(rhs.v);
    return *this *= qpow(rhs, mod - 2);
  }
  template <class T, must_int<T> = 0>
  friend modint qpow(modint a, T b) {
    modint r = 1;
    for (; b; b >>= 1, a *= a)
      if (b & 1) r *= a;
    return r;
  }
  friend modint operator+(modint lhs, const modint &rhs) { return lhs += rhs; }
  friend modint operator-(modint lhs, const modint &rhs) { return lhs -= rhs; }
  friend modint operator*(modint lhs, const modint &rhs) { return lhs *= rhs; }
  friend modint operator/(modint lhs, const modint &rhs) { return lhs /= rhs; }
  bool operator==(const modint &rhs) const { return v == rhs.v; }
  bool operator!=(const modint &rhs) const { return v != rhs.v; }
};/*}}}*/
using mint = modint<998244353>;
int n;
valarray<mint> a[510], p[510];
void guass() {
  for (int i = 0; i < n; i++) {
    for (int j = i + 2; j < n; j++) {
      while (raw(a[i + 1][i])) {
        mint k = raw(a[j][i]) / raw(a[i + 1][i]);
        a[j] -= a[i + 1] * k;
        for (int r = 0; r < n; r++) a[r][i + 1] += a[r][j] * k;
        swap(a[i + 1], a[j]);
        for (int r = 0; r < n; r++) swap(a[r][i + 1], a[r][j]);
      }
      swap(a[i + 1], a[j]);
      for (int r = 0; r < n; r++) swap(a[r][i + 1], a[r][j]);
    }
  }
}
valarray<mint> charpoly() {
  p[0].resize(n + 1);
  p[0][0] = -a[0][0], p[0][1] = 1;
  for (int i = 1; i < n; i++) {
    p[i] = p[i - 1].shift(-1) - a[i][i] * p[i - 1];
    mint pre = 1;
    for (int j = i - 1; j >= 0; j--) {
      pre *= a[j + 1][j];
      if (j) p[i] -= p[j - 1] * pre * a[j][i];
      else p[i][0] -= pre * a[j][i];
    }
    for (int j = 0; j <= n; j++) debug("%d%c", raw(p[i][j]), " \n"[j == n]);
  }
  return p[n - 1];
}
int main() {
#ifndef LOCAL
  cin.tie(nullptr)->sync_with_stdio(false);  
#endif
  cin >> n;
  for (int i = 0; i < n; i++) {
    a[i].resize(n);
    for (int j = 0; j < n; j++) cin >> a[i][j].v;
  }
  guass();
  for (int i = 0; i < n; i++) {
    for (int j = 0; j < n; j++) debug("%d%c", raw(a[i][j]), " \n"[j == n - 1]);
  }
  auto ret = charpoly();
  for (int i = 0; i <= n; i++) cout << ret[i] << " \n"[i == n];
  return 0;
}