《行列式+矩阵树定理》 随便学学

今天学了一下矩阵树定理感觉很 \(nb\) 。

参考了 \(cmd\) 的文章（明明就是大量抄袭）

前置知识：行列式

什么是行列式：

对于一个矩阵 \(A[1\dots n][1\dots n]\)

他的行列式为 \(det(A)=\sum\limits_P (-1)^{\mu(P)}\prod\limits_{i=1}^n A[i][p_i]\) 。（行列式是一个具体的数值）

其中 \(P\) 是 \(1\dots n\) 的一个排列， \(\mu (P)\) 是排列 \(P\) 的逆序对个数。

其实可以理解成每一行选一个，且在列上不重复的乘积的和，然后多了一个逆序对相关的东西做贡献。

直接按照定义计算时间复杂度是 \(O(n!\times n)\) 的，难以运用到实际中，考虑他的性质。

行列式的性质：

  • (1)单位矩阵的行列式为 \(1\) 。

这个比较显然，因为只有对角线这一种值可以获得。

  • (2)交换两行矩阵，行列式变号。

可以考虑不交换矩阵，想象成交换 \(p_i,p_j\) ，而交换两个数，造成逆序对差异是奇数，所以系数乘 \(-1\) ，所以变号。

  • (3)某一行乘 \(t\) ，那么行列式乘 \(t\)

这个显然，因为每一行只选一个。

  • (4)有一个公式：

\[\begin{vmatrix} a+a'&b+b' \\ c&d \end{vmatrix} = \begin{vmatrix} a+&b+ \\ c&d \end{vmatrix} + \begin{vmatrix} a'&b' \\ c&d \end{vmatrix} \]

这个也是显然的，因为容易知道满足分配率。

  • 有两个行相同，矩阵行列式为 \(0\)

因为如果交换了，根据 \((2)\) ，要变号，但是因为两行相同，所以交换了实际上行列式不变，所以只能行列式等于 \(0\) 。

  • 矩阵一行加另一行的倍数行列式不变。

容易由 \(3,4\) 得到。

\[\begin{vmatrix} a&b \\ c+ak&d+bk \end{vmatrix} = \begin{vmatrix} a&b \\ c&d \end{vmatrix} + \begin{vmatrix} a&b \\ ak&bk \end{vmatrix} \]

又因为

\[\begin{vmatrix} a&b \\ ak&bk \end{vmatrix} = k \begin{vmatrix} a&b \\ a&b \end{vmatrix} =k\times 0=0 \]

得证。

高消求解行列式：

因为高消要用到的是 \([\) 交换，一行加另一行的倍数，某一行乘一个常数 \(]\) 三个操作。而这三个操作的影响我们都在上面整出来了，就直接做高消，然后根据影响逆回去即可。

点击查看代码 
#include<bits/stdc++.h>
typedef long long LL;

using namespace std;
const int MAXN=610;
int n,P;
int a[MAXN][MAXN];
int gauss() {
	LL res=1,w=1;
	for(int i=1;i<=n;++i) {
		for(int j=i+1;j<=n;++j) {
			while(a[i][i]) {
				int div=a[j][i]/a[i][i];
				for(int k=i;k<=n;++k) {
					a[j][k]=(a[j][k]-1ll*div*a[i][k]%P+P)%P;
				}
				swap(a[i],a[j]); w=-w;
			}
			swap(a[i],a[j]); w=-w;
		}
	}
	for(int i=1;i<=n;++i) res=1ll*a[i][i]*res%P;
	res=1ll*w*res;
	return (res+P)%P;
}
int main () {
	scanf("%d%d",&n,&P);
	for(int i=1;i<=n;++i) {
		for(int j=1;j<=n;++j) {
			scanf("%d",&a[i][j]);
		}
	}
	printf("%d\n",gauss());
	return 0;
}

\(Part\ 1\) ：引入

这是一个能够在 \(O(n^3)\) 时间复杂度内求出一个图的生成树数量的算法，虽然我也不能证明，但就是可以。

具体的，实际上求的是这样一个东西 \(\sum\limits_{T} \prod\limits_{e\in E_T}e\)

