题解：P6178 【模板】Matrix-Tree 定理

· · 题解

Blog

1. 基础概念

根据图论里的知识，我们已经知道了邻接矩阵 A 的定义。如果是有向图，那么边 u_i\to v_i 只需要在 A_{u_i, v_i} 处加上 w_i 即可；如果是无向图，那么边 u_i\leftrightarrow v_i 只需要在 A_{u_i, v_i},A_{v_i, u_i} 处全都加上 w_i 即可。

这里我们引入一个新的概念：度数矩阵 D。

对于无向图而言，度数矩阵的定义为 D_{i, j} = \left\{\begin{matrix}\mathrm{deg}_i & (i=j)\\0 & (i\ne j)\end{matrix}\right.。

而对于有向图而言，度数矩阵则分为入度矩阵和出度矩阵。具体而言，入度矩阵的定义为 D^{\text{in}}_{\,i, j} = \left\{\begin{matrix}\mathrm{deg^{\text{in}}}_i & (i=j)\\0 & (i\ne j)\end{matrix}\right.，出度矩阵的定义为 D^{\text{out}}_{\,i, j} = \left\{\begin{matrix}\mathrm{deg^{\text{out}}}_i & (i=j)\\0 & (i\ne j)\end{matrix}\right.。

而对于带权图，此处的度数定义为对应边的权值之和。

由此可以引出图的 Laplace 矩阵（拉普拉斯矩阵，又称 Kirchhoff 基尔霍夫矩阵）：

最后还有一个矩阵：图的关联矩阵 M。这是一个 e\times n 的矩阵，其中 e 表示边数，n 表示点数。它依然分为入度关联矩阵、出度关联矩阵，具体定义如下：

M^{\mathrm{in}}_{\,i, j} = \left\{\begin{matrix}\sqrt{w_i} & (j = v_i)\\0 & (j \ne v_i)\end{matrix}\right. \\ M^{\mathrm{out}}_{\,i, j} = \left\{\begin{matrix}\sqrt{w_i} & (j = u_i)\\0 & (j \ne u_i)\end{matrix}\right.

其中 u_i, v_i 分别表示第 i 条边 e_i 的起点、终点。定义中会出现根号的原因是，我们后面需要用到平方的性质把边权还原为 w_i。

2. Matrix-Tree 定理

定理 1（Matrix-Tree 的无向图形式）：对于一张给定的带权无向图 G，其生成树个数 t(G) 为 \det L(G)_{[n]\setminus \{k\}, [n]\setminus \{k\}}。其中 [n] = [1, n]\cap\Z, k\in [n]。也就是说，图 G 的生成树个数等于其 Laplace 矩阵的任意一个 n-1 阶主子式的值。

定理 2（Matrix-Tree 的根向树形图形式）：对于一张给定的带权有向图 G，其以 k 为根的根向树形图个数 t^{\mathrm{root}}(G, k) 为 \det L^{\mathrm{out}}(G)_{[n]\setminus \{k\}, [n]\setminus \{k\}}。其中 [n] = [1, n]\cap\Z, k\in [n]。也就是说，图 G 以 k 为根的根向树形图个数等于其出度 Laplace 矩阵删去第 k 行、第 k 列的 n-1 阶主子式的值。

此处“根向树形图”又称内向树。同理下文中“叶向树形图”又称外向树。

定理 3（Matrix-Tree 的叶向树形图形式）：对于一张给定的带权有向图 G，其以 k 为根的叶向树形图个数 t^{\mathrm{leaf}}(G) 为 \det L^{\mathrm{in}}(G)_{[n]\setminus \{k\}, [n]\setminus \{k\}}。其中 [n] = [1, n]\cap\Z, k\in [n]。也就是说，图 G 以 k 为根的叶向树形图个数等于其入度 Laplace 矩阵删去第 k 行、第 k 列的 n-1 阶主子式的值。

3. Matrix-Tree 定理证明

容易发现定理 1、定理 3 全部都可以归纳至定理 2 的情况，因此我们只需要对定理 2 进行证明即可。

下文中，引理 1\sim 4 是用来证明 Cauchy-Binet 公式的，引理 5\sim 6 与证明 Matrix-Tree 的证明有直接关系。

引理 1：对于 n\times m 的矩阵 A 和 m\times n 的矩阵 B，有 |I_n - AB| = |I_m - BA|。

证明：
对分块矩阵的行列式 \begin{vmatrix} I_m & B \\ A & I_n\end{vmatrix} 进行变换。具体而言，分别将第一行倍加到第二行、将第二行倍加到第一行：

\begin{vmatrix} I_m & B \\ A & I_n\end{vmatrix} = \begin{vmatrix} I_m & B \\ O & I_n - AB\end{vmatrix} = \begin{vmatrix} I_m-BA & O \\ A & I_n\end{vmatrix}

将该主对角型行列式算出后得 |I_n - AB| = |I_m - BA|。

引理 2：对于 n\times m 的矩阵 A 和 m\times n 的矩阵 B，有 \lambda^m|\lambda I_n - AB| = \lambda^n|\lambda I_m - BA|。

证明：
对于 \lambda = 0 的情况，显然是成立的。而对于 \lambda\ne 0 的情况，\lambda^{-1} 显然存在，因此有：

\begin{aligned}\lambda^m|\lambda I_n - AB| &= \lambda^{n+m}|I_n - \lambda^{-1}AB| \\ &= \lambda^{n+m}|I_m - \lambda^{-1}BA| \\ &= \lambda^{n}|\lambda I_m - BA|\end{aligned}

其中第二个等号用到了上面引理 1 换元后的结论。

引理 3：对于式子 |\lambda I_m-P| = \sum_{i = 0}^m(-1)^ic_i\lambda^{m-i} ，满足 c_i 为 P 的所有 i 阶主子式之和（当 i = 0 的时候 c_i=1）。

证明：
注意到 \lambda I_m-P 是形如下面的矩阵：

\begin{pmatrix} \lambda-p_{1,1} & -p_{1, 2} & \cdots & -p_{1, m}\\ -p_{2, 1} & \lambda-p_{2,2} & \cdots & -p_{2, m}\\ \vdots & \vdots & \ddots & \vdots \\ -p_{m, 1} & -p_{m, 2} & \cdots & \lambda-p_{m, m}\end{pmatrix}

考虑暴力展开行列式时行顺排的过程，钦定乘积中 \lambda 的幂次为 m-i，那么就相当于钦定了 m-i 个元素必须要在主对角线上，且乘积中选择 \lambda 一项；其余元素无论在哪一列，乘积中一律选择 -p 一项。那么可以得到式子：

(-1)^ic_i=\sum_{\sigma \in S_{m}}\sum_{J_i}\prod_{j = 1}^m [j\notin J_i](-p_{j,\sigma(j)})

其中 J_i 表示所有满足 |J| = m-i 且所有元素均为排列 \sigma 中的不动点的集合 J。
交换求和顺序并化简后可得式子：

(-1)^ic_i=(-1)^{i}\sum_{J_i}\sum_{\sigma \in S_{m}}\prod_{j = 1}^m [j\notin J_i]p_{j,\sigma(j)}

注意到 \sum_{\sigma \in S_{m}}\prod_{j = 1}^m [j\notin J_i]p_{j,\sigma(j)} 部分很像行列式的暴力展开形式，于是可以继续化简为：

(-1)^ic_i = (-1)^{i}\sum_{t \subseteq [m] \wedge |t|=m-i}\det P_{[m]\setminus t,[m]\setminus t}

得到结论：

c_i = \sum_{t\subseteq [m] \wedge |t|=m-i}\det P_{[m]\setminus t,[m]\setminus t}

即：

c_i = \sum_{S\subseteq [m] \wedge |S|=i}\det P_{S,S}

这个结论同样可以引入 Kronecker 符号更加严谨地进行证明，此处仅为简单解释一下原理。

引理 4（Cauchy-Binet 公式）：对于 n\times m 的矩阵 A 和 m\times n 的矩阵 B，有：

\det AB = \sum_{|S|\subseteq [m]\wedge |S| = n}\det A_{[n], S}\det B_{S, [n]}

证明：
根据前文的引理 2，有 \lambda^m|\lambda I_n - AB| = \lambda^n|\lambda I_m - BA|，移项后得：

\lambda^{m-n}|\lambda I_n - AB| = |\lambda I_m - BA|

比较左右两侧 \lambda^{m-n} 的系数。容易发现等式左侧 \lambda^{m-n} 的系数相当于要求 |\lambda I_n - AB| 中 \lambda 的幂次为 0，因此等式左侧的系数为 (-1)^n\det AB；而等式右侧 \lambda^{m-n} 的系数可以套用引理 3 的结论，所以系数为 (-1)^nc_n = (-1)^n\sum_{S\subseteq [m] \wedge |S|=i}\det P_{S,S}。
显然等式左右两侧 \lambda 对应幂次的系数相等，因此有：

(-1)^n\det AB = (-1)^n\sum_{S\subseteq [m] \wedge |S|=n}\det (BA)_{S,S}

化简得：

\det AB = \sum_{S\subseteq [m] \wedge |S|=n}\det (BA)_{S,S}

因为留下的是主子式，所以有：

\det (BA)_{S,S} = \det{B_{S, [n]}A_{[n],S}} = \det{B_{S, [n]}} \times \det{A_{[n], S}}=\det{A_{[n], S}} \times \det{B_{S, [n]}}

因此得到结论：

\det AB = \sum_{S\subseteq [m] \wedge |S|=n}\det{A_{[n], S}} \times \det{B_{S, [n]}}

引理 5：对于关联矩阵 M，有如下三个式子成立：

\begin{cases}D^{\mathrm{out}} = (M^{\mathrm{out}})^{\mathrm{T}} M^{\mathrm{out}} \\D^{\mathrm{in}} = (M^{\mathrm{in}})^{\mathrm{T}} M^{\mathrm{in}} \\A = (M^{\mathrm{out}})^{\mathrm{T}} M^{\mathrm{in}}\end{cases}

证明很容易，直接考虑矩阵乘法的过程，模拟就能发现两个 \sqrt{w_i} 相乘还原为了 w_i，而不存在边的点对一定是 0。

根据这个结论我们还可以得出：

\begin{cases}L^{\mathrm{out}}=D^{\mathrm{out}}-A=(M^{\mathrm{out}})^{\mathrm{T}} M^{\mathrm{out}}-(M^{\mathrm{out}})^{\mathrm{T}} M^{\mathrm{in}}=(M^{\mathrm{out}})^{\mathrm{T}}(M^{\mathrm{out}}-M^{\mathrm{in}}) \\ L^{\mathrm{in}}=D^{\mathrm{in}}-A=(M^{\mathrm{in}})^{\mathrm{T}} M^{\mathrm{in}}-(M^{\mathrm{out}})^{\mathrm{T}} M^{\mathrm{in}}=(M^{\mathrm{in}}-M^{\mathrm{out}})^{\mathrm{T}}M^{\mathrm{in}} \\ \end{cases}

引理 6：对于有向图 G 的一个子图 (W, S)，若满足 |W| = |S| \le n，则子图 T = (W, S) 是一个以 V\setminus W 之中任意一点作为根的根向树形图，当且仅当 \det (M^{\mathrm{out}})^{\mathrm{T}}\det (M^{\mathrm{out}}-M^{\mathrm{in}}) \ne 0。且当 \det (M^{\mathrm{out}})^{\mathrm{T}}\det (M^{\mathrm{out}}-M^{\mathrm{in}}) \ne 0 的时候，必然有 \det (M^{\mathrm{out}})^{\mathrm{T}}\det (M^{\mathrm{out}}-M^{\mathrm{in}}) = \prod_{e\in S}w(e)，即行列式的值等于树的权值 w(T)。

证明：
首先我们可以将子图中所有的边权全部转化为 1。原因在于这两个行列式中的每一行都可以提出一个公因数 \sqrt{w_i}，且这些公因数的乘积等于树的权值 w(T)。
然后先对 \det (M^{\mathrm{out}})^{\mathrm{T}} 进行分析。显然里面没有任何一列全部为 0。但是如果存在任意一行全部为 0 的时候（即一个点没有出边），该行列式的值一定为 0，此时一定无法组成根向树形图。但是我们无法得到其他信息了，只能得到结论：每个点的出边数为 1 是形成根向树形图的必要不充分条件。
又因为当 \det (M^{\mathrm{out}})^{\mathrm{T}}\ne 0 的时候，每个点的出边数必然为 1。所以我们只需要判断图中是否存在环即可。
再对 \det (M^{\mathrm{out}}-M^{\mathrm{in}}) 进行分析。对每一行（每条边）进行考虑，只要这条边直接指向的不是根（根不在集合 W 中），那么必然存在一个数是 1，另一个数为 -1。我们考虑消掉这个 -1，具体而言，假设当前的边是 a\to b，因为 b\ne \text{root}，所以必然存在另一条边 b\to c。此时第二条边所在行的第 b 列必然为 1，而第一条边所在行的第 b 列为 -1，利用行列式的倍加性质，可以使得第一条边所在行的第 b 列变成 0，而第 c 列变成 -1。如此往复，直到“第二条边”直接指向根节点，第 c 列就会变成 0 了。此时这一行只有第 a 列为 1。这个构造的过程很像不断往父亲处跳，直到跳到根节点的过程。
按照这种方法模拟下去，如果某一行被完全消成了 0，就说明这张图存在环，行列式一定为 0。否则每一行只有起点所在列是 1，其余全部为 0，行列式显然为 1。

接下来就可以来证明 Matrix-Tree 定理了，我们以定理 2 为例，来推式子：

\begin{aligned}\det L^{\mathrm{out}}(G)_{[n]\setminus \{k\}, [n]\setminus \{k\}} &= \det L^{\mathrm{out}}(G)_{W, W} \\ &= \det{(M^{\mathrm{out}})^{\mathrm{T}}(M^{\mathrm{out}}-M^{\mathrm{in}})} \\ & = \sum_{S\subseteq [m] \wedge |S|=n-1} \det{(M^{\mathrm{out}})^{\mathrm{T}}}_{W, S}\det{(M^{\mathrm{out}}-M^{\mathrm{in}})}_{S, W} \end{aligned}

其中 W = [n] \setminus \{k\}。第二个等号用到了引理 5 的推论，第三个等号用到了引理 4 的 Cauchy-Binet 公式。

根据引理 6，求和中的式子 \det{(M^{\mathrm{out}})^{\mathrm{T}}}_{W, S}\det{(M^{\mathrm{out}}-M^{\mathrm{in}})}_{S, W}=w(T) 当且仅当子图 (W, S) 构成了一个以 k 为根的根向树形图，其余情况下该式子均为 0。所以相当于对所有以 k 为根的根向树形图权值求和了。就此定理 2 得证。

接下来的定理 1 其实可以由定理 2 证明。因为无向图的生成树并没有钦定根，所以任选一个 n-1 阶主子式相当于在有向图中随便钦定一个根，得到的数目就是无向图中生成树的个数。而定理 3 相当于就是把边反向了一下，用同样的方式容易证明。

4. Cayley 公式

定理 4（Cayley 公式）：包含 n 个节点的有标号无向无根树一共有 n^{n-2} 个。

证明：
原命题相当于给定了一张 n 个点的无向完全图，求生成树个数。
考虑套用矩阵树定理，首先写出矩阵的形态：

\begin{pmatrix} n-1 & -1 & \cdots & -1 \\ -1 & n-1 & \cdots & -1 \\ \vdots & \vdots & \ddots & \vdots\\ -1 & -1 & \cdots & n-1\end{pmatrix}_{n\times n}

把第一行第一列去掉，得到：

\begin{pmatrix} n-1 & -1 & \cdots & -1 \\ -1 & n-1 & \cdots & -1 \\ \vdots & \vdots & \ddots & \vdots\\ -1 & -1 & \cdots & n-1\end{pmatrix}_{(n-1)\times (n-1)}

现在我们只需要求出它的行列式即可。
这显然是行列式求值里“列和相等”的模型。具体而言，我们注意到每一列元素之和相等，于是把所有行全部利用倍加性质加到第一行上，得：

\begin{vmatrix} 1 & 1 & \cdots & 1 \\ -1 & n-1 & \cdots & -1 \\ \vdots & \vdots & \ddots & \vdots\\ -1 & -1 & \cdots & n-1\end{vmatrix}_{(n-1)\times (n-1)}

接下来再用第一行倍加到其他行上，把 -1 项消掉，得：

\begin{vmatrix} 1 & 1 & \cdots & 1 \\ 0 & n & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & n\end{vmatrix}_{(n-1)\times (n-1)}

这是一个主对角型行列式，只要将对角线上的元素乘起来就能得到行列式的值，即 n^{n-2}。

5. BEST 定理

定理 5（有向欧拉回路的判定法则）：若有向图 G 强连通，且对于 G 中的每个点均有 \mathrm{deg}^{\mathrm{in}} = \mathrm{deg}^{\mathrm{out}}。

根据这个定理，我们可以延伸出 BEST 定理的作用：求有向图中欧拉回路的个数。

需要注意的是，普通无向图中欧拉回路的计数是 NPC，无法在多项式时间复杂度内解决。

定理 6（BEST 定理）：若有向图中存在欧拉回路，那么它的欧拉回路个数为 \mathrm{ec}(G) = t^{\mathrm{root}}(G, k)\times \prod_{i=1}^n (\deg(i) - 1)!。其中 k 为 1\sim n 中的任意一个正整数。

其中，\deg(v) 表示点 v 的出度 / 入度。因为当有向图存在欧拉回路的时候，必然有 \deg^{\mathrm{in}} = \deg^{\mathrm{out}}。同理，t^{\mathrm{root}}(G, k) 和 t^{\mathrm{leaf}}(G, k) 在标准的 BEST 定理上也是等价的。

证明：
将定理的形式转化为 \deg(k)\times \mathrm{ec}(G) = t^{\mathrm{root}}(G, k)\times \deg(k)!\times \prod_{1\le i \le n, i\ne k} (\deg(i) - 1)!。
等式右侧的意义为：从欧拉回路中选出每个点（点 k 除外）最后走的那条边，这些边一定构成了一个以 k 为根的根向树形图。然后点 k 和其他点都随便钦定出边的顺序，按照每个点出边的顺序走就一定能构成欧拉回路。
等式左侧还有个 \deg(k) 的原因是，循环同构的欧拉回路算作同一个欧拉回路，需要用除法去掉。例如 1\to 2\to 1\to 3\to 1 和 1\to 3\to 1\to 2\to 1 是同一个欧拉回路。
接下来证明除 k 以外，其余点最后走的那条边一定构成根向树形图。
因为一共只有 n-1 条边，所以我们只需要证明子图中不存在环即可。
考虑反证法，如果子图中存在环，那么想要构成欧拉回路，环上必然存在一条边，排在所谓的“最后一条边”后面，并且直接或者间接地指向 k。那么就说明了这些边并不是最后走的一条边，与题设矛盾，因此子图中不存在环。
注意这只是一个不严谨的证明，想要看严谨证明的可以参考 OIwiki 中的双射构造。

需要注意的是，如果我们不把循环同构的回路当做同一个欧拉回路，那么等式左侧就不要加那个 \deg(k)，那么公式就变为了 \mathrm{ec}(G) = t^{\mathrm{root}}(G, k)\times \deg(k)!\times \prod_{1\le i \le n, i\ne k} (\deg(i) - 1)!。

6. 例题

打星号的不是我自己想出来的题。

6.1 P6178 【模板】Matrix-Tree 定理

最模板的矩阵树定理，直接套用定理 1 和定理 3 中的结论即可。因为需要使用行列式求值的板子，所以时间复杂度为 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 305;
const ll mod = 1e9 + 7;
ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return res;
}
ll Det(int n, ll a[N][N])
{
    bool flag = 0; ll ans = 1;
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n && a[i][i] == 0; j++)
            if(a[j][i])
                swap(a[i], a[j]), flag ^= 1;
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * qpow(a[i][i], mod - 2) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = ((a[j][k] - c * a[i][k]) % mod + mod) % mod;
        }
        ans = (ans * a[i][i]) % mod;
    }
    return (flag ? (mod - ans) % mod : ans);
}
int n, m, t;
ll a[N][N], b[N][N];
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> m >> t;
    for(int i = 1; i <= m; i++)
    {
        ll u, v, w;
        cin >> u >> v >> w;
        if(!t)
        {
            a[u][u] += w; 
            a[v][v] += w; 
            a[u][v] -= w;
            a[v][u] -= w;
        }
        else
        {
            a[v][v] += w; 
            a[u][v] -= w;
        }
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = (a[i + 1][j + 1] % mod + mod) % mod;
    cout << Det(n - 1, b);
    return 0;
}

6.2 SP104 HIGH - Highways

同样是矩阵树定理无向图形式的模板题，但是因为这题没有模数，所以需要使用辗转相除法求行列式的值。时间复杂度依然是 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 15;
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
        {
            while(a[j][i])
            {
                ll c = a[i][i] / a[j][i];
                for(int k = i; k <= n; k++)
                    a[i][k] -= c * a[j][k];
                Swap(i, j);
            }
        }
        res = res * a[i][i];
    }
    return (flag ? -res : res);
}
ll n, m, a[N][N], b[N][N];
void solve()
{
    cin >> n >> m;
    memset(a, 0, sizeof(a));
    while(m--)
    {
        ll u, v;
        cin >> u >> v;
        a[u][u]++; a[v][v]++;
        a[u][v]--; a[v][u]--;
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = a[i + 1][j + 1];
    cout << Det(n - 1, b) << "\n";
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    int t;
    cin >> t;
    while(t--) solve();
    return 0;
}

6.3 P4111 [HEOI2015] 小 Z 的房间

题意转化为给定一张无向网格图（部分点被删掉），求生成树个数。

注意到数据范围很小，在满数据点下点数也只有 9\times 9=81，因此可以直接套用矩阵树定理求解。

时间复杂度 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 105;
const ll mod = 1e9;
const int gox[] = {0, 0, 1, -1};
const int goy[] = {1, -1, 0, 0};
int n, m, a[N][N], id[N][N], cnt;
ll b[N][N];
bool legal(int x, int y)
{
    return (1 <= x && x <= n && 1 <= y && y <= m);
}
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
        {
            while(a[j][i])
            {
                ll c = a[i][i] / a[j][i];
                for(int k = i; k <= n; k++)
                    a[i][k] = (a[i][k] - c * a[j][k]) % mod;
                Swap(i, j);
            }
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> m;
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j <= m; j++)
        {
            char x;
            cin >> x;
            if(x == '.') id[i][j] = ++cnt;
        }
    }
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j <= m; j++)
        {
            if(!id[i][j]) continue;
            int u = id[i][j];
            for(int k = 0; k < 4; k++)
            {
                int tx = i + gox[k], ty = j + goy[k];
                if(!legal(tx, ty)) continue;
                int v = id[tx][ty];
                if(!v) continue;
                if(u < v) continue;
                a[u][u]++; a[v][v]++;
                a[u][v]--; a[v][u]--;
            }
        }
    }
    for(int i = 1; i < cnt; i++)
        for(int j = 1; j < cnt; j++)
            b[i][j] = a[i + 1][j + 1];
    cout << Det(cnt - 1, b);
    return 0;
}

