Matrix-Tree 定理
矩阵树定理，看名字就知道是把矩阵是树联系起来的一个定理，可以用于求图的任意/最大/最小生成树个数。
前置概念——行列式
他问我，一个 \(3\) 行 \(4\) 列的行列式怎么算，把我吓死了，你那叫懂吗？——汤神。
一般地，对于一个 \(n\times n\) 的方阵 \(\bf A\)，我们定义它的行列式为：
其中 \(\mathcal{P}\) 是一个 \(1\sim n\) 的排列，\(\mathrm{sgn}(\mathcal{P})\) 为 \(-1\) 当且仅当 \(\mathcal{P}\) 中逆序对个数为奇数，为 \(1\) 当且仅当 \(\mathcal{P}\) 中逆序对个数为偶数，\(\bf A\) 的行列式除 \(\mathrm{det}(\mathbf{A})\) 外，还记作 \(|\mathbf{A}|\)。注意到按照定义直接做复杂度是阶乘级别的，显然无法接受，考虑如何较快求出某个方阵的行列式。
首先要了解行列式的性质：
- \(|\mathbf{I}|=1\)。这个比较好想，因为除了 \(\mathcal{P}=(1,2,3\cdot\cdot\cdot,n)\) 以外的情况均为 \(0\)，而这种情况对应的值为 \(1\)。类似地，上三角矩阵和下三角矩阵行列式均为对角线乘积。
- 交换矩阵的两行，行列式变号。注意到交换行对于每个排列只会影响逆序对个数，而交换两个数对逆序对的影响必定是奇数（考虑这俩之间变一个，而其他的逆序对都是成对变的），所以最终结果就是整体变号。
- 行的线性性：
- 如果某一行乘以 \(t\)，则行列式也乘以 \(t\)。毕竟对于每个排列，每行选且仅选一个。
- 满足以下恒等式：\[\begin{vmatrix}a+x&b+y\\c&d\end{vmatrix}=\begin{vmatrix}a&b\\c&d\end{vmatrix}+\begin{vmatrix}x&y\\c&d\end{vmatrix} \]依然考虑每行选且仅选一个，乘法分配律就好。
- 有某两行一样的矩阵，其行列式为 \(0\)。考虑交换这两行，根据第二个性质行列式应该变号，但因为交换并没有实际意义，行列式应该不变，则值只能是 \(0\)。
- 用矩阵的一行加上另一行的倍数，行列式不变，即：\[\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\d&e&f\end{vmatrix}=\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\d+ta&e+tb&f+te\end{vmatrix} \]这个证明很巧妙，考虑这个矩阵：\[\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\ta&tb&tc\end{vmatrix} \]根据 3(1), 4 可以知道，它的行列式为 \(0\)，根据 3(2) 有：\[\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\d+ta&e+tb&f+te\end{vmatrix}=\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\d&e&f\end{vmatrix}+\begin{vmatrix}a&b&c\\\vdots&\vdots&\vdots\\ta&tb&tc\end{vmatrix} \]得证。
刚刚的性质中，我们提到了上三角矩阵，它的行列式是比较好求的。考虑高斯消元把原方阵消成一个上三角矩阵，然后就能求这个上三角矩阵的行列式了，而又注意到高斯消元的整个过程中，我们只会交换两行，给某行乘一个常数和用一行减去另一行，而这些操作对于行列式的影响是已知的。所以我们就可以通过高斯消元 \(\mathcal{O}(n^3)\) 求出方阵的行列式了。
int det()
{
int ans = 1, f = 1;
for (int j = 1; j <= n; ++j)
{
for (int i = j; i <= n; ++i)
{
if (!a[i][j]) continue;
if (i != j) std::swap(a[i], a[j]), f *= -1;
break;
}
if (!a[j][j]) return 0;
ans = 1ll * ans * a[j][j] % mod; int inv = ksm(a[j][j], mod - 2);
for (int k = j; k <= n; ++k) a[j][k] = 1ll * a[j][k] * inv % mod;
for (int i = j + 1; i <= n; ++i)
{
int t = a[i][j];
for (int k = j; k <= n; ++k)
a[i][k] = (a[i][k] - 1ll * t * a[j][k] % mod + mod) % mod;
}
}
return (ans * f + mod) % mod;
}
但现在有个问题，这里面会用到逆元。如果题目中要求取模的数不是质数怎么办？用 double 乱搞固然是一种解法，更稳妥的方法是放弃高斯消元，选择复杂度稍劣的辗转相除消元。能在不需要逆元的情况下 \(\mathcal{O}(n^3+n^2\log n)\) 的复杂度求出行列式。
int det()
{
int ans = 1, f = 1;
for (int i = 1; i <= n; ++i)
{
for (int j = i + 1; j <= n; ++j)
while (a[j][i])
{
int t = a[i][i] / a[j][i];
for (int k = i; k <= n && t; ++k)
a[i][k] = (a[i][k] - 1ll * t * a[j][k] % mod + mod) % mod;
std::swap(a[i], a[j]); f *= -1;
}
ans = 1ll * ans * a[i][i] % mod;
}
return (ans * f + mod) % mod;
}
正片——Matrix-Tree 定理
约定本文中提到的图均满足允许重边，不允许自环。
对于一张无向无权图，定义其邻接矩阵为 \(\bf A\)，度数矩阵为 \(\bf D\)，其中 \(\mathbf{D}_{i,i}=\mathrm{deg}_i,\mathbf{D}_{i,j}=0(i\ne j)\)。则其 \(\rm Kirchhoff\) 矩阵 \(\bf K\) 定义为：
令 \(\bf K_i\) 为 \(\bf K\) 矩阵去掉第 \(i\) 行第 \(i\) 列的结果，则原图的生成树个数为：
其中 \(i\) 任意。证明比较复杂，OI wiki 上的要用 LGV 引理，不会，所以就不证了。
考虑无向有权图，如果我们定义邻接矩阵 \(\mathbf{A}_{i,j}=w_{i,j}\)，度数矩阵 \(\mathbf{D}_{i,i}=\sum_{(i,v)\in E}w_{i,v}\)（即所有与该点相邻的边权值和），套用上面的思路就可以得到所有生成树边权积之和。考虑求出这个的实际意义，比较好想到的就是带重边的情况，把重边个数理解为边权，就能求出生成树个数了。
扩展——Matrix-Tree 定理的更多形态
刚刚说了无向图，事实上对于有向图，Matrix-Tree 定理也是能派上用场的。考虑定义入度矩阵 \(\mathbf{D^{in}}\)，其中 \(\mathbf{D^{in}}_{i,i}=\mathrm{in}_i,\mathbf{D^{in}}_{i,j}=0(i\ne j)\)，定义出度矩阵 \(\mathbf{D^{out}}\)，其中 \(\mathbf{D^{out}}_{i,i}=\mathrm{out}_i,\mathbf{D^{out}}_{i,j}=0(i\ne j)\)。它们两个分别与邻接矩阵 \(\bf A\) 作差可以得到两种 \(\rm Kirchhoff\) 矩阵：
- \(\mathbf{K^{leaf}}=\mathbf{D^{in}}-\mathbf{A}\)
- \(\mathbf{K^{root}}=\mathbf{D^{out}}-\mathbf{A}\)
则有结论：
- 有向图以 \(i\) 为根的外向生成树（所有边均指向叶子的方向）个数为：\[\mathrm{det}(\mathbf{K^{leaf}_i}) \]其中 \(\mathbf{K^{leaf}_i}\) 表示去掉第 \(i\) 行第 \(i\) 列，注意在无向图中这个 \(i\) 随便选是因为无向图无所谓谁是根。
- 有向图以 \(i\) 为根的内向生成树（所有边均指向根的方向）个数为：\[\mathrm{det}(\mathbf{K^{root}_i}) \]其中 \(\mathbf{K^{root}_i}\) 定义类似。
还有一个比较冷门的定理，还没做到相关的题，所以就扔公式跑路了：
设 \(G=(V,E)\) 是有向欧拉图，则 \(G\) 的不同欧拉回路的个数 \(\mathrm{ec}(G)\) 为：
\[\mathrm{ec}(G)=\mathbf{K^{root}_i}\prod_{v\in V}(\mathrm{deg}_v-1)! \]注意在欧拉图中，对于任意的两个结点 \(i,j\)，均有 \(\mathbf{K^{root}_i}=\mathbf{K^{root}_j}\)。
板子题 代码
#include <cstdio>
#include <algorithm>
const int N = 310, mod = 1e9 + 7; int *a[N], _a[N][N], n;
int ksm(int a, int b)
{
int ret = 1;
while (b)
{
if (b & 1) ret = 1ll * ret * a % mod;
a = 1ll * a * a % mod; b >>= 1;
}
return ret;
}
int det()
{
int ans = 1, f = 1;
for (int j = 2; j <= n; ++j)
{
for (int i = j; i <= n; ++i)
{
if (!a[i][j]) continue;
if (i != j) std::swap(a[i], a[j]), f *= -1;
break;
}
if (!a[j][j]) return 0;
ans = 1ll * ans * a[j][j] % mod; int inv = ksm(a[j][j], mod - 2);
for (int k = j; k <= n; ++k) a[j][k] = 1ll * a[j][k] * inv % mod;
for (int i = j + 1; i <= n; ++i)
{
int t = a[i][j];
for (int k = j; k <= n; ++k)
a[i][k] = (a[i][k] - 1ll * t * a[j][k] % mod + mod) % mod;
}
}
return (ans * f + mod) % mod;
}
int main()
{
int m, t; scanf("%d%d%d", &n, &m, &t);
for (int i = 1; i <= n; ++i) a[i] = _a[i];
for (int i = 1, x, y, z; i <= m; ++i)
{
scanf("%d%d%d", &x, &y, &z);
if (!t) (a[x][y] += (mod - z)) %= mod, (a[y][x] += (mod - z)) %= mod, (a[x][x] += z) %= mod, (a[y][y] += z) %= mod;
else (a[x][y] += (mod - z)) %= mod ,(a[y][y] += z) %= mod;
}
printf("%d\n", det()); return 0;
}
应用1——Matrix-Tree 定理求最大/最小生成树计数
对于无向有权图，我们想求出它的最小生成树个数。注意到在不同的最小生成树中，相同边权的边个数是相同的，且去掉这些相同边权的边，图的连通性是相同的。证明可以感性通过 \(\rm kruskal\) 的过程理解。
所以我们就仍然可以延续 \(\rm kruskal\) 的思路，把边排序，从小到大依次考虑每一种边权。对于每一种边权，我们先把端点已经连通的边去掉，把剩下的加入原图和一个新图中，注意原图要用并查集判断连通性保证一直是森林。求出这个新图中的生成树个数就是在所有最小生成树方案中这个边权对应的方案，求出所有边权对应的方案乘起来即得最终答案。有一些小细节放代码里面了。
板子题 代码：
#include <cstdio>
#include <vector>
#include <cstring>
#include <algorithm>
const int N = 110, M = 1e3 + 10, mod = 31011;
int w[M]; std::vector<std::pair<int, int> > v[M];
struct edge{ int u, v, w; }E[M]; int *a[N], _a[N][N], id[N], f[N], cnt;
int getf(int x) { return x == f[x] ? x : f[x] = getf(f[x]); }
int det(int n)
{
int ans = 1, f = 1;
for (int i = 1; i <= n; ++i)
{
for (int j = i + 1; j <= n; ++j)
while (a[j][i])
{
int t = a[i][i] / a[j][i];
for (int k = i; k <= n; ++k)
a[i][k] = (a[i][k] - 1ll * t * a[j][k] % mod + mod) % mod;
std::swap(a[i], a[j]); f *= -1;
}
ans = 1ll * ans * a[i][i] % mod;
}
return (ans * f + mod) % mod;
}
int main()
{
int n, m, ans = 1; scanf("%d%d", &n, &m);
for (int i = 1; i <= m; ++i) scanf("%d%d%d", &E[i].u, &E[i].v, &E[i].w), w[i] = E[i].w;
std::sort(w + 1, w + m + 1); int l = std::unique(w + 1, w + m + 1) - w - 1;
for (int i = 1; i <= m; ++i)
{
int id = std::lower_bound(w + 1, w + l + 1, E[i].w) - w;
v[id].emplace_back(E[i].u, E[i].v);
}
for (int i = 1; i <= n; ++i) f[i] = i;
for (int t = 1; t <= l; ++t)
{
memset(id, 0, sizeof (id)); cnt = 0; //注意新图要单独标号，不然如果带上未加入新图的点就永远不可能形成生成树了
for (auto e : v[t])
{
int u = getf(e.first), v = getf(e.second);
if (u == v) continue;
if (!id[u]) id[u] = ++cnt; if (!id[v]) id[v] = ++cnt;
}
for (int i = 1; i <= n; ++i) for (int j = 1; j <= n; ++j) _a[i][j] = 0;
for (int i = 1; i <= n; ++i) a[i] = _a[i]; //新图每次要清空
for (auto e : v[t])
{
int u = getf(e.first), v = getf(e.second);
if (u == v) continue; u = id[u]; v = id[v];
(a[u][v] += (mod - 1)) %= mod; (a[v][u] += (mod - 1)) %= mod;
(a[v][v] += 1) %= mod; (a[u][u] += 1) %= mod;
}
for (auto e : v[t])
{
int u = getf(e.first), v = getf(e.second);
if (u == v) continue; f[u] = v;
}
//这里其实变相的增加了一个新点 cnt+1 并把它与所有点连边，但因为求生成树要删掉一行一列，所以就把 cnt+1 行列删掉了。
for (int i = 1; i <= n; ++i) if (id[i] && f[i] == i) ++a[id[i]][id[i]];
ans = 1ll * ans * det(cnt) % mod;
}
//最后判一下连通性看看原图到底有没有生成树
for (int i = 1; i <= n; ++i) if (getf(i) != getf(1)) { ans = 0; break; }
printf("%d\n", ans); return 0;
}
最大生成树只需要更改一下循环顺序即可。
应用2——边权的 \(k\) 次方和
对于无向有权图，我们能求出的是：
其中 \(\mathcal{T}\) 为生成树。但如何求这个式子呢：
看起来一脸不可求，但事实上，基本所有跟生成树计数的问题都要往 Matrix-Tree 定理上面靠，而靠的方法就是把式子转化成我们能求出的那个样子。考虑用 EGF 暴力拆掉这个 \(k\) 次方：
注意到如果把 \(k!\) 提出来就是上面能求的形式了，只不过我们要把每个边的边权设置为：
因为最后只要 \(z^k\) 前面的系数，所以运算的时候可以只保留到第 \(k\) 次。确定好边权后用 Matrix-Tree 求最原始的式子，然后乘上 \(k!\) 即可。用朴素的多项式运算可以做到 \(\mathcal{O}(n^3k^2)\)，当然用 NTT 那一套可以优化到 \(\mathcal{O}(n^3k\log k)\)。
板子题 代码：
#include <cstdio>
#include <cstring>
#include <algorithm>
const int N = 40, mod = 998244353; typedef long long ll;
int fac[N], ifac[N], w[N][N], n, k;
inline int ksm(int a, int b)
{
int ret = 1;
while (b)
{
if (b & 1) ret = (ll)ret * a % mod;
a = (ll)a * a % mod; b >>= 1;
}
return ret;
}
struct poly
{
int x[N];
poly getinv()
{
poly a, ret; memset(ret.x, 0, sizeof (ret.x)); ret.x[0] = ksm(x[0], mod - 2);
for (int i = 1; i <= k; ++i) a.x[i] = (ll)x[i] * ret.x[0] % mod;
for (int i = 1; i <= k; ++i) for (int j = 1; j <= i; ++j)
ret.x[i] = (ret.x[i] + (ll)(mod - a.x[j]) * ret.x[i - j]) % mod;
return ret;
}
}; poly *a[N], _a[N][N], O, t, res;
poly operator+(poly a, poly b)
{
poly ret;
for (int i = 0; i <= k; ++i) ret.x[i] = (a.x[i] + b.x[i]) % mod;
return ret;
}
poly operator-(poly a, poly b)
{
poly ret;
for (int i = 0; i <= k; ++i) ret.x[i] = (a.x[i] - b.x[i] + mod) % mod;
return ret;
}
poly operator*(poly a, poly b)
{
poly ret; memset(ret.x, 0, sizeof (ret.x));
for (int i = 0; i <= k; ++i)
for (int j = 0; j <= k - i; ++j) (ret.x[i + j] += (ll)a.x[i] * b.x[j] % mod) %= mod;
return ret;
}
int det()
{
res.x[0] = 1; int f = 1;
for (int j = 2; j <= n; ++j)
{
for (int i = j; i <= n; ++i)
{
if (!a[i][j].x[0]) continue;
if (i != j) std::swap(a[i], a[j]), f *= -1;
break;
}
res = res * a[j][j]; t = a[j][j].getinv();
for (int k = j; k <= n; ++k) a[j][k] = a[j][k] * t;
for (int i = j + 1; i <= n; ++i)
{
t = a[i][j];
for (int k = j; k <= n; ++k) a[i][k] = a[i][k] - t * a[j][k];
}
}
if (f == -1) res = O - res;
return (ll)res.x[k] * fac[k] % mod;
}
int main()
{
scanf("%d%d", &n, &k); fac[0] = ifac[0] = 1;
for (int i = 1; i <= n; ++i) a[i] = _a[i];
for (int i = 1; i <= k; ++i) fac[i] = (ll)fac[i - 1] * i % mod;
ifac[k] = ksm(fac[k], mod - 2);
for (int i = k - 1; i >= 1; --i) ifac[i] = (ll)ifac[i + 1] * (i + 1) % mod;
for (int i = 1; i <= n; ++i) for (int j = 1; j <= n; ++j) scanf("%d", &w[i][j]);
for (int i = 1; i <= n; ++i)
for (int j = i + 1; j <= n; ++j)
{
int mul = 1;
for (int t = 0; t <= k; ++t)
{
int v = (ll)mul * ifac[t] % mod;
a[i][j].x[t] = a[j][i].x[t] = (mod - v) % mod;
(a[i][i].x[t] += v) %= mod; (a[j][j].x[t] += v) %= mod;
mul = (ll)mul * w[i][j] % mod;
}
}
printf("%d\n", det()); return 0;
}