就是求所有生成树的边的乘积。

一般情况：

（先只考虑无向图的情况）我们记 \(D\) 为度数矩阵，\(F\) 为邻接矩阵，记 \(K=D-F\)

那么 \(K\) 随便扣掉一行一列之后，他的行列式就是上面的这个东西。（所有边权记为 \(1\) ，那么实际上求的就是生成树个数，而如果有重边，直接加就好了）

上面的 \(K\) 有一个性质，不管随便扣掉哪行哪列，剩下的矩阵的行列式都一样。

例题：P4336 [SHOI2016] 黑暗前的幻想乡

容斥加矩阵树定理板子。

点击查看代码 
#include<bits/stdc++.h>
#define mp(x,y) make_pair(x,y)
#define mod 1000000007
using namespace std;
int TMP3301;
inline int read(){
    int res=0;
    char c;
    bool zf=0;
    while(((c=getchar())<'0'||c>'9')&&c!= '-');
    if(c=='-')zf=1;
    else res=c-'0';
    while((c=getchar())>='0'&&c<='9')res=(res<<3)+(res<<1)+c-'0';
    if(zf)return -res;
    return res;
}
int n;
int dta[21][21];
typedef pair<int,int>P;
vector<P>q[21];
inline void _swap(int x,int y){
	for(register int i=1;i<=n;++i){
		#define swap(x,y) TMP3301=x,x=y,y=TMP3301;
		swap(dta[i][x],dta[i][y]);
		#undef swap
	}
	return;
}
inline int det(){
	int ans=1;
	for(register int i=1;i<=n;++i){
		int p=-1;
		for(register int j=i;j<=n;++j)
			if(dta[i][j]){
				p=j;
				break;
			}
		if(!~p){
			return 0;
		}
		if(p!=i){
			_swap(i,p);
			ans*=-1;
		}
		for(register int j=i+1;j<=n;++j){
			while(dta[j][i]){
				int tmp=dta[j][i]/dta[i][i];
				for(register int k=i;k<=n;++k){
					dta[j][k]-=1ll*tmp*dta[i][k]%mod;
					dta[j][k]=(dta[j][k]%mod+mod)%mod;
				} 
				if(!dta[j][i]){
					break;
				}
				swap(dta[i],dta[j]);
				ans*=-1;
			}
		}
		ans=1ll*ans*dta[i][i]%mod;
	}
	return ans;
}
signed main(){
	n=read()-1;
	for(register int i=1;i<=n;++i){
		int cnt=read();
		while(cnt--){
			int x=read(),y=read();
			q[i].push_back(mp(x,y));
		}
	}
	int ans=0;
	for(register int i=1;i<(1<<n);++i){
		memset(dta,0,sizeof(dta));
		int cnt=0;
		for(register int j=1;j<=n;++j){
			if(!(i&(1<<(j-1)))){
				continue;
			}
			++cnt;
			for(register int k=0;k<q[j].size();++k){
				int x=q[j][k].first,y=q[j][k].second;
				++dta[x][y],++dta[y][x],--dta[x][x],--dta[y][y];
			}
		}
		ans=(ans+((cnt&1)?-1:1)*det())%mod;
	}
	printf("%d\n",(ans%mod+mod)%mod);
	return 0;
}

加权拓展：

加权或者有重边其实也可以直接处理，容易理解，根据乘法原理。

例题：P3317 [SDOI2014] 重建

化式子。

\(\sum\limits_{T} \prod\limits_{e\in E_T} p_e\prod\limits_{e\not \in E_T}(1-p_e)=\sum\limits_{T} \prod\limits_{e\in E_T} p_e\frac {\prod\limits_{e}(1-p_e)}{\prod \limits_{e\in E[T]}(1-p_e)}=\sum\limits_T\prod\limits_{e}(1-p_e)\prod\limits_{e\in T} \frac {p_e}{(1-p_e)}\)

然后对于 \(p_e=1\) 的，有 \(\frac 1{1-p_e}=\infty,\infty\approx \frac 1{eps}\)

点击查看代码 
#include<bits/stdc++.h>
#define eps 1e-9
typedef double DB;
typedef long long LL;