6.4 P3317 [SDOI2014] 重建

容易将“概率”转化为树边的权值，这样一棵树出现的概率即为权值之积，可以使用矩阵树定理求解。

但是本题并不是模板的矩阵树定理。对于非树边，依然要乘上它不出现在图中的概率。这里我们考虑转化权值，具体而言，令常数 C = \prod_{1\le i < j \le n}(1-p_{i, j})，将树边的权值转化为 \dfrac{p_{i, j}}{1-p_{i, j}}。最后求出所有数权值之积的和之后，将答案整体乘上 C 即可。用乘法分配律容易证明这是正确的。

一个细节是，概率 p_{i, j} 可能为 1，这样分母就会变成 0 了。有如下两种解决方法：

时间复杂度 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 55;
const ldb eps = 1e-10;
int n;
ldb g[N][N], a[N][N], b[N][N], ans = 1;
ldb Det(int n, ldb a[N][N])
{
    bool flag = 0; ldb res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(abs(a[j][i]) > eps)
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        for(int j = i + 1; j <= n; j++)
        {
            ldb c = a[j][i] / a[i][i];
            for(int k = i; k <= n; k++)
                a[j][k] -= c * a[i][k];
        }
        res *= a[i][i];
    }
    return (flag ? -res : res);
}
int main()
{
    scanf("%d", &n);
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j <= n; j++)
        {
            scanf("%Lf", &g[i][j]);
            if(i > j) continue;
            a[i][i] += g[i][j] / (1 - g[i][j] + eps);
            a[j][j] += g[i][j] / (1 - g[i][j] + eps);
            a[i][j] -= g[i][j] / (1 - g[i][j] + eps);
            a[j][i] -= g[i][j] / (1 - g[i][j] + eps);
            ans *= (1 - g[i][j] + eps);
        }
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = a[i + 1][j + 1];
    printf("%.10Lf", ans * Det(n - 1, b));
    return 0;
}

