Determinant of $A+Bz$
// created on 23.09.22
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\) 要么列的选择方式，可以得到：
于是以 \(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;
}
