学习笔记:Matrix-Tree 定理

· · 算法·理论

往期回顾:

学习笔记:状压 DP
学习笔记:数位 DP

Matrix-Tree 定理,矩阵的又一个绝妙应用。

::::align{right} ——前言 ::::

Matrix-Tree 定理是一个比较困难的东西,主要难点在于理解模板,以及其背后的原理。虽然我们的老师不谈证明,但是我认为理解一下为什么是比较重要的,所以我会给出简要的证明。

1. 前置知识:行列式

1.1 定义

行列式是一种建立在方阵(即长宽相等的矩阵)上的运算,对于矩阵 A,我们通常记其行列式为 |A| 或者 \det(A)

行列式的定义式是一个比较复杂的式子,长这样:

\det(A)=\sum_{p}(-1)^{\tau(p)}\prod_{i=1}^{n} A_{i,p_i}

其中求和符号表示枚举 1\sim n 的全排列,p 就是当前枚举的排列。\tau(p) 表示 p 的逆序对个数。容易发现,按照定义式直接做行列式求值是 O(n\log n\times n!) 的,可能可以精细实现做到 O(n\times n!),但是这个方法太劣了,需要优化。

1.2 性质

1.2.1 性质 1

行列式有着一些巧妙的性质。为了探究这些性质,我们先从排列的奇偶性入手,因为排列的奇偶性决定了后面那一堆乘积的正负号。奇偶性的定义如下:

如果 \tau(p) 为奇数,则称 p 为一个奇排列,否则为一个偶排列。

然后,奇偶排列是可以互相转换的。有以下的定理:

交换排列的两个数,排列的奇偶性改变。 ::::info[证明] 先从简单情况入手。假设交换的是相邻两个数,那么排列的逆序对数量必定 +1 或者 -1,那么奇偶性必定会改变。

接着想一般情况怎么证。交换离的比较远的两个数(设位置为 ij,其中 i<j),我们可以看做是在做下面的交换:先把 j 向左交换 j-i 次,再把 i 向右交换 j-i-1 次,那么总次数就是 2(j-i)-1 次。这是一个奇数。奇偶性改变奇数次,那么就等同于奇偶性改变一次。 ::::

于是就可以证行列式的第一个性质了。

交换行列式的两行,行列式变号。即下面式子:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}=-\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}

::::info[证明] 首先,我们注意到后面的 \prod^{n}_{i=1}A_{i,p_i} 这里,对于第 m 行和第 k 行都只会贡献一次。于是,这个排列后面乘上去的取值不变,因为交换了之后只是乘的顺序变了。

在前面的 (-1)^{\tau(p)},由于交换了两个数,所以排列奇偶性改变,也就是行列式整体变号。于是得证。 ::::

1.2.2 性质 2

接着,还有以下的一些性质:

每一行乘上一个系数,行列式整体乘以这个系数。即:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ ka_{m,1}&ka_{m,2}&\ldots&ka_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}=k\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}

::::info[证明] 依旧注意到 \prod^{n}_{i=1}A_{i,p_i} 对于每一行只会贡献一次,于是我们可以把 k 提出来,提到最外面你就发现是 k 倍后面矩阵的行列式了。形式化地,

\begin{aligned}&\sum_{p}(-1)^{\tau(p)}\prod^{n}_{i=1}A_{i,p_i}\\=&\sum_{p}(-1)^{\tau(p)}A_{m,p_m}\prod^{n}_{i=1,i\neq m}A_{i,p_i}\\=&\sum_{p}(-1)^{\tau(p)}ka_{m,p_m}\prod^{n}_{i=1,i\neq m}A_{i,p_i}\\=&k\sum_{p}(-1)^{\tau(p)}\prod^{n}_{i=1}a_{i,p_i}\end{aligned}

::::

1.2.3 性质 3

如果有两行的值相同,那么这个矩阵的行列式为 0。即下面式子:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}=0

::::info[证明] 我们可以交换这两行,所以行列式的符号改变;但是交换之后行列式又没有实质性的改变,所以行列式值不变。一个数等于其相反数,行列式的值只能是 0。 ::::

1.2.4 性质 4