6.5 P4208 [JSOI2008] 最小生成树计数

MST 经典套路题。保证边权的重数不超过 10 是何意味呢？

考虑 Kruskal 的过程，对于边权相等的边，显然我们可以以任意顺序尝试将其加入最小生成树。因此，我们可以假设当前的边权为 w，且边权 < w 的边已经全部尝试过加入最小生成树：

这个过程可以开一个桶和并查集实现。

最后对图 G 做一个矩阵树定理即可。

时间复杂度 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 105;
const ll mod = 31011;
int n, m;
map<int, vector<pi> > tot;
ll a[N][N], b[N][N];
struct DSU{
    int fa[N];
    void init()
    {
        for(int i = 1; i <= n; i++) fa[i] = i;
    }
    int findf(int x)
    {
        if(fa[x] != x) fa[x] = findf(fa[x]);
        return fa[x];
    }
    void combine(int x, int y)
    {
        int fx = findf(x), fy = findf(y);
        fa[fx] = fy;
    }
} dsu;
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i + 1; j <= n; j++)
        {
            while(a[j][i])
            {
                ll c = a[i][i] / a[j][i];
                for(int k = i; k <= n; k++)
                    a[i][k] = (a[i][k] - a[j][k] * c) % mod;
                Swap(i, j);
            }
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> m;
    while(m--)
    {
        ll u, v, w;
        cin >> u >> v >> w;
        tot[w].push_back({u, v});
    }
    dsu.init();
    for(auto itm : tot)
    {
        for(auto eg : itm.se)
        {
            int u = dsu.findf(eg.fi), v = dsu.findf(eg.se);
            if(u == v) continue;
            a[u][u]++; a[v][v]++;
            a[u][v]--; a[v][u]--;
        }
        for(auto eg : itm.se)
        {
            int u = dsu.findf(eg.fi), v = dsu.findf(eg.se);
            if(u == v) continue;
            dsu.combine(u, v);
        }        
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = a[i + 1][j + 1];
    cout << Det(n - 1, b);
    return 0;
}

6.6 P2144 [FJOI2007] 轮状病毒

高精度能不能滚出 OI。

显然可以套用矩阵树定理，把 Laplace 矩阵写出来（n = 5 时的一个例子）：

\begin{pmatrix} n-1 & -1 & -1 & -1 & -1 & -1 \\ -1 & 3 & -1 & & & -1 \\ -1 & -1 & 3 & -1 & & \\ -1 & & -1 & 3 & -1 & \\ -1 & & & -1 & 3 & -1\\ -1 & -1 & & & -1 & 3\end{pmatrix}

容易发现第一行第一列是最不规则的，所以可以选择去掉第一行第一列：

\begin{pmatrix} 3 & -1 & & & -1 \\ -1 & 3 & -1 & & \\ & -1 & 3 & -1 & \\ & & -1 & 3 & -1\\ -1 & & & -1 & 3\end{pmatrix}

我们只需要求这个矩阵的行列式即可。

接下来就是比较机械的过程了。注意到这很像行列式模型里的“川型行列式”，只是左下角和右上角多了一个 -1。于是我们仿照川型行列式的思路，利用行列式展开定理来计算。

首先对第一行展开，行列式的值为：

3\times \begin{vmatrix} 3 & -1 & & \\ -1 & 3 & -1 & \\ & -1 & 3 & -1\\ & & -1 & 3\end{vmatrix}+\begin{vmatrix} -1 & -1 & & \\ & 3 & -1 & \\ & -1 & 3 & -1\\ -1 & & -1 & 3\end{vmatrix}+(-1)^{i+2}\times \begin{vmatrix} -1 & 3 & -1 & \\ & -1 & 3 & -1 \\ & & -1 & 3 \\ -1 & & & -1\end{vmatrix}

注意到其中的第一项就是川型行列式的形式。于是我们设 B_i 表示 i 阶川型行列式的值，那么第一项就被转化为了 3\times B_{i-1}。

接下来对第二项的第一列展开：

\begin{aligned}\begin{vmatrix} -1 & -1 & & \\ & 3 & -1 & \\ & -1 & 3 & -1\\ -1 & & -1 & 3\end{vmatrix}&=(-1)\times\begin{vmatrix} 3 & -1 & \\ -1 & 3 & -1\\ & -1 & 3\end{vmatrix}+(-1)^{i+1}\times \begin{vmatrix} -1 & & \\ 3 & -1 & \\ -1 & 3 & -1\\\end{vmatrix}\\&=(-1)\times B_{i-2}+(-1)^{i+1}\times (-1)^{i-2}\\&=-B_{i-2}-1\end{aligned}

最后对第三项的第一列展开：

\begin{aligned}(-1)^{i+2}\times\begin{vmatrix} -1 & 3 & -1 & \\ & -1 & 3 & -1 \\ & & -1 & 3 \\ -1 & & & -1\end{vmatrix}& = (-1)^{i+3}\times \begin{vmatrix} -1 & 3 & -1 \\ & -1 & 3 \\ & & -1\end{vmatrix}+(-1)^{2i+3}\times \begin{vmatrix} 3 & -1 & \\ -1 & 3 & -1 \\ & -1 & 3 \\ \end{vmatrix}\\&=(-1)^{2i+1}+(-1)^{2i+3}\times B_{i-2}\\&=-B_{i-2}-1\end{aligned}

于是答案即为：

\mathrm{Ans}_i = 3B_{i-1}-2B_{i-2}-2

其中 B_i 川型行列式可以通过递推求解，也可以通过二阶线性递推求出通项公式计算。此处不再赘述。

本题答案较大，需要写高精度。时间复杂度 O(n\lg^2 V)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 105;
ll n;
struct Bigint{
    string num;
    Bigint () {
        num = "0";
    }
    Bigint (ll x) {
        num = "";
        while(x)
        {
            num += (x % 10) + '0';
            x /= 10;
        }
    }
    int getdig(int id) const {
        if(id >= num.length()) return 0;
        return (num[id] - '0');
    }
    Bigint operator + (const Bigint & t) const{
        Bigint tmp;
        tmp.num = "";
        ll jw = 0;
        for(int i = 0; i < max(num.length(), t.num.length()) || jw > 0; i++)
        {
            jw += getdig(i) + t.getdig(i);
            tmp.num += (jw % 10) + '0';
            jw /= 10;
        }
        return tmp;
    }
    Bigint operator - (const Bigint & t) const{
        Bigint tmp;
        tmp.num = "";
        ll jw = 0;
        for(int i = 0; i < max(num.length(), t.num.length()); i++)
        {
            int upv = jw + getdig(i), dwv = t.getdig(i);
            if(upv < dwv) jw = -1, tmp.num += '0' + upv + 10 - dwv;
            else jw = 0, tmp.num += '0' + upv - dwv;
        }
        return tmp;
    }
    Bigint operator * (const ll x) const{
        Bigint tmp;
        tmp.num = "";
        ll jw = 0;
        for(int i = 0; i < num.length() || jw > 0; i++)
        {
            jw += getdig(i) * x;
            tmp.num += '0' + (jw % 10);
            jw /= 10;
        }
        return tmp;
    }
    void outp()
    {
        string tmp = num;
        reverse(tmp.begin(), tmp.end());
        cout << tmp;
    }
} b[N];
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n;
    b[1] = Bigint(3ll); b[2] = Bigint(8ll);
    for(int i = 3; i <= n; i++) b[i] = (b[i - 1] * 3) - b[i - 2];
    if(n == 1) cout << 1;
    else if(n == 2) cout << 5;
    else
    {
        Bigint res = (b[n - 1] * 3) - (b[n - 2] * 2) - Bigint(2ll);
        res.outp();
    }
    return 0;
}

