Matrix-Tree是无向图生成树计数的利器，本文需要线性代数的前缀知识。
        另外：第100篇博客纪念XD
        UPD: 修改了一些个人认为比较严重的问题。

        给定一个图，如何统计这个图生成树的数目？这就是Matrix-Tree解决的问题。
        在了解这个定理之前，先介绍高斯消元。

高斯消元

        我们都知道矩阵可以用来求解线性方程组。但是对于一个给定的矩阵，如何去求解？高斯消元法利用了矩阵初等行变换的性质，将矩阵化为行阶梯型矩阵，再求解方程组。

        对于上面的矩阵，可以将其化为下面的形式：

        进行这样操作的算法就是高斯消元法，该算法时间复杂度为$O(n^3)$，对于模板题，它的算法代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
inline void Guass() {
    for (int i = 1; i <= n; i++) {//枚举列数，总列数为n
        int maxn = i;//假设当前行（第i列对应第i行）
        for (int j = i + 1; j <= n; j++)if (fabs(op[j][i]) > fabs(op[maxn][i]))maxn = j;//找向下的行中对应列元素绝对值最大的一个，绝对值大有利于减小误差
        if (i != maxn)swap(op[i], op[maxn]);//不同的话，交换两行，swap可以直接交换两个数组
        if (fabs(op[maxn][i]) < EPS) {//绝对值最大的为0，说明不满秩，方程组没有唯一解，结束程序。EPS用来消浮点误差。
            cout << "No Solution";
            exit(0);
        }
        double div = op[i][i];//找到除法因子
        for (int j = i; j <= n + 1; j++)op[i][j] /= div;//这一行先全部除以该值，此时首元素应为1
        for (int j = i + 1; j <= n; j++) {//其余行减掉该行的倍数
            div = op[j][i];
            for (int z = 1; z <= n + 1; z++)op[j][z] -= op[i][z] * div;
        }
    }
}

        高斯消元就是这样了，原理上的东西见线性代数教材。注意这里的高斯消元会将主对角线上的元素同一化为1。

Kirchhoff矩阵

        介绍Kirchhoff矩阵之前需要先了解邻接矩阵和度数矩阵。
        邻接矩阵中，$A[i][j]$表示结点i和结点j之间边的数目；度数矩阵中若$i≠j$则$D[i][j]=0$，否则为点i的度数（即连接该点的边的数目）。Kirchhoff矩阵定义为$K=D-A$。

        如上图，其邻接矩阵和度数矩阵分别为：

        于是可以得到Kirchhoff矩阵：

        Kirchhoff矩阵有一些性质如下：

  • 性质一：对于任意无向图，其Kirchhoff矩阵行列式必为0即有$|K|=0$
            易知矩阵K中每一行元素之和必为0，因此行列式为0。

  • 性质二：对于不连通的无向图，其Kirchhoff矩阵的任意n-1阶主子式必为0
            首先认识到调整行列式中行的顺序不会影响行列式是否为0这个问题的结果，因此可以调整行列式中行的顺序关系，使属于同一个连通分量中的点在相邻的两行。这样$|K|$就成为下面的对角形式：

            图不连通时，去除一行一列后再求n-1阶主子式，必然使得其中至少一个连通分量是完整的。根据性质一，其行列式为0，因此乘积为0。

  • 性质三：如果图是一棵树，则其Kirchhoff矩阵的任意n-1阶主子式必为1
            证明略。

Matrix-Tree（矩阵树）定理

        【Matrix-Tree】矩阵树定理：对于一个无向图，其生成树数目为其Kirchhoff矩阵任意n-1阶主子式。
        矩阵树定理给出了生成树计数的方法。那么如何求行列式？一种方法是将行列式化为上三角形式，然后主对角线上的元素的乘积即为行列式的值，对角化过程可以用高斯消元实现。
        看一个模板题。很明显是生成树计数问题，直接用矩阵树定理即可。注意要用高精度计算，这里用python配decimal消精度误差偷个懒。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#!/usr/bin/env python3
from decimal import Decimal, getcontext

getcontext().prec = 45
n = int(input())
if n == 1 or n == 2:
    if n == 1:
        print(1)
    else:
        print(5)
    exit(0)
op = []
ans = Decimal('1')
for i in range(0, n): # 构造Kirchhoff矩阵
    op.append([])
    for j in range(0, n + 1):
        if i == j:
            op[i].append(Decimal('3'))
        elif j == (i - 1 + n) % n or j == (i + 1 + n) % n:
            op[i].append(Decimal('-1'))
        else:
            op[i].append(Decimal('0'))
for i in range(0, n): # 高斯消元
    p = i
    for j in range(i + 1, n):
        if op[j][i] > op[p][i]:
            p = j
    if p != i:
        op[p], op[i] = op[i], op[p]  # Swap
    for j in range(i + 1, n): # 这里和上面高斯消元不同，省略了当前行首元素化为1的步骤
        div = op[j][i] / op[i][i]
        for z in range(i, n):
            op[j][z] -= div * op[i][z]
