【数学】Matrix-Tree 定理

题目描述

给定一张 $n$ 个结点 $m$ 条边的带权图（可能为无向图，可能为有向图）。

定义其一个生成树 $T$ 的权值为 $T$ 中所有边权的乘积。

求其所有不同生成树的权值之和，对 ${10}^9+7$ 取模。

---

注意：

1. 本题中，有向图的生成树指的是 以 $1$ 为根的外向树；

2. 两棵生成树 $T_1,T_2$ 不同，当且仅当存在存在一条边 $e$，满足 $e\in T_1,e\notin T_2$。

$1\le n\le 300$ 。

算法描述（暂缺证明）

考虑我们定义这张无向图的邻接矩阵 $A$ ，$A_{i,j}=1$ 当且仅当 $i,j$ 有边。度数矩阵 $D$ ，$D_{i,i}$ 等于 $de{g}_i$ ，其他值等于 $0$ 。

这张图的 基尔霍夫矩阵 为 $K=D-A$ 。

那么这个基尔霍夫矩阵的一个 $n-1$ 阶主子式（去掉任意一行，一列）的行列式就是这张图的生成树个数。

可以用于重边情况。

至于行列式，直接 $\mathrm{\Theta }(n^3)$ 高斯消元消成上三角后主对角线乘起来即可。

#include<bits/stdc++.h>
using namespace std;
const int N = 305,MOD = 1e9 + 7;
typedef long long ll;
ll D[N][N],A[N][N],K[N][N];
int n,m,t;
inline ll ksm(ll base,int pts)
{
	ll ret = 1;
	for(;pts > 0;pts >>= 1,base = base * base % MOD)
		if(pts & 1)
			ret = ret * base % MOD;
	return ret;
}
inline ll Gauss()
{
	ll res = 1;
	for(int i = 1;i <= n;i++)
	{
		int tmp = i;
		for(;tmp <= n;tmp++) if(K[tmp][i]) break;
		if(tmp > n) return 0;
		swap(K[tmp],K[i]);  
		if(tmp != i) res = (MOD - res);
		ll inv = ksm(K[i][i],MOD - 2);
		res = res * K[i][i] % MOD;
		for(int j = i + 1;j <= n;j++)
		{
			ll rate = K[j][i] * inv % MOD;
			for(int k = 1;k <= n;k++) K[j][k] = (K[j][k] - (rate * K[i][k] % MOD) + MOD) % MOD;
		}
	}
	return res;
}
int main()
{
	cin>>n>>m>>t;
	for(int i = 1,x,y,z;i <= m;i++)
	{
		cin>>x>>y>>z;
		if(x == y) continue;
		if(t == 0)
		{
			D[x][x] = (D[x][x] + z) % MOD; D[y][y] = (D[y][y] + z) % MOD;
			A[x][y] = (A[x][y] + z) % MOD; A[y][x] = (A[y][x] + z) % MOD;
		}
		else
		{
			D[y][y] = (D[y][y] + z) % MOD; A[x][y] = (A[x][y] + z) % MOD;
		}
	}
	for(int i = 1;i <= n;i++) for(int j = 1;j <= n;j++) K[i][j] = (D[i][j] - A[i][j] + MOD) % MOD;
	n--;
	for(int i = 1;i <= n;i++)
		for(int j = 1;j <= n;j++)
			K[i][j] = K[i + 1][j + 1]; 
	ll ans = Gauss();
	cout<<ans;
	return 0;
 } 

拓展应用

1.有向图：

对于有向图，如果求外向树，就将度数矩阵改为入度，反之，内向树改为出度即可。

还有一个结论就是以 $x$ 为根，就删去 $K$ 的第 $x$ 行和第 $x$ 列。

2.意义：

我们发现，重边可以通过直接在 $D,A$ 中加度数来解决，我们又知道，对于一个树 $T$ ，生成树数量是所有边的重复次数乘积，我们猜想，如果一条边的权值不是重复次数，而是其他东西，还可以通过直接改变度数来算吗？

事实上是可以的，也就是说，矩阵树其实算的是这样一个东西：

$$\sum _T\prod _{i\in e}{w}_i$$

其中 $w_i$ 是边 $i$ 的权值，直接将 $D,A$ 改为相应权值就好了。

例题：

P3317 [SDOI2014] 重建 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)

这题要求

$$\sum _T(\prod _{u\in e}{w}_u\prod _{u\in other}(1-w_u))$$

就是

$$\prod _u(1-w_u)\sum _T\prod _{u\in e}{w}_u\times \frac{1}{1-w_u}$$

将边权设为 $w_u\times \frac{1}{1-w_u}$ 即可。

#include<bits/stdc++.h>
using namespace std;
const int N = 55;
const double eps = 1e-7;
struct Matrix{
	double a[N][N];
}D,A,K,P;
int n;
inline double getdet(Matrix x)
{
	double ret = 1;
	for(int i = 1;i <= n - 1;i++)
	{
		int tmp = i;
		while(tmp <= n - 1 && fabs(x.a[tmp][i]) < eps) tmp++;
		if(tmp == n) return 0;
		if(tmp ^ i)
		{
			for(int j = 1;j <= n - 1;j++) swap(x.a[tmp][j],x.a[i][j]);
			ret = -ret;
		}
		ret *= x.a[i][i];
		for(int j = i + 1;j <= n - 1;j++)
		{
			double rate = x.a[j][i] / x.a[i][i];
			for(int k = 1;k <= n - 1;k++) x.a[j][k] = x.a[j][k] - x.a[i][k] * rate;
		}
	}	
	return ret;
}
int main()
{
	cin>>n;
	for(int i = 1;i <= n;i++)
		for(int j = 1;j <= n;j++)
			cin>>A.a[i][j];
	for(int i = 1;i <= n;i++)
		for(int j = 1;j <= n;j++)
			K.a[i][j] = D.a[i][j] = 0;
	double prod = 1;
	for(int i = 1;i <= n;i++)
		for(int j = i + 1;j <= n;j++)
			prod *= 1 + eps - A.a[i][j];
	for(int i = 1;i <= n;i++)
		for(int j = 1;j <= n;j++)
			A.a[i][j] = A.a[i][j] / (1 + eps - A.a[i][j]);
	for(int i = 1;i <= n;i++) 
		for(int j = i + 1;j <= n;j++)
			D.a[i][i] += A.a[i][j],D.a[j][j] += A.a[i][j];
	for(int i = 1;i <= n - 1;i++)
		for(int j = 1;j <= n - 1;j++)
			K.a[i][j] = D.a[i][j] - A.a[i][j];
	cout<<fixed<<setprecision(15)<<prod * getdet(K);
	return 0;
}