模板：矩阵树定理

wind__whisper

于 2022-03-17 15:24:06 发布

阅读量223

点赞数

分类专栏：模板数学图论文章标签：矩阵 c++ 图论

本文链接：https://blog.csdn.net/BUG_Creater_jie/article/details/123551294

版权

模板同时被 3 个专栏收录

90 篇文章 2 订阅

订阅专栏

数学

60 篇文章 0 订阅

订阅专栏

图论

57 篇文章 0 订阅

订阅专栏

文章目录

前言
解析

所谓矩阵树定理，就是用矩阵解决树问题的定理。

（逃）

前言

神奇科技。
之前一直没有写博客，觉得还是写一发比较好。
~~证明什么的是不可能会的~~
背下来背下来！

解析

无向图

定义一个矩阵 $A$ ：
$A_{i,i}=d_i\\A_{i,j}=-\#(i,j)(i\ne j)$
其中 $\#(i,j)$ 表示这条边的数量。
求出这个矩阵的任意一个余子式即可。

有向图

根向树

定义一个矩阵 $A$ ：
$A_{i,i}=d_i^+\\A_{i,j}=-\#(i,j)(i\ne j)$
去掉第 $i$ 行得到的余子式就是以 $i$ 为根的答案。

叶向树

定义一个矩阵 $A$ ：
$A_{i,i}=d_i^-\\A_{i,j}=-\#(i,j)(i\ne j)$
唯一的区别就是主对角线从出度变成了入度（~~为什么呢？因为根的思想太陈旧，已经out了~~）

code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
#define ok debug("OK\n")
using namespace std;

const int N=305;
const int mod=1e9+7;

inline ll read(){
	ll x(0),f(1);char c=getchar();
	while(!isdigit(c)) {if(c=='-')f=-1;c=getchar();}
	while(isdigit(c)) {x=(x<<1)+(x<<3)+c-'0';c=getchar();}
	return x*f;
}