using namespace std;
const int N=60;
int n;
DB a[N][N];
DB gauss() {
	DB res=1;
	for(int i=1;i<n;++i) {
		int p=0;
		for(int j=i;j<n;++j) {
			if(fabs(a[i][i])>-eps) {
				p=j;
				break;
			}
		}
		if(!p) continue;
		swap(a[i],a[p]);
		for(int j=i+1;j<n;++j) {
			DB div=a[j][i]/a[i][i];
			for(int k=i;k<n;++k) {
				a[j][k]=a[j][k]-1ll*div*a[i][k];
			}
		}
	}
	for(int i=1;i<n;++i) {
		res=a[i][i]*res;
	}
	return res;
}
int main () {
	scanf("%d",&n);
	DB ans=1;
	for(int i=1;i<=n;++i) {
		for(int j=1;j<=n;++j) {
			DB x;
			scanf("%lf",&x);
			if(i==j) continue;
			if(x>1-eps) x-=eps;
			if(i<j) ans=ans*(1-x);
			a[i][j]=x/(1-x);
		}
	}
	for(int i=1;i<=n;++i) {
		for(int j=1;j<=n;++j) {
			if(i==j) continue;
			a[i][i]+=a[i][j];
			a[i][j]=-a[i][j];
		}
	}
	printf("%.4lf",gauss()*ans);
	return 0;
}

有向拓展：

实际上只要改一下度数矩阵就好了。

对于外向树而言。\(D_{i,i}=\sum\limits_{j=1} F_{j,i}\) （即为到这个点的边权和）

而对于内向树而言 \(D_{i,j}=\sum\limits_{j=1} F_{i,j}\) （即为从这个点出发的边权和）

而对于有向边，自然也就有根，而我们扣掉哪一行，就是以哪个点为根，记结论就好了。

例题：P6178 【模板】Matrix-Tree 定理

这个直接看模板题就好了。

点击查看代码 
#include<bits/stdc++.h>
typedef long long LL;

using namespace std;
const int MAXN=610,P=1e9+7;
int n,m,op;
LL a[MAXN][MAXN],u[MAXN][MAXN],v[MAXN][MAXN];
LL add(LL x,LL y) {
	return (x+y>=P?x+y-P:x+y);
}
int gauss() {
	LL res=1,w=1;
	for(int i=1;i<n;++i) {
		for(int j=i+1;j<n;++j) {
			while(a[i][i]) {
				int div=a[j][i]/a[i][i];
				for(int k=i;k<n;++k) {
					a[j][k]=add(a[j][k],P-1ll*div*a[i][k]%P);
				}
				swap(a[i],a[j]); w=-w;
			}
			swap(a[i],a[j]); w=-w;
		}
	}
	for(int i=1;i<n;++i) {
		res=1ll*a[i][i]*res%P;
	}
	res=1ll*w*res;
	memset(a,0,sizeof(a));
	return (res+P)%P;
}
int main () {
	scanf("%d%d%d",&n,&m,&op);
	for(int i=1;i<=m;++i) {
		int u,v,w;
		scanf("%d%d%d",&u,&v,&w);
		if(u==1) u=n;
		else if(u==n) u=1;
		if(v==1) v=n;
		else if(v==n) v=1;
		if(op==1) {
			a[v][v]=add(a[v][v],w);
			a[u][v]=add(a[u][v],P-w);
		}
		else {
			a[u][u]=add(a[u][u],w);
			a[v][v]=add(a[v][v],w);
			a[u][v]=add(a[u][v],P-w);
			a[v][u]=add(a[v][u],P-w);
		}
	}
	printf("%d\n",gauss());
	return 0;
}

求生成树边权和：

\(\sum\limits_{T}\sum\limits_{e\in T}c_e\)

只要将矩阵的元素换成一次函数，一次项是边权，常数项是 \(1\) 。然后最后乘起来之后，答案就是一次项的系数，具体怎么理解，相当于看一条边在所有生成树中出现的次数然后加起来。

例题：P6624 [省选联考 2020 A 卷] 作业题