for i in range(0, n):
    ans *= op[i][i]
print("{:.0f}".format(ans))

最小生成树计数

        矩阵树可以求生成树的数量，但是如果求最小生成树又如何求呢？
        首先需要知道这样一个性质：最小生成树中，相同权值的边的数量是一定的。根据这个性质，我们先建出一个最小生成树，然后枚举权值，从树上剔除这些相同权值的边，可以得到若干联通块。进而将这些联通块看成点，跑一次矩阵树就能找出它的生成树数量，易知这些都是最小生成树。枚举边权后，根据乘法原理就可以求出方案数。
        模板。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#include <bits/stdc++.h>

#define MOD 31011
using namespace std;

struct Edge0 {
    int t1, t2, v;

    bool operator<(Edge0 s) {
        return v < s.v;
    }
} e[1005];

struct Edge {
    int to, v, next;
} edge[205];

double op[105][105];
int n, m, fa[105], id[105], cnt = 1, head[105], ID, now, ans = 1;
set<int> ssp;

inline void add(int x, int y, int v) {
    edge[cnt].v = v, edge[cnt].to = y, edge[cnt].next = head[x], head[x] = cnt++;
    edge[cnt].v = v, edge[cnt].to = x, edge[cnt].next = head[y], head[y] = cnt++;
}

int findF(int x) {
    return fa[x] == x ? x : fa[x] = findF(fa[x]);
}

inline double Guass(int n) {
    double a = 1;
    for (int i = 1; i <= n; i++) {
        int maxn = i;
        for (int j = i + 1; j <= n; j++)
            if (fabs(op[j][i]) > fabs(op[maxn][i]))maxn = j;
        if (i != maxn)swap(op[i], op[maxn]);
        double div = op[i][i];
        a *= div;
        for (int j = i; j <= n; j++)op[i][j] = op[i][j] / div;
        for (int j = i + 1; j <= n; j++) {
            div = op[j][i];
            for (int z = 1; z <= n; z++)op[j][z] -= op[i][z] * div;
        }
    }
    return a;
}

void DFS(int x) {
    id[x] = ID;
    for (int i = head[x]; i; i = edge[i].next) {
        if (id[edge[i].to] == 0 && edge[i].v != now)DFS(edge[i].to);
    }
}

int main() {
    scanf("%d%d", &n, &m);
    for (int i = 1; i <= n; i++)fa[i] = i;
    for (int i = 1; i <= m; i++)scanf("%d%d%d", &e[i].t1, &e[i].t2, &e[i].v);
    sort(e + 1, e + m + 1);
    for (int i = 1; i <= m; i++) {
        if (findF(e[i].t1) != findF(e[i].t2)) {
            fa[findF(e[i].t1)] = findF(e[i].t2), ssp.insert(e[i].v), add(e[i].t1, e[i].t2, e[i].v);
        }
    }
    for (int i = 2; i <= n; i++)
        if (findF(i) != findF(i - 1)) {
            cout << 0;
            return 0;
        }
    for (int it : ssp) {
        memset(id, 0, sizeof(id)), ID = 0, now = it;
        for (int i = 1; i <= n; i++)for (int j = 1; j <= n; j++)op[i][j] = 0;
        for (int i = 1; i <= n; i++)if (id[i] == 0)++ID, DFS(i);
        for (int i = 1; i <= m; i++) {
            if (e[i].v == now) {
                if (id[e[i].t1] == id[e[i].t2])continue;
                ++op[id[e[i].t1]][id[e[i].t1]], ++op[id[e[i].t2]][id[e[i].t2]];
                --op[id[e[i].t1]][id[e[i].t2]], --op[id[e[i].t2]][id[e[i].t1]];
            }
        }
        ans = ans * (int) (Guass(ID - 1) + 0.5) % MOD;
    }
    cout << (ans + MOD) % MOD;
    return 0;
}

行列式取模

        有些题目需要答案对某个数取模，这样我们需要一个可以取模的行列式板子。放到这里了：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
int det(int n)
{
    int i, j, k, res = 1, x, y, tmp, f = 1;
    for (i = 1; i <= n; ++i)
        for (j = 1; j <= n; ++j) B[i][j] = (B[i][j] + mod) % mod;
    for (i = 1; i <= n; ++i) {
        for (j = i + 1; j <= n; ++j) {
            x = B[i][i], y = B[j][i];
            while (y) {
                tmp = x / y;
                x %= y;
                swap(x, y);
                for (k = i; k <= n; ++k)
                    B[i][k] = ((B[i][k] - 1ll * tmp * B[j][k]) % mod + mod) % mod;
                for (k = i; k <= n; ++k) swap(B[i][k], B[j][k]);
                f = -f;
            }
        }
        if (!B[i][i]) return 0;
        res = 1ll * res * B[i][i] % mod;
    }
    if (f == -1) res = mod - res;
    return (res % mod + mod) % mod;
}