6.7 P4455 [CQOI2018] 社交网络

矩阵树定理求叶向树形图板子，注意是对入度 Laplace 矩阵求行列式。

时间复杂度 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 305;
const ll mod = 1e4 + 7;
ll n, m, a[N][N], b[N][N];
ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return res;
}
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(a[j][i])
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        if(!a[i][i]) return 0;
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * qpow(a[i][i], mod - 2) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[i][k]) % mod;
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> m;
    while(m--)
    {
        int u, v;
        cin >> u >> v;
        a[u][u]++; a[v][u]--;
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = a[i + 1][j + 1];
    cout << Det(n - 1, b);
    return 0;
}

6.8 * P2109 [NOI2007] 生成树计数

这个题，有点深刻了。

很显然可以先用矩阵树定理来刻画问题，我们只需要求下列矩阵的任意一个 n-1 阶主子式即可（以 n = 8, k = 2 为例）。

\begin{bmatrix} 2 & -1 & -1 & & & & & \\ -1 & 3 & -1 & -1 & & & & \\ -1 & -1 & 4 & -1 & -1 & & & \\ & -1 & -1 & 4 & -1 & -1 & & \\ & & -1 & -1 & 4 & -1 & -1 & \\ & & & -1 & -1 & 4 & -1 & -1\\ & & & & -1 & -1 & 3 & -1\\ & & & & & -1 & -1 & 2\end{bmatrix}_{8\times 8}

