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;
}