det(A+Bx) 及其推广
矩阵的特征多项式
对于一个矩阵 \(\mathbf A\)，其特征多项式为 \(f_{\mathbf A}(\lambda) = \det(\lambda \mathbf I - \mathbf A)\)。
注意到 \(\det(\lambda \mathbf I - \mathbf A) = \det(\mathbf P (\lambda \mathbf I - \mathbf A) \mathbf P^{-1}) = \det(\lambda \mathbf I - \mathbf P\mathbf A \mathbf P^{-1})\)，这就是说我们可以通过相似变换将 \(\mathbf A\) 相似成上海森堡矩阵。
对于上海森堡矩阵求行列式，这个是经典问题，可以 \(O(n^3)\) DP 求出答案。
std::vector<Z> charPoly(std::vector<std::vector<Z>> A) {
const size_t n = A.size();
for (size_t i = 0; i < n; i ++) {
for (size_t j = 0; j < n; j ++) {
A[i][j] = -A[i][j];
}
}
for (size_t i = 0; i + 2 < n; i ++) {
size_t pivot = i + 1;
for (; pivot < n && !A[pivot][i]; ++ pivot);
if (pivot == n) {
continue;
}
if (pivot > i + 1) {
for (size_t j = i; j < n; j ++) {
std::swap(A[i + 1][j], A[pivot][j]);
}
for (size_t j = 0; j < n; j ++) {
std::swap(A[j][i + 1], A[j][pivot]);
}
}
const Z inv = A[i + 1][i].inv();
for (size_t j = i + 2; j < n; j ++) {
if (A[j][i]) {
const Z t = A[j][i] * inv;
for (size_t k = i; k < n; k ++) {
A[j][k] -= t * A[i + 1][k];
}
for (size_t k = 0; k < n; k ++) {
A[k][i + 1] += t * A[k][j];
}
}
}
}
std::vector<std::vector<Z>> dp(n + 1);
dp[0] = {Z(1)};
for (size_t i = 0; i < n; i ++) {
dp[i + 1].assign(i + 2, Z(0));
for (size_t k = 0; k <= i; k ++) {
dp[i + 1][k + 1] = dp[i][k];
}
for (size_t k = 0; k <= i; k ++) {
dp[i + 1][k] += A[i][i] * dp[i][k];
}
Z prod = 1;
for (size_t j = i; j; j --) {
prod *= -A[j][j - 1];
const Z t = prod * A[j - 1][i];
for (size_t k = 0; k < j; k ++) {
dp[i + 1][k] += t * dp[j - 1][k];
}
}
}
return dp[n];
}
\(\det(\mathbf A + \mathbf Bx)\)
想法是把它转化成特征多项式问题，也就是说我们需要把 \(\mathbf B\) 消成 \(\mathbf I\)。
直接高消的话会有一点小问题，就是如果有一行的 \(\mathbf B\) 全是 \(0\)，那你就消不动了。解决方法是把这一行乘 \(x\)，然后用 \(B\) 的前 \(i - 1\) 行把这一行的前 \(i - 1\) 列消掉，最后把答案除以 \(x\)。一直进行直到此行不全 \(0\)，或除的 \(x\) 次数 \(>n\)。
高消完成之后剩下的就是特征多项式问题，总时间复杂度 \(O(n^3)\)。
\(\det(\mathbf A_0 + \mathbf A_1 x + \mathbf A_2 x^2 + \cdots + \mathbf A_k x^k)\)
拉格朗日插值
首先我们知道答案一定是 \(nk\) 次多项式，所以我们可以带入 \(nk\) 个点值然后拉插。
时间复杂度 \(O(n^4k + n^3k^2)\)。
转化为特征多项式
同样地，先把 \(\mathbf A_k\) 消成 \(\mathbf I\)，这样问题就转化为了求 \(\det(x^k \mathbf I - \mathbf F_1 x^{k - 1} - \cdots - \mathbf F_k)\)。
如果记：
那么我们就有 \(\det(x^m \mathbf I - \mathbf F_1 x^{m - 1} - \cdots - \mathbf F_k) = \det(x \mathbf I - \mathbf A)\)。
\(\text{Proof}\)：首先我们知道：
\[x\mathbf{I} - \mathbf{A} = \begin{bmatrix}x\mathbf{I} & -\mathbf{I} & & & \\0 & x\mathbf{I} & -\mathbf{I} & & \\\vdots & \vdots & \ddots & \ddots & \\0 & 0 & \cdots & x\mathbf{I} & -\mathbf{I} \\-\mathbf{F}_k & -\mathbf{F}_{k-1} & \cdots & -\mathbf{F}_2 & x\mathbf{I} - \mathbf{F}_1\end{bmatrix} \]可以将矩阵分块为：
\[x\mathbf{I} - \mathbf{A} = \begin{bmatrix}\mathbf{P} & \mathbf{Q} \\\mathbf{R} & \mathbf{S}\end{bmatrix} \]\(\mathbf P\) 是左上角的 \((k - 1)n \times (k - 1)n\) 的分块矩阵：
\[\mathbf{P} = \begin{bmatrix}x\mathbf{I} & -\mathbf{I} & & \\0 & x\mathbf{I} & \ddots & \\\vdots & \ddots & \ddots & -\mathbf{I} \\0 & \cdots & 0 & x\mathbf{I}\end{bmatrix} \]\(\mathbf Q\) 是右上角的 \((k - 1)n \times n\) 的分块矩阵：
\[\mathbf{Q} = \begin{bmatrix}0 \\\vdots \\0 \\-\mathbf{I}\end{bmatrix} \]\(\mathbf R\) 是左下角 \(n \times (k - 1) n\) 的分块矩阵：
\[\mathbf{R} = \begin{bmatrix}-\mathbf{F}_k & -\mathbf{F}_{k-1} & \cdots & -\mathbf{F}_2\end{bmatrix} \]\(\mathbf S\) 是右下角的 \(n \times n\) 的分块矩阵：
\[\mathbf{S} = x\mathbf{I} - \mathbf{F}_1 \]显然 \(\mathbf P\) 是可逆的，那么行列式就可以表示为：
\[\det(x \mathbf I - \mathbf A) = \det(\mathbf P) \cdot \det(\mathbf S - \mathbf R \mathbf P^{-1}\mathbf Q) \]首先可以简单得到：
\[\mathbf{P}^{-1} = \begin{bmatrix}x^{-1}\mathbf{I} & x^{-2}\mathbf{I} & x^{-3}\mathbf{I} & \cdots \\0 & x^{-1}\mathbf{I} & x^{-2}\mathbf{I} & \cdots \\\vdots & \ddots & \ddots & \vdots \\0 & \cdots & 0 & x^{-1}\mathbf{I}\end{bmatrix} \]然后可以计算得出：
\[\begin{aligned}\mathbf{R} \mathbf{P}^{-1} \mathbf{Q} &= \mathbf{F}_k x^{-k} + \mathbf{F}_{k-1} x^{-(k-1)} + \cdots + \mathbf{F}_2 x^{-2}\\ \mathbf{S} - \mathbf{R} \mathbf{P}^{-1} \mathbf{Q} &= x\mathbf{I} - \mathbf{F}_1 - \mathbf{F}_k x^{-k} - \mathbf{F}_{k-1} x^{-(k-1)} - \cdots - \mathbf{F}_2 x^{-2}\end{aligned} \]由于 \(\mathbf P\) 是分块上三角矩阵，那么就有 \(\det(\mathbf P) = x^{n(k - 1)}\)，这样就能得到：
\[\begin{aligned}\det(x \mathbf I - \mathbf A) &= x^{n(k - 1)} \cdot \det\left(x\mathbf{I} - \sum_{1 \le i \le k} \mathbf{F}_i x^{-(i-1)}\right) \\&= \det\left(x^k \mathbf{I} - \mathbf{F}_1 x^{k-1} - \mathbf{F}_2 x^{k-2} - \cdots - \mathbf{F}_k\right) \\\end{aligned} \]
这样时间复杂度是 \(O((nk)^3)\) 的。
std::vector<Z> detPoly(std::vector<std::vector<std::vector<Z>>> A) {
assert(A.size());
const int m = A.size() - 1, n = A[0].size();
Z prod = Z(1); int offset = 0;
for (int h = 0; h < n; h ++) {
for (; ; ) {
if (A[m][h][h]) {
break;
}
for (int j = h + 1; j < n; j ++) {
if (A[m][h][j]) {
prod = -prod;
for (int d = 0; d <= m; d ++) {
for (int i = 0; i < n; i ++) {
std::swap(A[d][i][h], A[d][i][j]);
break;
}
}
}
}
if (A[m][h][h]) {
break;
}
if (++ offset > n * m) {
return std::vector<Z>(n * m + 1, Z(0));
}
for (int d = m; -- d >= 0; ) {
for (int j = 0; j < n; j ++) {
A[d + 1][h][j] = A[d][h][j];
}
}
for (int j = 0; j < n; j ++) {
A[0][h][j] = Z(0);
}
for (int i = 0; i < h; i ++) {
const Z t = A[m][h][i];
for (int d = 0; d <= m; d ++) {
for (int j = 0; j < n; j ++) {
A[d][h][j] -= t * A[d][i][j];
}
}
}
}
prod *= A[m][h][h];
const Z inv = A[m][h][h].inv();
for (int d = 0; d <= m; d ++) {
for (int j = 0; j < n; j ++) {
A[d][h][j] *= inv;
}
}
for (int i = 0; i < n; i ++) {
if (h != i) {
const Z t = A[m][i][h];
for (int d = 0; d <= m; d ++) {
for (int j = 0; j < n; j ++) {
A[d][i][j] -= t * A[d][h][j];
}
}
}
}
}
std::vector<std::vector<Z>> B(n * m, std::vector<Z>(n * m, Z(0)));
for (int i = 0; i < (m - 1) * n; i ++) {
B[i][n + i] = Z(1);
}
for (int d = 0; d < m; d ++) {
for (int i = 0; i < n; i ++) {
for (int j = 0; j < n; j ++) {
B[(m - 1) * n + i][d * n + j] = -A[d][i][j];
}
}
}
const std::vector<Z> f = charPoly(B);
std::vector<Z> g(n * m + 1, Z(0));
for (int i = 0; i + offset <= n * m; i ++) {
g[i] = prod * f[i + offset];
}
return g;
}
浙公网安备 33010602011771号