长得很像川型行列式对吧，掐头去尾行列式展开一下就好了。但是你很快就会发现 k=5 的时候需要展开非常多项，并且要进行很多次分类讨论，代码难度不止一般的大。

我们考虑将它化为主对角型行列式的形式。具体而言，把矩阵分为三部分：

然后我们把第一部分与第二部分进行交换，得到如下的矩阵：

\begin{bmatrix} -1 & -1 & 4 & -1 & -1 & & & \\ & -1 & -1 & 4 & -1 & -1 & & \\ & & -1 & -1 & 4 & -1 & -1 & \\ 2 & -1 & -1 & & & & & \\ -1 & 3 & -1 & -1 & & & & \\ & & & -1 & -1 & 4 & -1 & -1\\ & & & & -1 & -1 & 3 & -1\\ & & & & & -1 & -1 & 2\end{bmatrix}

接下来是这个做法最神奇的部分，我们可以利用前 n-2k-1 行对第 n-2k\sim n-k-1 行进行消元，并且这个消元过程恰好可以使用“矩阵乘法”这个工具进行刻画！

其原因在于，前 n-2k-1 行里每一行的非零部分全都是相同的、连续的。这也就意味着我们可以不断消去非零的第一列，直到这一行非零元素全部移动至最后 2k+1 列。