矩阵的一行加上某个序列,行列式等于原矩阵的行列式加上把那一行替换为序列之后矩阵的行列式。也就是下面这个式子:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ b_1+a_{k,1}&b_2+a_{k,2}&\ldots&b_3+a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}=\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}+\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ b_1&b_2&\ldots&b_3 \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}

::::info[证明] 还是用到最基本的东西:后面的求积式中每一行只会贡献一遍。于是可以使用乘法分配律拆开,整体拆开成两个部分,也就得到了后面两个矩阵的行列式。形式化地,

\begin{aligned}&\sum_{p}(-1)^{\tau(p)}\prod^{n}_{i=1}A_{i,p_i}\\=&\sum_{p}(-1)^{\tau(p)}(b_{p_k}+a_{k,p_k})\prod^{n}_{i=1,i\neq k}A_{i,p_i}\\=&\sum_{p}(-1)^{\tau(p)}b_{p_k}\prod^{n}_{i=1,i\neq k}A_{i,p_i}+\sum_{p}(-1)^{\tau(p)}a_{k,p_k}\prod^{n}_{i=1,i\neq k}A_{i,p_i}\end{aligned}

::::

1.2.5 性质 5

将行列式的某一行加上另一行的某一倍,行列式值不变。即:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ xa_{m,1}+a_{k,1}&xa_{m,2}+a_{k,2}&\ldots&xa_{m,n}+a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}=\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}

::::info[证明] 使用之前的性质,变形即可。

