【数学】Matrix-Tree 定理
题目描述
给定一张 \(n\) 个结点 \(m\) 条边的带权图(可能为无向图,可能为有向图)。
定义其一个生成树 \(T\) 的权值为 \(T\) 中所有边权的乘积。
求其所有不同生成树的权值之和,对 \(10^9+7\) 取模。
注意:
-
本题中,有向图的生成树指的是 以 \(1\) 为根的外向树;
-
两棵生成树 \(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\) ,生成树数量是所有边的重复次数乘积,我们猜想,如果一条边的权值不是重复次数,而是其他东西,还可以通过直接改变度数来算吗?
事实上是可以的,也就是说,矩阵树其实算的是这样一个东西:
其中 \(w_i\) 是边 \(i\) 的权值,直接将 \(D,A\) 改为相应权值就好了。
例题:
P3317 [SDOI2014] 重建 - 洛谷 | 计算机科学教育新生态 (luogu.com.cn)
这题要求
就是
将边权设为 \(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;
}