具体而言，假设第 i 行的前 2k 个元素组成的行向量为：

[a_1, a_2, a_3, a_4, 0, \cdots, 0]

那么我们令第一行乘上 a_1 后倍加到第 i 行上，那么之后的行向量变为了：

[0, a_2-a_1, a_3+2k\times a_1, a_4-a_1, -a_1, 0, \cdots, 0 ]

此时非零的依然只有 2k 个元素，因此我们可以对这 2k 个非零元素组成的行向量写出转移矩阵：

[a_2-a_1, a_3+2k\times a_1, a_4-a_1, -a_1] = [a_1, a_2, a_3, a_4]\times \begin{bmatrix} -1 & 2k & -1 & -1\\ 1 & & & \\ & 1 & & \\ & & 1 &\end{bmatrix}

按照这样的步骤进行下去，最后得到的矩阵即为：

\begin{bmatrix} -1 & -1 & 4 & -1 & -1 & & & \\ & -1 & -1 & 4 & -1 & -1 & & \\ & & -1 & -1 & 4 & -1 & -1 & \\ & & & -24 & 41 & -7 & -10 & \\ & & & 25 & -39 & 5 & 9 & \\ & & & -1 & -1 & 4 & -1 & -1\\ & & & & -1 & -1 & 3 & -1\\ & & & & & -1 & -1 & 2\end{bmatrix}

因为求的是 n-1 阶主子式，我们可以选择把最后一行、最后一列删去。

注意到此时前 n-2k-1 行全都变成了主对角型行列式，而只有后 2k 行是不规则的，且 k 非常小，只有最后 2k 列有值，于是我们可以对后 2k 行做高斯消元至主对角型的形式。再把前面主对角线上的 (-1) 乘上即可。

直接做的时间复杂度是 O(k^4\log n)。但是我们可以注意到这 k 行的转移矩阵是相同的，且初始矩阵为行向量，所以可以直接预处理出转移矩阵的 n-2k-1 次幂，然后左乘行向量（行向量只有第一行有值，所以做矩阵乘法的复杂度为 O(k^2)），最后就得到了 O(k^3\log n) 的做法。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 15;
const ll mod = 65521;
ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return res;
}
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(a[j][i])
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        if(a[i][i] == 0) return 0;
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * qpow(a[i][i], mod - 2) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[i][k]) % mod;
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
namespace Sol1 {
    ll n, m, a[N][N], b[N][N];
    void solve1(ll nnn, ll mmm)
    {
        n = nnn; m = mmm;
        for(int i = 1; i <= n; i++)
        {
            for(int j = i + 1; j <= min(n, i + m); j++)
            {
                a[i][i]++; a[j][j]++;
                a[i][j]--; a[j][i]--;
            }
        }
        for(int i = 1; i < n; i++)
            for(int j = 1; j < n; j++)
                b[i][j] = a[i + 1][j + 1];
        cout << Det(n - 1, b);
    }
}
struct Matrix{
    ll a[N][N];
    Matrix () { memset(a, 0, sizeof(a)); }
    Matrix operator * (const Matrix & t) const{
        Matrix res;
        for(int i = 0; i < N; i++)
            for(int k = 0; k < N; k++)
                for(int j = 0; j < N; j++)
                    res.a[i][j] = (res.a[i][j] + a[i][k] * t.a[k][j]) % mod;
        return res;
    }
};
Matrix Mqpow(Matrix a, ll b)
{
    Matrix res;
    for(int i = 0; i < N; i++) res.a[i][i] = 1;
    while(b)
    {
        if(b & 1) res = res * a;
        b >>= 1;
        a = a * a;
    }
    return res;
}
namespace Sol2{
    ll n, m;
    ll a[N], b[N][N];
    void getprod(ll (&a)[N], Matrix b)
    {
        ll res[N];
        memset(res, 0, sizeof(res));
        for(int i = 1; i <= 2 * m; i++)
            for(int j = 1; j <= 2 * m; j++)
                res[i] = (res[i] + a[j] * b.a[j - 1][i - 1]) % mod;
        memcpy(a, res, sizeof(a));
    }
    void solve2(ll nnn, ll mmm)
    {
        n = nnn; m = mmm;
        bool flag = (((n - 2 * m - 1) * (m + 1)) & 1); ll res = 1;
        Matrix dp;
        for(int i = 0; i < 2 * m; i++) dp.a[0][i] = -1;
        dp.a[0][m - 1] = 2 * m;
        for(int i = 0; i < 2 * m - 1; i++) dp.a[i + 1][i] = 1;
        dp = Mqpow(dp, n - 2 * m - 1);
        for(int i = 1, j = m; i <= m; i++, j++)
        {
            memset(a, 0, sizeof(a));
            a[i] = j;
            b[2 * m - i + 1][2 * m - i + 1] = j + 1;
            for(int k = 1; k <= m + i; k++)
            {
                if(k == i) continue;
                a[k] = -1;
                b[2 * m - i + 1][2 * m - k + 1] = -1;
            }
            getprod(a, dp);
            memcpy(b[i], a, sizeof(a));
        }
        res = Det(2 * m, b);
        if(flag) res = -res;
        cout << (res % mod + mod) % mod;
    }
};
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    ll m, n;
    cin >> m >> n;
    if(n <= 2 * m + 1) Sol1::solve1(n, m);
    else Sol2::solve2(n, m);
    return 0;
}

6.9 * U659114 矩阵树定理？

本题是我在 u 群询问得到的做法，做法并不是我原创的，但是数据是我造的。

第一步是显然的，使用矩阵树定理，问题转化为如何求一个 n 阶矩阵的所有 n-1 阶主子式。

我们考虑强化一下这个问题：如何求这个矩阵的伴随矩阵 A^{\ast} 呢？

根据伴随矩阵的公式 A^{\ast} = \det A\times A^{-1}，我们只需要对 A 分别求逆和求行列式，就能得到 A^{\ast} 的值。但是这仅限 A 满秩时的情况。如果 A 不满秩，说明 A^{-1} 不存在，我们就必须另辟蹊径了。

另外一种情况是 \mathrm{rank}(A) \le n-2。此时我们删去 A 的一行一列，得到的新矩阵 A' 的秩一定小于 n-1，于是行列式的值一定是 0，得到的伴随矩阵就是全 0 的。具体原因可以从“秩的值等于向量组中独立向量的个数”来理解，删去一行一列后秩的数目并不会变大，而 n 只会减小 1，此时矩阵依然不满秩。

这说明我们只需要考虑 \mathrm{rank}(A) = n - 1 时的情况即可。

首先，线性代数里有一个经典的公式：AA^{\ast}=A^{\ast}A = \det A \times I，因为此时 A 不满秩，所以 \det A = 0，进一步得出 AA^{\ast}=A^{\ast}A = 0。

