526互联

【数学】Matrix-Tree 定理

发布时间 2023-11-21 22:02:14作者: The_Last_Candy

题目描述

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

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

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

注意：

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

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

\(1 \leq n \leq 300\) 。

算法描述（暂缺证明）

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

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

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

可以用于重边情况。

至于行列式，直接 \(\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;
}