另一種特征多項式算法

在 OI 中，比較普及的求解矩陣的特征多項式的算法是這個，在閱讀一些文獻后，這里給出另一種可實現的做法，不過從實測結果來看不是很有優勢。

對於給定的矩陣 \(A\) 和向量 \(v\)，我們設 \(p\) 是最大的正整數使得 \(\{v,Av,\dots,A^{p-1}v\}\) 線性無關。事實上，\(\mathcal K_r(A,v)=\operatorname{span}\{v,Ap,\dots,A^{p-1}v\}\) 被稱為 Krylov 子空間。

如果 \(p=n\)，我們設 \(K = \begin{pmatrix} v & Av & A^2v & \cdots & A^{p-1}v \end{pmatrix}\)，那么就有 \(AK=\begin{pmatrix} Av & A^2v & A^3v & \cdots & A^pv \end{pmatrix}\)。那么進一步就有

\[K^{-1}AK = \begin{pmatrix} & & & & *\\ 1 & & & & *\\ & 1 & & & \vdots\\ & & \ddots & & \vdots\\ & & & 1 & * \end{pmatrix}\]

這種矩陣的特征多項式是很顯然的。但是另一方面也很顯然，對於任意的 \(A\)，這樣的向量肯定是不一定存在的，比如 \(A\) 的最小多項式不是特征多項式的時候。那我們這個時候怎么辦呢？

我們令 \(V\) 為當前所有向量張成的線性空間，選取 \(V\) 外一個向量 \(u\)，繼續嘗試加入 \(u, Au,\dots\)，直到又得到一個在原來空間內的向量，這時我們再次選一個新的 \(V\) 以外的向量繼續做……最后，我們得到一個形如這樣一個線性無關的向量序列：

\[v_1, Av_1,\dots,A^{p_1-1}v_1, v_2, \dots, A^{p_2-1}v_2, \dots, v_k,\dots,A^{p_k-1}v_k \]

它同時還滿足 \(A^{p_i}v_i\) 能被前面的向量表出。我們這些向量並排得到矩陣 \(K\)，那么這時，\(K^{-1}AK\) 就形如下圖這樣：

顯然，整個矩陣的特征多項式也就是對角線上每個塊的特征多項式之乘積了。

具體來說，我們有如下流程：

1. 維護一個線性無關組 \(S\)，初始為空。特征多項式 \(p(t)\) 初始為 \(1\)。
2. 按照 \(i=1,\dots,n\) 的順序考慮單位向量 \(e_i\)，設為 \(v\)，不斷給 \(S\) 插入 \(v,Av,\dots\)，直到得到一個 \(A^pv\) 已經能被 \(S\) 線性表出，通過高斯消元得到系數 \(A^pv - c_1 A^{p-1}v - \cdots - c_p v \in \operatorname{span}\{\mathrm{previous \, vectors}\}\)。此時將 \(p(t)\) 乘以 \(t^p - c_1 t^{p-1} - \cdots - c_p\)。

代碼實現

#include <bits/stdc++.h>

using namespace std;

using ull = unsigned long long;

const int P = 998244353;

int norm(int x) { return x >= P ? x - P : x; }
int reduce(int x) { return x < 0 ? x + P : x; }
int neg(int x) { return x ? P - x : 0; }
void add(int& x, int y) { if ((x += y - P) < 0) x += P; }
void sub(int& x, int y) { if ((x -= y) < 0) x += P; }
void fam(int& x, int y, int z) { x = (x + y * (ull)z) % P; }
int mpow(int x, unsigned k) {
	if (k == 0) return 1;
	int ret = mpow(x * (ull)x % P, k >> 1);
	if (k & 1) ret = ret * (ull)x % P;
	return ret;
}
int quo2(int x) { return ((x & 1) ? x + P : x) >> 1; }

const int _ = 510;

int mat[_][_], bas[_][_], coe[_][_];
int vec[_], elim[_], tmp[_], pol[_];

int main() {
	ios::sync_with_stdio(false); cin.tie(NULL);

	int N; cin >> N;
	for (int i = 1; i <= N; ++i) for (int j = 1; j <= N; ++j) cin >> mat[i][j];

	pol[0] = 1;
	int dim = 0;
	for (int i = 1; i <= N; ++i) {
		memset(vec, 0, sizeof(vec)); vec[i] = 1;
		int lst = dim;
		while (true) {
			int id = dim + 1;
			memset(tmp, 0, sizeof(tmp));
			tmp[id] = 1;
			memcpy(elim, vec, sizeof(elim));
			bool fl = false;
			for (int j = 1; j <= N; ++j)
				if (elim[j]) {
					if (bas[j][j]) {
						int c = reduce(-elim[j]);
						for (int k = 1; k <= N; ++k) fam(elim[k], c, bas[j][k]);
						for (int k = 1; k <= N; ++k) fam(tmp[k], c, coe[j][k]);
					} else {
						fl = true;
						int c = mpow(elim[j], P - 2);
						for (int k = 1; k <= N; ++k) elim[k] = elim[k] * (ull)c % P;
						for (int k = 1; k <= N; ++k) tmp[k] = tmp[k] * (ull)c % P;
						memcpy(bas[j], elim, sizeof(elim));
						memcpy(coe[j], tmp, sizeof(tmp));
						break;
					}
				}
			if (!fl) {
				for (int j = dim; j; --j) for (int k = 1; k <= min(dim - lst, j); ++k)
					fam(pol[j], pol[j - k], tmp[dim + 1 - k]);
				break;
			}
			for (int j = 1; j <= N; ++j) {
				tmp[j] = 0;
				for (int k = 1; k <= N; ++k) fam(tmp[j], mat[j][k], vec[k]);
			}
			memcpy(vec, tmp, sizeof(tmp));
			++dim;
		}
		if (dim == N) break;
	}
	reverse(pol, pol + N + 1);
	for (int i = 0; i <= N; ++i) cout << pol[i] << " \n"[i == N];

	return 0;
}