先从 AA^{\ast} = 0 入手，把伴随矩阵拆为向量组，于是对于 A^{\ast} 的每一列，都有 Ax = 0，说明 A^{\ast} 的每一列都在线性映射 A 的零空间里面。

我们把 A^{\ast} 左乘 A 看作对 A^{\ast} 进行线性变换，根据“秩-零度化定理”：\dim V = \mathrm{rank}(T) + \mathrm{nullity}(T)，显然定义域 V 的维度是 n，而 \mathrm{rank}(T) = n-1，说明零空间 \ker T 的维度为 1，即线性映射 A 的零空间是一条直线，可以被表示为 \lambda x 的形式。

同理，根据另一个式子 A^{\ast}A = 0，转化得到 A^T(A^{\ast})^T = 0，之后按照上文的思路得到：线性映射 A^T 的零空间是一条直线，可以被表示为 \lambda x 的形式。

因此伴随矩阵的所有列向量都在 Ax = 0 的解空间里，所有行向量都在 yA=0 的解空间里。因为解空间可以表示为 \lambda x, \lambda y 的形式，所以伴随矩阵就是形如 c\times x \times y 的形式。其中 c 为任意一个常数。

最后的问题在于，我们如何求出这个常数 c 呢？我们只需要求出伴随矩阵中任意一个非零的值 A^{\ast}_{\,i, j}，那么 c = (x_j\times y_i)^{-1}\times A^{\ast}_{\,i, j}。求出非零的位置可以在高斯消元构造行最简形矩阵的过程中一并实现。

时间复杂度 O(n^3)，常数很小，所以可以跑过 n=10^3。std 在未卡常的情况下最大点跑了约 1.5s。

注意这个做法是依赖模数是质数的，因为最后一步需要求逆元。如果模数是合数，我暂时只能想到 O(n^3\log V) 的做法，就是对模数质因数分解，然后对每个质因数做一遍 O(n^3) 的算法，最后用 exCRT 合并答案。不知道是否存在复杂度更优秀的做法。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 1005;
const ll mod = 998244353;
ll ans[N], gx, gy, bel[N], A[N], B[N];
ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return res;
}
ll inv(ll a)
{
    return qpow(a, mod - 2);
}
void GaussX(ll n, const ll a_[N][N])
{
    ll a[N][N];
    memcpy(a, a_, sizeof(a));
    ll rk = 0;
    for(int i = 1; i <= n; i++)
    {
        for(int j = rk + 1; j <= n; j++)
        {
            if(a[j][i])
            {
                if((rk + 1) ^ j) swap(a[rk + 1], a[j]);
                break;
            }
        }
        if(a[rk + 1][i] == 0)
        {
            gy = i;
            continue;
        }
        for(int j = rk + 2; j <= n; j++)
        {
            ll c = a[j][i] * inv(a[rk + 1][i]) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[rk + 1][k]) % mod;
        }
        if(a[rk + 1][i] != 0)
        {
            rk++;
            bel[rk] = i;
        }
    }
    if(rk < n - 1)
    {
        for(int i = 1; i <= n; i++) cout << "0 ";
        exit(0);
    }
    gx = a[n][n + 1];
    A[gy] = 1;
    for(int i = n - 1; i >= 1; i--)
    {
        ll tmp = 0;
        for(int j = bel[i] + 1; j <= n; j++) tmp = (tmp - a[i][j] * A[j]) % mod;
        a[i][bel[i]] = (a[i][bel[i]] % mod + mod) % mod;
        A[bel[i]] = (tmp * inv(a[i][bel[i]]) % mod + mod) % mod;
    }
}
void GaussY(ll n, const ll a_[N][N])
{
    ll a[N][N];
    memcpy(a, a_, sizeof(a));
    for(int i = 1; i <= n; i++)
        for(int j = i + 1; j <= n; j++)
            swap(a[i][j], a[j][i]);
    ll rk = 0;
    ll tmpgy = 0;
    for(int i = 1; i <= n; i++)
    {
        for(int j = rk + 1; j <= n; j++)
        {
            if(a[j][i])
            {
                if((rk + 1) ^ j) swap(a[rk + 1], a[j]);
                break;
            }
        }
        if(a[rk + 1][i] == 0)
        {
            tmpgy = i;
            continue;
        }
        for(int j = rk + 2; j <= n; j++)
        {
            ll c = a[j][i] * inv(a[rk + 1][i]) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[rk + 1][k]) % mod;
        }
        if(a[rk + 1][i] != 0)
        {
            rk++;
            bel[rk] = i;
        }
    }
    B[tmpgy] = 1;
    for(int i = n - 1; i >= 1; i--)
    {
        ll tmp = 0;
        for(int j = bel[i] + 1; j <= n; j++) tmp = (tmp - a[i][j] * B[j]) % mod;
        a[i][bel[i]] = (a[i][bel[i]] % mod + mod) % mod;
        B[bel[i]] = (tmp * inv(a[i][bel[i]]) % mod + mod) % mod;
    }
}
ll Det(ll n, const ll a_[N][N])
{
    ll a[N][N];
    memcpy(a, a_, sizeof(a));
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(a[j][i])
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        if(a[i][i] == 0) continue;
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * inv(a[i][i]) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[i][k]) % mod;
        }
        res = (res * a[i][i]) % mod;
    }    
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
ll n, m, a[N][N], b[N][N], adj[N][N];
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    cin >> n >> m;
    for(int i = 1; i <= m; i++)
    {
        ll u, v, w;
        cin >> u >> v >> w;
        a[u][u] += w;
        a[u][v] -= w;
    }
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j <= n; j++)
            a[i][j] = (a[i][j] % mod + mod) % mod;   
        a[i][n + 1] = i;     
    }
    GaussX(n, a);
    GaussY(n, a);
    for(int i = 1, icnt = 0; i <= n; i++)
    {
        if(i == gx) continue;
        icnt++;
        for(int j = 1, jcnt = 0; j <= n; j++)
        {
            if(j == gy) continue;
            jcnt++;
            b[icnt][jcnt] = a[i][j];
        }
    }
    ll d = Det(n - 1, b);
    if((gx ^ gy) & 1) d = -d;
    d = (d % mod + mod) % mod;
    ll dk = d * inv(A[gy] * B[gx] % mod) % mod;
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            adj[i][j] = (A[i] * B[j] % mod * dk % mod);
    for(int i = 1; i <= n; i++) cout << adj[i][i] << " ";
    return 0;
}

6.10 P5807 【模板】BEST 定理 / Which Dreamed It

这就是上文中“循环同构不算同一个欧拉回路”的例子。直接套用 BEST 定理的公式计算即可。注意需要判断图中是否存在欧拉回路，并且有的点可能是孤立点。

