Determinant of $A+Bz$

// created on 23.09.22

目录

• Determinant of \(A+Bz\)

Determinant of \(A+Bz\)

You are given two \(n\times n\) matrices \(A\) and \(B\) .
For each \(k=0,1,\cdots,n\) , you need to calculate \(f_k=[z^k]\det(A+Bz)\bmod 998244353\) .
\(n\leq 500\) 。

矩阵 \(A\) 的特征多项式为：\(f_A(\lambda)=|\lambda I-A|\) 。使得 \(f_A(\lambda)=0\) 的 \(\lambda\) 称为矩阵 \(A\) 的特征值。

先尝试解决矩阵的特征多项式问题。一个重要的性质是，特征多项式在基变更下不变：若存在可逆方阵 \(C\) 使得 \(B=P^{-1}AP\) 则 \(f_A(\lambda)=f_B(\lambda)\) 。

我们通过对高斯消元作出修改达成目的。一般的高斯消元可以看出，给出 \(A\) 求出 \(PA\)，后者是一个上三角矩阵。但此时差了一个 \(P^{-1}\) 。

注意到而高斯消元中，我们要做的只有交换两行，以及将一行乘 \(k\) 加到另一行。它们对应的初等矩阵的逆，在整个式子右侧乘上的结果是很好预测的（都是对列操作）。

于是，我们在 \(P\rightarrow xP\) 时同时执行 \(P^{-1}\rightarrow P^{-1}x^{-1}\)，并且实时维护 \(PAP^{-1}\) 。我们的目标变为消出上海森堡矩阵（即 \(j-i<-1\) 的 \((PAP^{-1})_{i,j}\) 都为 \(0\)，需要和上三角矩阵区分）。显然易见，这总是可以完成的。

接下来的任务就变成求上海森堡矩阵的特征多项式了（记 \(B=PAP^{-1}\)）。

记 \(p_i(\lambda)\) 为保留前 \(i\) 行前 \(i\) 列时，子方阵的特征多项式。\(i\rightarrow i+1\) 时，讨论 \(i+1\) 要么列的选择方式，可以得到：

\[p_i(\lambda)=(x-A_{i,i})p_{i-1}(x)-\sum_{j=1}^{i-1}A_{i-j,i}p_{i-j-1}(\lambda)\prod_{k=i-j+1}^{i}A_{k,k-1} \]

于是以 \(O(n^3)\) 的代价我们求出了 \(f_A(\lambda)=|\lambda I-A|\) 。

提交记录：Submission #186123 - QOJ.ac 。

现在的问题是解决 \(|A+Bz|\) 。假设 \(B\) 满秩，高斯消元消成 \(I=MB\)，那么 \(\det(A+Bz)=\det(M^{-1})\det(MA+MBz)\)，问题变为求特征多项式。否则，我们消元的过程中，如果消不下去了，就将这一列的 \(A\) 乘 \(z\) 搬过来（然后最后除掉这个 \(z\)）。如果仍然消不下去，那么 \(|A+Bz|\) 就是 \(0\) 。

提交记录：Submission #186128 - QOJ.ac 。

template <class T> inline vector <int> char_poly (int n, T _a) {
	static int a[N][N], p[N][N];
	lep (i, 1, n) lep (j, 1, n) a[i][j] = _a[i][j] ? mod - _a[i][j] : 0;
	lep (i, 1, n - 1) {
		int t = i + 1; for ( ; t <= n && ! a[t][i]; ++ t);
		if (t > n) continue ;
		if (t != i + 1) {
			lep (j, i, n) swap (a[t][j], a[i + 1][j]);
			lep (j, 1, n) swap (a[j][t], a[j][i + 1]);
		}
		int inv = qpow (a[i + 1][i], mod - 2);
		lep (j, i + 2, n) if (a[j][i]) {
			int buf = mul (a[j][i], inv);
			lep (k, i, n) sub (a[j][k], mul (a[i + 1][k], buf));
			lep (k, 1, n) pls (a[k][i + 1], mul (a[k][j], buf));
		}
	}
	p[n + 1][0] = 1;
	rep (i, n, 1) {
		lep (j, 0, n + 1 - i) p[i][j] = j ? p[i + 1][j - 1] : 0;
		for (int j = i, t = 1; j <= n; t = mul (t, a[j + 1][j]), ++ j) {
			int coef = mul (mul (t, a[i][j]), (((j - i) & 1) ? mod - 1 : 1));
			lep (k, 0, n - j) pls (p[i][k], mul (p[j + 1][k], coef));
		}
	}
	return vector <int> (p[1], p[1] + n + 1);
}
template <class T> inline vector <int> det (int n, T _a, T _b) {
	static int a[N][N], b[N][N];
	lep (i, 1, n) lep (j, 1, n) a[i][j] = _a[i][j], b[i][j] = _b[i][j];
	int offset = 0, ret = 1;
	lep (i, 1, n) {
		for (int j = 1, t; t = b[j][i], j < i; ++ j) {
			lep (k, 1, n) sub (a[k][i], mul (a[k][j], t));
			b[j][i] = 0;
		}
		int t = i; for ( ; t <= n && ! b[t][i]; ++ t);
		if (t > n) {
			++ offset;
			lep (j, 1, n) b[j][i] = a[j][i], a[j][i] = 0;
			for (int j = 1, t; t = b[j][i], j < i; ++ j) {
				lep (k, 1, n) sub (a[k][i], mul (a[k][j], t));
				b[j][i] = 0;
			}
			t = i; for ( ; t <= n && ! b[t][i]; ++ t);
			if (t > n) return vector <int> (n + 1, 0);
		}
		if (t > i) {
			ret = mod - ret;
			lep (j, i, n) swap (b[i][j], b[t][j]);
			lep (j, 1, n) swap (a[i][j], a[t][j]);
		}
		ret = mul (ret, b[i][i]);

		int inv = qpow (b[i][i], mod - 2);
		lep (j, i, n) b[i][j] = mul (b[i][j], inv);
		lep (j, 1, n) a[i][j] = mul (a[i][j], inv);

		lep (j, i + 1, n) if (b[j][i]) {
			int buf = b[j][i];
			lep (k, i, n) sub (b[j][k], mul (buf, b[i][k]));
			lep (k, 1, n) sub (a[j][k], mul (buf, a[i][k]));
		}
	}
	lep (i, 1, n) lep (j, 1, n) a[i][j] = a[i][j] ? mod - a[i][j] : 0;
	auto q = char_poly (n, a);
	vector <int> p (n + 1);
	lep (i, 0, n - offset) p[i] = mul (q[i + offset], ret);
	return p;
}