Upper-Hessenberg 矩阵

\[ H_n = \begin{bmatrix} a_{0, 0} & a_{0, 1} & a_{0, 2} & a_{0, 3} & \cdots & a_{0, n - 1} & a_{0, n}\\ a_{1, 0} & a_{1, 1} & a_{1, 2} & a_{1, 3} & \cdots & a_{1, n - 1} & a_{1, n}\\ 0 & a_{2, 1} & a_{2, 2} & a_{2, 3} & \cdots & a_{2, n - 1} & a_{2, n}\\ 0 & 0 & a_{3, 2} & a_{3, 3} & \cdots & a_{3, n - 1} & a_{3, n}\\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & a_{n - 1, n-2} & a_{n - 1, n - 1} & a_{n - 1, n}\\ 0 & 0 & 0 & 0 & 0 & a_{n, n - 1} & a_{n, n} \end{bmatrix} \]

设 \(h_n = \det H_n\)​，则 \[ \begin{aligned} h_n &= a_{n, n}h_{n-1} - a_{n, n-1} \begin{vmatrix} a_{0, 0} & a_{0, 1} & a_{0, 2} & a_{0, 3} & \cdots & a_{0, n}\\ a_{1, 0} & a_{1, 1} & a_{1, 2} & a_{1, 3} & \cdots & a_{1, n}\\ 0 & a_{2, 1} & a_{2, 2} & a_{2, 3} & \cdots & a_{2, n}\\ 0 & 0 & a_{3, 2} & a_{3, 3} & \cdots & a_{3, n}\\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & 0 & 0 & 0 & a_{n - 1, n-2} & a_{n - 1, n}\\ \end{vmatrix}\\ &= a_{n, n}h_{n-1} - a_{n, n-1}\left( a_{n-1, n}h_{n-2} - a_{n-1, n-2} \begin{vmatrix} a_{0, 0} & a_{0, 1} & a_{0, 2} & \cdots & a_{0, n}\\ a_{1, 0} & a_{1, 1} & a_{1, 2} & \cdots & a_{1, n}\\ 0 & a_{2, 1} & a_{2, 2} & \cdots & a_{2, n}\\ 0 & 0 & a_{3, 2} & \cdots & a_{3, n}\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & 0 & a_{n - 2, n-3} & a_{n - 2, n}\\ \end{vmatrix}\\ \right)\\ &= \cdots\cdots\cdots\cdots\\ &= \sum_{i = 0}^{n-1}h_{n-1-i}\cdot(-1)^ia_{n-i, n}\prod_{j=0}^{i-1}a_{n-j, n-j-1} \end{aligned} \]

复杂度 \(O(n^3)\)。

std::vector<Z> charPoly(std::vector<std::vector<Z>> a) {
  int n = a.size();
  for (int j = 0; j < n - 2; j++) {
    for (int i = j + 1; i < n; i++) {
      if (!a[i][j]) continue;
      std::swap(a[i], a[j + 1]);
      for (int k = 0; k < n; k++) {
        std::swap(a[k][i], a[k][j + 1]);
      }
      break;
    }
    if (a[j + 1][j]) {
      auto inv = a[j + 1][j].inv();
      for (int i = j + 2; i < n; i++) {
        auto c = a[i][j] * inv;
        for (int k = 0; k < n; k++) a[i][k] -= a[j + 1][k] * c;
        for (int k = 0; k < n; k++) a[k][j + 1] += a[k][i] * c;
      }
    }
  }
  std::vector<std::vector<Z>> h(n + 1);
  h[0] = {1};
  for (int i = 0; i < n; i++) {
    Z prod = 1;
    h[i + 1].resize(h[i].size() + 1);
    for (int j = 0; j < h[i].size(); j++) {
      h[i + 1][j + 1] = h[i][j];
      h[i + 1][j] -= h[i][j] * a[i][i];
    }
    for (int j = 0; j < i; j++) {
      prod *= a[i - j][i - j - 1];
      auto c = a[i - j - 1][i] * prod;
      for (int k = 0; k < h[i - j - 1].size(); k++) {
        h[i + 1][k] -= h[i - j - 1][k] * c;
      }
    }
  }
  return h[n];
}

最后修改于 2022-03-08