inline ll ksm(ll x,ll k){
	ll res(1);
	while(k){
		if(k&1) res=res*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return res;
}

int n,m;
ll a[N][N];
ll calc(int n){
	ll ans(1);
	for(int i=1;i<=n;i++){
		for(int j=i;j<=n;j++){
			if(a[j][i]){
				if(j!=i) swap(a[i],a[j]);
				break;
			}
		}
		for(int j=i+1;j<=n;j++){
			ll d=a[j][i]*ksm(a[i][i],mod-2)%mod;
			for(int k=i;k<=n;k++) a[j][k]=(a[j][k]+mod-a[i][k]*d%mod)%mod;
		}
	}
	for(int i=1;i<=n;i++) ans=ans*a[i][i]%mod;
	return ans;
}
void work0(){
	for(int i=1;i<=m;i++){
		int x=read(),y=read(),w=read();
		if(x==y) continue;
		(a[x][x]+=w)%=mod;
		(a[y][y]+=w)%=mod;
		(a[x][y]+=mod-w)%=mod;
		(a[y][x]+=mod-w)%=mod;
	}
	printf("%lld\n",calc(n-1));
}
void work1(){
	for(int i=1;i<=m;i++){
		int x=read(),y=read(),w=read();
		if(x==y) continue;
		(a[y][y]+=w)%=mod;
		(a[x][y]+=mod-w)%=mod;
	}
	for(int i=1;i<n;i++){
		for(int j=1;j<n;j++) a[i][j]=a[i+1][j+1];
	}
	printf("%lld\n",calc(n-1));
}

signed main(){
#ifndef ONLINE_JUDGE
//	freopen("a.in","r",stdin);
//	freopen("a.out","w",stdout);
#endif
	n=read();m=read();int op=read();
	if(op==0) work0();
	else work1();
	return 0;
}
/*
3
2 1 1
1 1
*/

带权图

如果一棵生成树的权值定义为边权乘积，直接把边理解成边权条1权边在对应的位置加边权即可。

如果一棵生成树的权值定义为边权加和，需要在矩阵内维护一个一次函数，最终在 $mod x^2$ 意义下求出的行列式的一次项就是答案。
因为这就相当于单独考虑一条边的边权，看它能加入多少棵生成树中。

code

#include<bits/stdc++.h>
using namespace std;
#define ll long long
#define ull unsigned long long
#define debug(...) fprintf(stderr,__VA_ARGS__)
#define ok debug("OK\n")
using namespace std;

const int N=35;
const int M=1050;
const int S=2e5+100;
const int inf=1e9;
const int mod=998244353;

inline ll read(){
	ll x(0),f(1);char c=getchar();
	while(!isdigit(c)) {if(c=='-')f=-1;c=getchar();}
	while(isdigit(c)) {x=(x<<1)+(x<<3)+c-'0';c=getchar();}
	return x*f;
}

inline ll ksm(ll x,ll k){
	ll res(1);
	while(k){
		if(k&1) res=res*x%mod;
		x=x*x%mod;
		k>>=1;
	}
	return res;
}

int n,m;
int mx;
int prime[S],vis[S],tot;
ll mu[S];
void init(int n){
	mu[1]=1;
	for(int i=2;i<=n;i++){
		if(!vis[i]){
			prime[++tot]=i;
			mu[i]=-1;
		}
		for(int j=1;j<=tot&&prime[j]<=n/i;j++){
			int now=prime[j];
			vis[i*now]=1;
			if(i%now==0){
				mu[i*now]=0;break;
			}
			mu[i*now]=-mu[i];
		}
	}
	for(int i=1;i<=n;i++) mu[i]=(mu[i]+mod)%mod;
	return;
}

struct node{
	ll a,b;
};
node operator + (const node x,const node y){return (node){(x.a+y.a)%mod,(x.b+y.b)%mod};}
node operator - (const node x,const node y){return (node){(x.a+mod-y.a)%mod,(x.b+mod-y.b)%mod};}
node operator * (const node x,const node y){return (node){(x.a*y.b+x.b*y.a)%mod,x.b*y.b%mod};}
node operator / (const node x,const node y){
	ll niv=ksm(y.b,mod-2);
	return (node){(x.a*y.b-x.b*y.a%mod+mod)%mod*niv%mod*niv%mod,x.b*niv%mod};
}
void operator += (node &x,const node y){x=x+y;}
void operator -= (node &x,const node y){x=x-y;}
void operator *= (node &x,const node y){x=x*y;}
void operator /= (node &x,const node y){x=x/y;}

int u[M],v[M],w[M];
node a[N][N];
ll g[S],f[S];
void print(int n){
	puts("\n------\n");
	for(int i=1;i<=n;i++){
		for(int j=1;j<=n;j++) printf("(%lld %lld)  ",a[i][j].a,a[i][j].b);
		puts("");
	}
}
ll calc(int n){
	node ans=(node){0,1};
	for(int i=1;i<=n;i++){
		for(int j=i+1;j<=n;j++){
			node d=a[j][i]/a[i][i];
			for(int k=i;k<=n;k++) a[j][k]-=(a[i][k]*d);
		}
	}
	for(int i=1;i<=n;i++){
		//printf("  i=%d a=%lld ans=%lld\n",i,a[i][i],ans[i])
		ans=ans*a[i][i];
	}
	return ans.a;
}
void work(int x){
	memset(a,0,sizeof(a));
	int cnt(0);
	for(int i=1;i<=m;i++){
		if(w[i]%x) continue;
		a[u[i]][u[i]]+=(node){w[i],1};
		a[v[i]][v[i]]+=(node){w[i],1};
		a[u[i]][v[i]]-=(node){w[i],1};
		a[v[i]][u[i]]-=(node){w[i],1};
		cnt++;
	}
	if(cnt<n-1) return;
	//printf("\n---x=%d\n",x);
	//print(n);
	g[x]=calc(n-1);
	//printf("ans=%lld\n",g[x]);
	return;
}

signed main() {
#ifndef ONLINE_JUDGE
//	freopen("a.in","r",stdin);
//	freopen("a.out","w",stdout);
#endif
	n=read();m=read();
	for(int i=1;i<=m;i++){
		u[i]=read();v[i]=read();w[i]=read();
		mx=max(mx,w[i]);
	}
	init(mx);
	for(int x=1;x<=mx;x++) work(x);
	ll ans(0);
	for(int i=1;i<=mx;i++){
		for(int j=1;j*i<=mx;j++) (f[i]+=g[i*j]*mu[j])%=mod;
		(ans+=i*f[i])%=mod;
	}
	printf("%lld\n",ans);
	return 0;
}
/*
3
2 1 1
1 1
*/