\begin{aligned}&\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ xa_{m,1}+a_{k,1}&xa_{m,2}+a_{k,2}&\ldots&xa_{m,n}+a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}\\=&\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ xa_{m,1}&xa_{m,2}&\ldots&xa_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}+\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}\\=&x\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}+\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}\\=&x\times 0+\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\ a_{2,1}&a_{2,2}&\ldots&a_{2,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{m,1}&a_{m,2}&\ldots&a_{m,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{k,1}&a_{k,2}&\ldots&a_{k,n} \\ \vdots&\vdots&\ddots&\vdots \\ a_{n,1}&a_{n,2}&\ldots&a_{n,n}\end{pmatrix}\end{aligned}

于是得证。 ::::

1.2.6 性质 6

对于一个上三角矩阵,其行列式为主对角线元素的乘积。即:

\det\begin{pmatrix}a_{1,1}&a_{1,2}&\ldots&a_{1,n}\\0&a_{2,2}&\ldots&a_{2,n}\\\vdots&\vdots&\ddots&\vdots\\0&0&\ldots&a_{n,n}\end{pmatrix}=\prod^{n}_{i=1}a_{i,i}

::::info[证明] 注意到只有 1,2,3,\ldots,n 这一个排列才对于总行列式有贡献,且逆序对为 0,所以贡献是正的。这一个排列的贡献就是主对角线元素的乘积。

如果是其他排列,显然会选到 0,乘起来就没有贡献了。 ::::

1.3 求法

知道了上面这些性质,求出行列式也就很容易了。注意到我们在上面提供了这些操作:交换两行,一行加上另一行的倍数,以及最后对于上三角矩阵的处理方法。

于是我们考虑将矩阵转化为上三角矩阵,那么具体该怎么做呢?显然是高斯消元啊!

注意到一件事:行列式的值一般很大,所以需要取模。假如模数是质数,可以直接使用费马小定理求逆元,但是假如模数不是质数,就可能没有逆元,你不炸了

其实还是可以做的。我们考虑类似欧几里得算法的辗转相除,我们每次将一个式子的主元和另一个式子的主元算出系数(不用 double,向下取整),然后两个式子相减,最后可以反复相减得到答案。只不过这种方法好像常数大一点,跑不过费马小定理求逆元。

接下来,我们证明辗转相除法的时间复杂度是 O(n^3) 的。

::::info[证明] 高斯消元本身是 O(n^3) 的。所以我们考虑辗转相除部分一共会除几次。

有一个很显然的性质,就是 ab 求余(b<a),得到的结果必定小于 \frac{a}{2}。因为假如 b<\frac{a}{2},结果小于 b,所以结果小于 \frac{a}{2}。假如 b>\frac{a}{2},那么结果为 a-b,这个是小于 \frac{a}{2} 的。如果 b=\frac{a}{2},那结果不就是 0 吗……

因此,每一次辗转相除,都使用 O(n) 的时间复杂度让 O(n) 个数字变成了原来的 \frac{1}{2} 以下。总共有 O(n^2) 个数字,每一个数字只需要 O(\log V) 次辗转相除就会变成 1,所以这一部分至多是 O(n^2\log V) 的,不是主要时间复杂度。

故取高斯消元时间复杂度,为 O(n^3)。 ::::

最后,使用模板题给出代码,由于这道题是任意模数的,所以只能辗转相除:

::::success[code]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int N = 600 + 5;
int n, mod, a[N][N];

int det(int a[][N], int mod) {
    int w = 1, res = 1;
    for(int i = 1; i <= n; i ++ ) {
        for(int j = i + 1; j <= n; j ++ ) {
            while(a[i][i]) { //辗转相除
                int d = a[j][i] / a[i][i]; //系数
                for(int k = i; k <= n; k ++ ) a[j][k] = (a[j][k] - d * a[i][k] % mod + mod) % mod;
                swap(a[j], a[i]); //交换两行
                w = -w; //行列式要变号
            }
            swap(a[j], a[i]); //最后要换回去
            w = -w;
        }
    }
    for(int i = 1; i <= n; i ++ ) res = res * a[i][i] % mod; //最后答案就是主对角线元素的乘积
    return (res * w + mod) % mod; //还要乘上符号
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n >> mod;
    for(int i = 1; i <= n; i ++ ) for(int j = 1; j <= n; j ++ ) cin >> a[i][j];
    cout << det(a, mod) << endl;
    return 0;
}

::::

由于这不是重点,所以没有例题。只不过我又翻出来了两道模板:https://www.luogu.com.cn/problem/UVA684,https://www.luogu.com.cn/problem/SP2832。

2. Matrix-Tree 定理

2.1 理论

Matrix-Tree 定理,主要用于解决图的生成树计数问题。它可以直接解决的问题是什么呢?参见模板题。但是,模板题稍微有一点困难,我们考虑弱化问题:

给定一个无向图 G,求其生成树的个数。

矩阵树定理提供了一种美妙的,通过行列式解决这种问题的方法。矩阵树定理是什么?其实无外乎以下一句话:

图的 kirchhoff 矩阵去掉最后一行和最后一列的行列式,即为该图的生成树个数。

其中,kirchhoff 矩阵定义为图的度数矩阵(对角线上为该点的入度,非对角线为 0)减去图的邻接矩阵。比如说,假如我们有图 G 的边集为 \{(1,2),(2,3),(1,3)\},那么其 kirchhoff 矩阵如下:

K=\begin{bmatrix}2&0&0\\0&2&0\\0&0&2\end{bmatrix}-\begin{bmatrix}0&1&1\\1&0&1\\1&1&0\end{bmatrix}=\begin{bmatrix}2&-1&-1\\-1&2&-1\\-1&-1&2\end{bmatrix}

所以我们要求的是 \begin{bmatrix}2&-1\\-1&2\end{bmatrix} 的行列式。

接下来,就是对于矩阵树定理困难的证明了!(其实你直接记结论没人会骂你的)

我们考虑钦定一个点为根,并且把这一个点对应的那一行那一列换到 kirchhoff 矩阵的最后一行,最后一列,为了处理无向图方便,我们可以直接让 n 为根就不用换了。

对于度数矩阵的情况,就是每一个节点选择正好一个出边的方案数了。但是呢。选择正好一个出边,可能连出来的并不是一个树,可能会连出来一些基环树。

于是,减去邻接矩阵,其实就相当于一个容斥(就是容斥公式……)。可以做如下的思考(注意后面讨论的全是加上邻接矩阵):

于是,我们注意到了另一个重要的东西:一个排列,可以被表示为一个置换,即使用若干条边连接 ip_i,得到的若干个环。对于 kirchhoff 矩阵,这个排列正好对应的就是 \prod^{n}_{i=1}K_{i,p_i}

于是我们带上形成置换环的个数,记为 c(p),作为 -1 的指数(容斥系数),最后的式子就是这个:

\sum_{p}(-1)^{c(p)}\prod^{n}_{i=1}K_{i,p_i}

注意到这个式子长得和行列式的定义式很像啊!只不过逆序对数量变成了置换环的数量,所以这东西到底是否等价于行列式呢?

不是等价的,但是可以转化。对于每一个环,我们都可以交换环长减一个数字使其有序,所以一个排列最多的交换次数是 n-c(p)。最后排列有序了,所以逆序对数量为 0,环长为 1。每一次交换两个数的时候,排列奇偶性会改变(在之前的行列式板块已经证明),所以排列奇偶性和置换环奇偶性是相反的。

于是最后容斥的时候就减去邻接矩阵,而不是加上邻接矩阵,这样就可以再取一次反,就等价啦!

至此,Matrix-Tree 定理得证。剩下的就是一个行列式求值板子了,可以容易的做到 O(n^3)

但是,模板题是有向图啊!其实并不困难,注意到容斥的意义,就是去掉形成环的情况。由于是外向树,我们只需要改成去掉相反顺序的环,也就是容斥的时候邻接矩阵改为内向的边数就可以了。假如是内向树,我们就将邻接矩阵改为外向的边数。注意这种情况一定要把根的行列换到外边。

对于边有权的情况,可以使用以下的方法思考:将权值当做重边数量,那么对于同一个形态的生成树,其数量为各个边的重边条数之积。(显然的乘法原理)于是,这个正好符合题意,直接把边权加到邻接矩阵就好啦~

2.2 代码实现

由于这道题模数是质数,所以使用辗转相除和费马小定理都可以求行列式。所以,我给出两种方法的实现。后面的代码我都只写辗转相除了,因为更通用。但是费马小定理要更快一点,大概快一倍,卡常的时候可以一试。

::::success[辗转相除]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int N = 300 + 5;
const int MOD = 1e9 + 7;
int n, m, t, a[N][N];

int det(int a[N][N], int n) {
    int w = 1, res = 1;
    for(int i = 1; i <= n; i ++ ) {
        for(int j = i + 1; j <= n; j ++ ) {
            while(a[i][i]) {
                int d = a[j][i] / a[i][i];
                for(int k = i; k <= n; k ++ ) a[j][k] = (a[j][k] - d * a[i][k] % MOD + MOD) % MOD;
                swap(a[i], a[j]);
                w = -w;
            }
            swap(a[i], a[j]);
            w = -w;
        }
    }
    for(int i = 1; i <= n; i ++ ) res = res * a[i][i] % MOD;
    return (res * w + MOD) % MOD; 
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n >> m >> t;
    for(int i = 1; i <= m; i ++ ) {
        int u, v, w;
        cin >> u >> v >> w; //构造 kirchhoff 矩阵
        if(t == 0) {
            a[u][v] = (a[u][v] + MOD - w) % MOD;
            a[v][u] = (a[v][u] + MOD - w) % MOD;
            a[u][u] = (a[u][u] + w) % MOD;
            a[v][v] = (a[v][v] + w) % MOD; 
        } else {
            a[u][v] = (a[u][v] + MOD - w) % MOD; 
            a[v][v] = (a[v][v] + w) % MOD;
        }
    }
    if(t == 1) { //把根换到外面
        swap(a[1], a[n]);
        for(int i = 1; i <= n; i ++ ) swap(a[i][1], a[i][n]);
    }
    cout << det(a, n - 1) << endl; 
    return 0; 
}

::::

::::success[费马小定理]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int MOD = 1e9 + 7;
const int N = 300 + 5;
int n, m, t;
int *a[N], _a[N][N];

int qpow(int a, int b = MOD - 2) {
    int res = 1;
    for(; b > 0; b >>= 1, a = (a * a) % MOD)
        if(b & 1) res = (res * a) % MOD;
    return res;
}

int det(int **a) {
    int ans = 1;
    bool tr = 0;
    for(int j = 1; j < n; j ++ ) {
        for(int i = j; i < n; i ++ ) {
            if(a[i][j]) {
                if(i != j) {
                    swap(a[i], a[j]);
                    tr ^= 1;
                }
                break;
            }
        }
        if(a[j][j] == 0) return 0;
        ans = ans * a[j][j] % MOD;
        int T = qpow(a[j][j]);
        for(int k = j; k < n; k ++ ) a[j][k] = T * a[j][k] % MOD;
        for(int i = j + 1; i < n; i ++ ) {
            T = MOD - a[i][j];
            for(int k = j; k < n; k ++ ) a[i][k] = (a[i][k] + T * a[j][k]) % MOD;
        }
    }
    return tr ? (MOD - ans) % MOD : ans;
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n >> m >> t;
    for(int i = 1; i <= n; i ++ ) a[i] = _a[i];
    for(int i = 1, u, v, w; i <= m; i ++ ) {
        cin >> u >> v >> w;
        if(u == v) continue;
        if(t == 1) {
            a[u][v] = (a[u][v] - w + MOD) % MOD;
            a[v][v] = (a[v][v] + w) % MOD;
        } else {
            a[u][v] = (a[u][v] - w + MOD) % MOD;
            a[v][u] = (a[v][u] - w + MOD) % MOD;
            a[u][u] = (a[u][u] + w) % MOD;
            a[v][v] = (a[v][v] + w) % MOD;
        }
    }
    if(t == 1) {
        swap(a[1], a[n]);
        for(int i = 1; i <= n; i ++ ) swap(a[i][1], a[i][n]);
    }
    cout << det(a) << endl;
    return 0;
}

::::

2.3 例题

2.3.1 那些板子

太板了不用讲,但是多水几道紫题不是好事吗?建议做这些题的时候每一道题都写一遍板子,不要直接照抄。

2.3.2 重建

https://www.luogu.com.cn/problem/P3317

给定一个无向图,每一条边有边权 w,表示这条边存在的概率为 w。求这个图是一个树的概率。

哎不对,矩阵树定理不是用来算计数的吗?怎么开始算期望了?等等,这道题本质就是计数。我们先用数学语言形式化的表示矩阵树定理算的东西和这道题要算的东西。首先,【模板】Matrix-Tree 定理这道题求得是下面一个式子:

\sum_{T}\prod_{e\in T}p_e

也就是所有生成树边权的乘积之和。而这道题要算什么呢?可以理解为每一种生成树的概率相加,于是就是下面的式子:

\sum_{T}\prod_{e\in T}p_e\prod_{e\notin T}(1-p_e)

于是尝试变形化为矩阵树定理的标准形式。

\begin{aligned}&\sum_{T}\prod_{e\in T}p_e\prod_{e\notin T}(1-p_e)\\=&\sum_{T}\prod_{e\in T}p_e\dfrac{\prod(1-p_e)}{\prod_{e\in T}(1-pe)}\\=&\prod(1-p_e)\sum_{T}\prod_{e\in T}\dfrac{p_e}{1-p_e}\end{aligned}

于是把边权规定为 \frac{p_e}{1-p_e} 跑一遍矩阵树定理就可以了。注意 p_e=1 的时候会除 0,所以我们把 1-p_e 加上 eps 即可。注意本题有一点卡精度。求行列式的时候,直接除就可以了。

时间复杂度显然 O(n^3)

::::success[code]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int N = 50 + 5;
const long double eps = 1e-12;
long double g[N][N], a[N][N], res = 1;
int n;

long double det() {
    int w = 1;
    long double res = 1;
    for(int i = 1; i < n; i ++ ) {
        int mx = i;
        for(int j = i + 1; j < n; j ++ )
            if(fabs(a[j][i]) > fabs(a[mx][i])) mx = j;
        if(mx != i) swap(a[i], a[mx]);
        if(fabs(a[i][i]) < eps) return 0;
        for(int j = i + 1; j < n; j ++ ) {
            double t = a[j][i] / a[i][i];
            for(int k = i; k < n; k ++ ) a[j][k] -= a[i][k] * t;
        }
        res *= a[i][i];
    }
    return res;
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n;
    for(int i = 1; i <= n; i ++ ) for(int j = 1; j <= n; j ++ ) cin >> g[i][j];
    for(int i = 1; i <= n; i ++ ) for(int j = 1; j <= n; j ++ ) {
        if(i != j) {
            a[i][j] -= g[i][j] / (1 - g[i][j] + eps);
            a[i][i] += g[i][j] / (1 - g[i][j] + eps);
        }
        if(i < j) res *= (1 - g[i][j] + eps);
    }
    cout << res * det() << endl;
    return 0;
}

::::

2.3.3 巨额奖金

https://www.luogu.com.cn/problem/P2143

给定一个无向图,求其最小生成树的个数。

这个看着就很矩阵树定理啊,但是具体怎么处理呢?最小生成树有一个很重要的性质,不知道这个性质这题完全没法做。

对于一个图的所有最小生成树,边权相等的边的数量相等。

::::info[证明] 我们考虑 kruskal 算法的过程。对于一些权值相同的边,kruskal 算法肯定是把能够合并的连通块全部都合并了,能合并的连通块数量肯定是相等的,所以这种边的选择数量肯定是相等的,并且是在保证是一棵树的情况下最大的。 ::::

于是,我们就可以做如下的处理:先随便跑出一个最小生成树,然后讨论每一种权值的边的替换次数。这个可以怎么做呢?矩阵树定理。

首先,其它权值的边肯定都是不会变的,并且我们添加的边不能形成环,所以我们可以先把其它权值的点缩到一起,这样形成环就变为了自环,生成树计数就不会考虑到了。

对于权值相等的边,会形成一个图。因为我们要连接所有的连通块,所以就会形成一个生成树。直接在这个缩点后的图上跑生成树计数就可以了。最后对于每一种权值的边的替换情况,可以使用乘法原理,直接统计即可。

时间复杂度乍一看好像是 O(mn^3) 的,过不了。但是下面会证明时间复杂度是 O(n^3) 的。

我们思考一下每一次矩阵树定理都会对于多少个点跑。我们令这些数字为一个序列 c,那么可以进行如下分析:

O\left(\sum c_i^3\right)=O\left(\sum(c_i-1)^3\right)\leq O\left(\left(\sum c_i-1\right)^3\right)=O(n^3)

::::success[代码]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int N = 100 + 5;
const int M = 1000 + 5;
const int MOD = 31011;
int n, m, dt[N], cnt, a[N][N], ans = 1;
struct Edge { int u, v, w; bool ont; } e[M];
vector<int> wsy;

struct dsu {
    int k[N];
    void init() { for(int i = 1; i <= n; i ++ ) k[i] = i; }
    int find(int i) { return k[i] == i ? i : k[i] = find(k[i]); }
    void merge(int x, int y) { k[find(x)] = find(y); }
} u;

bool cmp(Edge a, Edge b) { return a.w < b.w; }
void kruskal() {
    sort(e + 1, e + m + 1, cmp);
    u.init();
    int cnt = 0;
    for(int i = 1; i <= m && cnt <= n - 2; i ++ ) {
        if(u.find(e[i].u) != u.find(e[i].v)) {
            u.merge(e[i].u, e[i].v);
            wsy.push_back(e[i].w);
            e[i].ont = true;
            cnt ++;
        }
    }
    sort(wsy.begin(), wsy.end());
    wsy.erase(unique(wsy.begin(), wsy.end()), wsy.end());
}

int det(int a[N][N]) {
    int w = 1, res = 1;
    for(int i = 1; i < cnt; i ++ ) {
        for(int j = i + 1; j < cnt; j ++ ) {
            while(a[i][i]) {
                int d = a[j][i] / a[i][i];
                for(int k = i; k <= n; k ++ ) a[j][k] = (a[j][k] - d * a[i][k] % MOD + MOD) % MOD;
                swap(a[j], a[i]);
                w = -w;
            }
            swap(a[j], a[i]);
            w = -w;
        }
    }
    for(int i = 1; i < cnt; i ++ ) res = res * a[i][i] % MOD;
    return (res * w + MOD) % MOD; 
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n >> m;
    for(int i = 1; i <= m; i ++ ) {
        cin >> e[i].u >> e[i].v >> e[i].w;
        e[i].ont = false;
    }
    kruskal();
    for(auto w : wsy) {
        u.init();
        memset(dt, 0, sizeof dt);
        memset(a, 0, sizeof a);
        cnt = 0;
        for(int i = 1; i <= m; i ++ )
            if(e[i].ont && e[i].w != w) u.merge(e[i].u, e[i].v);
        for(int i = 1; i <= n; i ++ ) {
            int f = u.find(i);
            if(!dt[f]) dt[f] = ++ cnt;
        }
        for(int i = 1; i <= m; i ++ )
            if(e[i].w == w) {
                int uu = dt[u.find(e[i].u)], v = dt[u.find(e[i].v)];
                a[uu][v] = (a[uu][v] - 1 + MOD) % MOD;
                a[v][uu] = (a[v][uu] - 1 + MOD) % MOD;
                a[uu][uu] = (a[uu][uu] + 1) % MOD;
                a[v][v] = (a[v][v] + 1) % MOD;
            }
        ans = (ans * det(a)) % MOD;
    }
    cout << ans << endl;
    return 0;
}

::::

有双倍经验。不是,JSOI 搬两年前的题是 hyw。那道题由于边权相同的边不超过 10 条,所以直接爆搜做生成树计数也没有问题。

2.3.4 黑暗前的幻想乡

https://www.luogu.com.cn/problem/P4336

给定 n-1 个边集,从每个边集中各选一条边,求可能形成的生成树个数。

这道题 n\leq 17 很小啊,考虑指数级做法。首先呢,生成树一共有 n-1 条边,其中每个边集各选一条边,那么换言之,就是假如有一个边集没有选,就不合法。我们需要统计的就是所有边集都选了的方案数。

考虑容斥。其实只要想到容斥这道题就非常简单了。 我们可以枚举所有的选择公司情况,算出来选择这些公司的方案数。(不要求每个公司都要选择)

最后按照标准容斥,直接做就可以了。容斥的系数是 (-1)^x,其中 x 表示未选择公司的个数。这个容斥系数很好理解,我感觉容斥系数我一般都是这么推的,可能是因为我现在做的容斥题都太简单了:

具体实现可以使用搜索或者状压,我选择的是状压,因为常数更小。

时间复杂度 O(n^32^n)

::::success[code]

#include <bits/stdc++.h>
#define int long long
#define endl '\n'
using namespace std;

const int N = 17 + 5;
const int MOD = 1e9 + 7;
int n, m[N], a[N][N], ans;
vector<pair<int, int> > q[N];

int det(int a[N][N], int n) {
    int w = 1, res = 1;
    for(int i = 1; i <= n; i ++ ) {
        for(int j = i + 1; j <= n; j ++ ) {
            while(a[i][i]) {
                int d = a[j][i] / a[i][i];
                for(int k = i; k <= n; k ++ ) a[j][k] = (a[j][k] - d * a[i][k] % MOD + MOD) % MOD;
                swap(a[i], a[j]);
                w = -w;
            }
            swap(a[i], a[j]);
            w = -w;
        }
    }
    for(int i = 1; i <= n; i ++ ) res = res * a[i][i] % MOD;
    return (res * w + MOD) % MOD;
}

signed main() {
    ios :: sync_with_stdio(0);
    cin.tie(0); cout.tie(0);
    cin >> n;
    for(int i = 1; i < n; i ++ ) {
        cin >> m[i];
        for(int j = 1; j <= m[i]; j ++ ) {
            int u, v;
            cin >> u >> v;
            q[i].emplace_back(u, v);
        }
    }
    for(int i = 0; i < (1 << (n - 1)); i ++ ) {
        int cnt = n;
        memset(a, 0, sizeof a);
        for(int j = 1; j < n; j ++ ) {
            if(!(i & (1 << (j - 1)))) continue;
            cnt -- ;
            for(auto [u, v] : q[j]) {
                a[u][u] = (a[u][u] + 1) % MOD;
                a[v][v] = (a[v][v] + 1) % MOD;
                a[u][v] = (a[u][v] - 1 + MOD) % MOD;
                a[v][u] = (a[v][u] - 1 + MOD) % MOD;
            }
        }
        ans = (ans + ((cnt & 1) ? 1 : -1) * det(a, n - 1) + MOD) % MOD; 
    }
    cout << ans << endl;
    return 0; 
}

::::

Matrix-Tree 定理,如同桥梁一般,把图论和矩阵两个看似毫不相干的领域连接了起来。有的时候,两个领域看似毫不相干,但是结合起来,便会产生新的,美妙的东西。OI 之中,这种东西从来都不会少。一道难题看似不可攻克,但是,它毕竟是无数个琐碎的知识点拼起来的。

::::align{right} ——后记 ::::