点击查看代码 
#include<bits/stdc++.h>
#define pi pair<LL,LL>
#define mp make_pair
#define fi first
#define se second
typedef long long LL;

using namespace std;
const int MAXN=2e5+10,M=110*110,N=40,P=998244353;
int n,m;
struct ddl {
	LL f,t,w;
}a[M];
LL prime[MAXN],cnt,phi[MAXN],sf[MAXN];
LL ksm(LL x,LL y) {
	LL ret=1;
	while(y) {
		if(y&1) ret=ret*x%P;
		x=x*x%P;
		y>>=1;
	}
	return ret;
}
void ycl() {
	phi[1]=1; sf[1]=1;
	for(int i=2;i<=MAXN-10;++i) {
		if(!sf[i]) {
			phi[i]=i-1;
			prime[++cnt]=i;
		}
		for(int j=1;j<=cnt;++j) {
			int t=i*prime[j];
			if(t>MAXN-10) break;
			sf[t]=1;
			if(i%prime[j]==0) {
				phi[t]=phi[i]*prime[j];
				break;
			}
			else phi[t]=phi[i]*phi[prime[j]];
		}
	}
}
LL add(LL x,LL y) {
	return (x+y>=P?x+y-P:x+y);
}
pi b[N][N];
pi operator +(pi a,pi b) {
	return mp(add(a.fi,b.fi),add(a.se,b.se));
}
pi operator -(pi a,pi b) {
	return mp(add(a.fi,P-b.fi),add(a.se,P-b.se));
}
pi operator *(pi a,pi b) {
	return mp(add(a.fi*b.se%P,b.fi*a.se%P),a.se*b.se%P);
}
pi operator /(pi a,pi b) {
	LL inv=ksm(b.se,P-2);
	return mp(add(a.fi*b.se%P,P-a.se*b.fi%P)*inv%P*inv%P,a.se*inv%P);
} 
void operator +=(pi &a,pi b) {
	a=a+b;
}
void operator -=(pi &a,pi b) {
	a=a-b;
}
void operator *=(pi &a,pi b) {
	a=a*b;
}
void operator /=(pi &a,pi b) {
	a=a/b;
}
LL gauss(int n) {
	pi ans=mp(0,1);
	bool rev=0;
	for(int i=1;i<=n;++i) {
		int p=0;
		for(int j=i;j<=n;++j) {
			if(b[j][i].se) {
				p=j; if(i!=j) rev^=1;
				break;
			}
		}
		if(!p) return 0;
		swap(b[i],b[p]);
		pi inv=mp(0,1)/b[i][i];
		for(int j=i+1;j<=n;++j) {
			pi ls=b[j][i]*inv;
			for(int q=i;q<=n;++q) {
				b[j][q]-=ls*b[i][q];
			}
		}
		ans=ans*b[i][i];
	}
	if(rev) ans=mp(0,0)-ans;
	return ans.fi;
}
LL solve(int x) {
	memset(b,0,sizeof(b));
	for(int i=1;i<=m;++i) {
		if(a[i].w%x!=0) continue;
		int t=a[i].f,tt=a[i].t;
		b[t][t]+=mp(a[i].w,1);
		b[tt][tt]+=mp(a[i].w,1);
		b[t][tt]-=mp(a[i].w,1);
		b[tt][t]-=mp(a[i].w,1);
	}
	return gauss(n-1);
}
int main () {
	ycl();
	scanf("%d%d",&n,&m);
	LL w=0;
	for(int i=1;i<=m;++i) {
		scanf("%lld%lld%lld",&a[i].f,&a[i].t,&a[i].w);
		w=max(w,a[i].w);
	}
	LL ans=0;
	for(int i=1;i<=w;++i) {
		int ls=0;
		for(int j=1;j<=m;++j) {
			if(a[j].w%i==0) ++ls;
		}
		if(ls<n-1) continue;
		ans=add(ans,phi[i]*solve(i)%P);
	}
	printf("%lld\n",ans);
	return 0;
}

总结：

其实这个东西就背背公式，套路就好了，没什么好说的。

posted @ 2023-10-22 22:32  daduoli  阅读(87)  评论(0)    收藏  举报