// created on 23.12.11
Pfaffian And Determinant of \(\sum A_iz^i\)
求斜对称矩阵的 Pfaffian 即 \(\mathrm{Pf}(A)\)，主要用到了性质 \(\mathrm{Pf}(A)|B|=\mathrm{Pf}(BAB^{T})\) 。
我们对矩阵进行消元。若第一列全 \(0\)，则矩阵不满秩，\(\mathrm{Pf}(A)=0\) 。否则找到任一第一列元素非 \(0\) 的行，交换到第二行。同时执行 \(B^T\)，将列交换。
我们用第二行第一列的非 \(0\) 元素将后面行的第一列全部消掉。在执行 \(B^T\) 时，由于斜对称性，也可以恰好消掉第一列所有元素。
接着用第一行第二列的非 \(0\) 元素将后面行的第二列全部消掉。在执行 \(B^T\) 时，由于斜对称性，也可以恰好消掉第二行所有元素。
此时通过 Pfaffian 的展开式 \(\mathrm{Pf}(A)=\sum\limits_{i\not=1}(-1)^{i}A_{1,i}\mathrm{Pf}(A\backslash\{1,i\})\) 有 \(\mathrm{Pf}(A)=A_{1,2}\mathrm{Pf}(A\backslash\{1,2\})\) 。于是我们得到了 Pfaffian 的 \(O(n^3)\) 解法。
// Pf (a), O (n ^ 3)
template <class T> inline int pfaffian (T a, int n) {
int prod = 1;
for (int h = 1; h < n; h += 2) {
lep (i, h + 1, n) if (a[i][h]) {
if (i != h + 1) {
prod = mod - prod, swap (a[i], a[h + 1]);
lep (j, h, n) swap ( a[j][i], a[j][h + 1] );
}
break ;
}
if (! a[h][h + 1]) return 0;
prod = mul (prod, a[h][h + 1]);
int buf = qpow (a[h][h + 1], mod - 2);
lep (j, h, n) {
a[h][j] = mul (a[h][j], buf);
a[j][h] = mul (a[j][h], buf);
}
lep (j, h + 2, n) if (buf = a[j][h]) {
lep (k, h, n) pls (a[j][k], mul (a[h + 1][k], buf));
lep (k, h, n) pls (a[k][j], mul (a[k][h + 1], buf));
}
lep (j, h + 2, n) if (buf = a[j][h + 1]) {
lep (k, h, n) sub (a[j][k], mul (a[h][k], buf));
lep (k, h, n) sub (a[k][j], mul (a[k][h], buf));
}
}
return prod;
}
接着考虑 \(\mathrm{Pf}(\sum A_iz^i)\) 的问题。在 网格图最大流计数 中，\([z^0]\mathrm{Pf}(\sum A_iz^i)=1\) 时必然存在的。于是问题实际转化为求解 \(\det(\sum A_iz_i)\) 。
实际上，只需要通过消元找到 \(MA_m=I\)，就可以使得 \(\det(\sum\limits_{i=0}^{m}A_iz^i)\) 转化为 \(\det(M^{-1})\det(M\sum\limits_{i=0}^{m-1}A_iz^i+Iz^m)\) 。在这个过程中，如果找不到某个 \(A_{m,i,i}\not=0\)，则将这一列乘 \(z\) 继续考虑。如果最后成为全 \(0\) 列意味着不满秩，原式为 \(0\) 。
最后问题就转化为 \(\det(\sum\limits_{i=0}^{m}A_iz^i+Iz^m)\) 形式的问题。我们设计分块矩阵 \(\begin{bmatrix}Iz& -I& 0&\cdots& 0&0\\0& Iz& -I&\cdots& 0&0\\\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\0&0&0&\cdots&Iz&-I\\A_0&A_1&A_2&\cdots& A_{m-2}&A_{m-1}+Iz\end{bmatrix}\)，其中 \(I\) 和 \(A_i\) 均为 \(n\) 阶方阵。由于 \(Iz\) 有逆，简单归纳可以发现以上矩阵的行列式就是原问题的答案。而做法是 trivial 的：先消元得到上海森堡矩阵，然后递推即可。这个过程和求解特征多项式是一模一样的。
于是我们得到了 \(\det(\sum\limits_{i=0}^{m}A_iz^i)\) 的 \(O((mn)^3)\) 解法。
// det (xI + A), O (n ^ 3)
template <class T> inline vector <int> char_poly (T a, int n) {
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));
}
}
vector <vector <int> > p (n + 2, vector <int> (n + 1, 0) );
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; ; 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));
if (j == n) break ;
}
}
return p[1];
}
// det (a0 + a1 x + a2 x ^ 2 + ... ak x ^ k), O ((nk) ^ 3)
template <class T> inline vector <int> det (T a, int k, int n) {
int prod = 1, offset = 0;
lep (h, 1, n) {
lep (i, 1, h - 1) {
int buf = a[k][i][h];
lep (d, 0, k) lep (j, 1, n) {
sub (a[d][j][h], mul (a[d][j][i], buf));
}
}
for ( ; ; ) {
if (a[k][h][h]) break ;
lep (i, h + 1, n) if (a[k][i][h]) {
prod = mod - prod;
lep (d, 0, k) lep (j, 1, n) swap (a[d][h][j], a[d][i][j]);
}
if (a[k][h][h]) break ;
if ( ++ offset > k * n ) return vector <int> (k * n + 1, 0);
rep (d, k, 1) lep (j, 1, n) a[d][j][h] = a[d - 1][j][h];
lep (j, 1, n) a[0][j][h] = 0;
lep (i, 1, h - 1) {
int buf = a[k][i][h];
lep (d, 0, k) lep (j, 1, n) {
sub (a[d][j][h], mul (a[d][j][i], buf));
}
}
}
if (! a[k][h][h]) return vector <int> (k * n + 1, 0);
prod = mul (prod, a[k][h][h]);
int inv = qpow (a[k][h][h], mod - 2);
lep (d, 0, k) lep (j, 1, n) a[d][h][j] = mul (a[d][h][j], inv);
lep (i, 1, n) if (i != h) {
int buf = a[k][i][h];
lep (d, 0, k) lep (j, 1, n) {
sub (a[d][i][j], mul (a[d][h][j], buf));
}
}
}
vector <vector <int> > b (k * n + 1, vector <int> (k * n + 1, 0) );
lep (i, 1, (k - 1) * n) b[i][n + i] = mod - 1;
lep (d, 0, k - 1) {
lep (i, 1, n) lep (j, 1, n) {
b[(k - 1) * n + i][d * n + j] = a[d][i][j];
}
}
auto ret = char_poly (b, k * n);
lep (i, offset, k * n) ret[i - offset] = mul (ret[i], prod);
lep (i, k * n - offset + 1, k * n) ret[i] = 0;
return ret;
}