时间复杂度 O(n^3)。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 105, V = 1e6 + 5;
const ll mod = 1e6 + 3;
ll n, a[N][N], b[N][N], degin[N], degout[N], inv[V], f[V];
int low[N], dfn[N], tot, cnt, scc[N], stk[N], tp;
bitset<N> instk;
vector<int> g[N];
void init()
{
    inv[1] = 1; f[0] = f[1] = 1;
    for(int i = 2; i < mod; i++)
    {
        f[i] = (f[i - 1] * i) % mod;
        inv[i] = (mod - mod / i) * inv[mod % i] % mod;
    }
}
void tarjan(int u)
{
    low[u] = dfn[u] = ++tot;
    stk[++tp] = u; instk[u] = 1;
    for(auto v : g[u])
    {
        if(!dfn[v])
        {
            tarjan(v);
            low[u] = min(low[u], low[v]);
        }
        else if(instk[v])
        {
            low[u] = min(low[u], dfn[v]);
        }
    }
    if(low[u] == dfn[u])
    {
        int x;
        ++cnt;
        do{
            x = stk[tp--];
            scc[x] = cnt;
            instk[x] = 0;
        } while(u != x);
    }
}
bool check()
{
    for(int i = 1; i <= n; i++)
        if(degin[i] ^ degout[i])
            return 0;
    for(int i = 1; i <= n; i++)
        if(!dfn[i])
            tarjan(i);   
    for(int i = 1; i <= n; i++)
    {
        if(degout[i] == 0) continue;
        if(scc[1] != scc[i]) return 0;
    } 
    return 1;
}
ll Det(ll n, ll a[N][N])
{
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(a[j][i])
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        if(!a[i][i]) return 0;
        a[i][i] = (a[i][i] % mod + mod) % mod;
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * inv[a[i][i]] % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[i][k]) % mod;
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
void solve()
{
    memset(a, 0, sizeof(a));
    memset(degin, 0, sizeof(degin));
    memset(degout, 0, sizeof(degout));
    memset(dfn, 0, sizeof(dfn));
    memset(low, 0, sizeof(low));
    memset(scc, 0, sizeof(scc));
    instk.reset();
    tp = tot = cnt = 0;
    cin >> n;
    for(int i = 1; i <= n; i++) g[i].clear();
    int smx = 0;
    for(int i = 1; i <= n; i++)
    {
        int x;
        cin >> x;
        smx += x;
        a[i][i] += x;
        degout[i] = x;
        if(x == 0) a[i][i] = 1, a[i][1] = -1;
        while(x--)
        {
            int v;
            cin >> v;
            a[i][v]--;
            degin[v]++;
            g[i].push_back(v);
        }
    }
    if(!check())
    {
        cout << "0\n";
        return;
    }
    if(smx == 0)
    {
        cout << "1\n";
        return;
    }
    for(int i = 1; i < n; i++)
        for(int j = 1; j < n; j++)
            b[i][j] = a[i + 1][j + 1];
    ll res = Det(n - 1, b);
    res = (res * f[degout[1]]) % mod;
    for(int i = 2; i <= n; i++)
    {
        if(degout[i] == degin[i] && degout[i] == 0) continue;
        res = (res * f[degout[i] - 1]) % mod;
    }
    cout << res << "\n";
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    init();
    int t;
    cin >> t;
    while(t--) solve();
    return 0;
}

6.11 P7531 [USACO21OPEN] Routing Schemes P

这有黑？这有黑？这有黑？这有黑？这有黑？

容易想到网络流里建立超级源点、超级汇点的思路，然后把超级汇点朝超级源点连一条边，那么问题就被转化为对有向图的欧拉回路计数，套用 BEST 定理即可解决。时间复杂度 O(Tn^3)。

注意因为欧拉回路中各路径的顺序并不会影响最后答案的计数，所以答案要除以 S!\times (S-1)!。其中 S 表示起点的个数。这个式子的含义是，除以 S! 代表去掉“从超级汇点走向超级源点选择不同边”的方案；而除以 (S-1)! 代表去掉从超级源点出发选择的不同路径的顺序，不是除以 S! 的原因是 BEST 定理已经去掉了循环同构的欧拉回路。

#include <bits/stdc++.h>
#define fi first
#define se second
#define eb(x) emplace_back(x)
#define pb(x) push_back(x)
#define lc(x) (tr[x].ls)
#define rc(x) (tr[x].rs)
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
typedef long double ldb;
typedef __int128 i128;
using pi = pair<int, int>;
const int N = 105;
const ll mod = 1e9 + 7;
int n, m, s, t;
ll a[N][N], b[N][N], f[N];
bitset<N> legal;
void init()
{
    f[0] = 1;
    for(int i = 1; i < N; i++)
        f[i] = (f[i - 1] * i) % mod;
}
void add(int u, int v, int w)
{
    a[u][u] += w;
    a[u][v] -= w;
}
ll qpow(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if(b & 1) res = (res * a) % mod;
        b >>= 1;
        a = (a * a) % mod;
    }
    return res;
}
ll Det(ll n, ll _a[N][N])
{
    ll a[N][N];
    memcpy(a, _a, sizeof(a));
    bool flag = 0; ll res = 1;
    auto Swap = [&] (int x, int y) -> void {
        swap(a[x], a[y]); flag ^= 1;
    };
    for(int i = 1; i <= n; i++)
    {
        for(int j = i; j <= n; j++)
        {
            if(a[j][i])
            {
                if(i ^ j) Swap(i, j);
                break;
            }
        }
        for(int j = i + 1; j <= n; j++)
        {
            ll c = a[j][i] * qpow(a[i][i], mod - 2) % mod;
            for(int k = i; k <= n; k++)
                a[j][k] = (a[j][k] - c * a[i][k]) % mod;
        }
        res = (res * a[i][i]) % mod;
    }
    if(flag) res = -res;
    return (res % mod + mod) % mod;
}
void solve()
{
    cin >> n >> m;
    s = n + 1; t = n + 2;
    memset(a, 0, sizeof(a));
    legal.reset();
    int cnts = 0;
    for(int i = 1; i <= n; i++)
    {
        char c;
        cin >> c;
        if(c == 'S') add(s, i, 1), cnts++;
        else if(c == 'R') add(i, t, 1);
    }
    add(t, s, cnts);
    for(int i = 1; i <= n; i++)
    {
        for(int j = 1; j <= n; j++)
        {
            char c;
            cin >> c;
            if(c - '0') add(i, j, 1);
        }
    }
    for(int i = 1; i <= n + 2; i++)
        if(a[i][i])
            legal[i] = 1;
    int tmpn = 0;    
    for(int i = 1, icnt = 0; i <= n + 2; i++)
    {
        if(!legal[i]) continue;
        icnt++; tmpn++;
        for(int j = 1, jcnt = 0; j <= n + 2; j++)
        {
            if(!legal[j]) continue;
            jcnt++;
            b[icnt][jcnt] = (a[i][j] % mod + mod) % mod;     
        }
    }
    ll res = Det(tmpn - 1, b);
    for(int i = 1; i <= n + 2; i++)
        if(a[i][i])
            res = (res * f[a[i][i] - 1]) % mod;
    res = (res * qpow(f[cnts], mod - 2)) % mod;
    res = (res * qpow(f[cnts - 1], mod - 2)) % mod;
    cout << res << "\n";
}
int main()
{
    //freopen("sample.in", "r", stdin);
    //freopen("sample.out", "w", stdout);
    ios::sync_with_stdio(0);
    cin.tie(0);
    cout.tie(0);
    init();
    int T;
    cin >> T;
    while(T--) solve();
    return 0;
}

参考资料