树形图求和：一道经典矩阵知识题

题目大意

nflsoj题目链接

正睿 题目链接：under权限：18省选十连测

给⼀张(n)个点(m)条有向带权边的有向图。

这张图的⼀个树形图定义为：从(m)条边中选出(n-1)条边，使得所有点都可以通过这些边走到(n)号点。

⼀个树形图的权值定义为这(n-1)条边的权值和。

求出所有树形图的权值和。

答案对(10^9+7)取模。

数据范围：(nleq 300), (mleq 10^5), ( ext{边权}leq 10^9)。

有向图Matrix Tree

我们在小学二年级时都学过无向图生成树计数：Matrix Tree定理。这个定理可以扩展到有向图有根树的计数。

在无向图中，我们定义图的拉普拉斯矩阵是度数矩阵-邻接矩阵。在有向图问题中，如果是求内向树数量，则拉普拉斯矩阵是出度矩阵-邻接矩阵，如果是求外向树数量，则拉普拉斯矩阵是入度矩阵-邻接矩阵。本题显然是前一种情况。

假设我们钦定了一个节点(u)为根（如果题目没有指定，我们就枚举这个(u)），则以(u)为根的生成树数量是拉普拉斯矩阵去掉第(u)行、第(u)列以后的矩阵的行列式. 本题中，(u=n)。

矩阵的行列式可以用带辗转相除法的高斯消元求。代码片段：

inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
struct Matrix{
	int a[305][305];
	Matrix(){memset(a,0,sizeof(a));}
};
int get_det(Matrix A,int n){
	//行列式
	int res=1;
	for(int i=1;i<=n;++i){
		for(int j=i+1;j<=n;++j){
			while(A.a[j][i]){
				int t=A.a[i][i]/A.a[j][i];
				for(int k=1;k<=n;++k)A.a[i][k]=mod2(A.a[i][k]-(ll)A.a[j][k]*t%MOD);
				swap(A.a[i],A.a[j]);res=mod2(-res);
			}
		}
		res=(ll)res*A.a[i][i]%MOD;
	}
	return res;
}

初步转化

本题中，我们考虑每条边对答案的贡献，是：它的权值( imes)它在多少个生成树上出现过. 设原图的生成树数量是(x_0)，去掉第(i)条边以后的生成树数量是(x_i)，则答案就是：

[sum_{i=1}^{m}(x_0-x_i)w_i ]

求(x_i)，可以在原图拉普拉斯矩阵去掉第(n)行、第(n)列得到的矩阵（设为矩阵(A)）的基础上，让(a_{u_i,v_i}+1)，(a_{u_i,u_i}-1)，相当于去掉了第(i)条边，对这个新矩阵求行列式即可。复杂度(O(mn^3))。

优化

矩阵的行列式具有如下性质：

设矩阵(A)的第(k)行满足(forall iin[1,n],a_{k,i}=b_i+c_i). 则：

[|A|= left|egin{matrix} a_{1,1} & a_{1,2} & cdots & a_{1,n}\ vdots & vdots & ddots & vdots\ a_{k,1} & a_{k,2} & cdots & a_{k,n}\ vdots & vdots & ddots & vdots\ a_{n,1} & a_{n,2} & cdots & a_{n,n} end{matrix} ight| = left|egin{matrix} a_{1,1} & a_{1,2} & cdots & a_{1,n}\ vdots & vdots & ddots & vdots\ b_{1} & b_{2} & cdots & b_{n}\ vdots & vdots & ddots & vdots\ a_{n,1} & a_{n,2} & cdots & a_{n,n} end{matrix} ight| + left|egin{matrix} a_{1,1} & a_{1,2} & cdots & a_{1,n}\ vdots & vdots & ddots & vdots\ c_{1} & c_{2} & cdots & c_{n}\ vdots & vdots & ddots & vdots\ a_{n,1} & a_{n,2} & cdots & a_{n,n} end{matrix} ight| ]

（注：(|A|)表示矩阵(A)的行列式）

于是，我们对(A)修改两个位置后求行列式，就相当于三个行列式加起来：

[|newA|=|A|+ left|egin{matrix} a_{1,1} & a_{1,2} & cdots & a_{1,v_i}& cdots & a_{1,n}\ vdots & vdots & ddots & vdots & ddots & vdots\ a_{u_i-1,1} & a_{u_i-1,2} & cdots & a_{u_i-1,v_i}& cdots & a_{u_i-1,n}\ 0 & 0 & cdots & 1& cdots & 0\ a_{u_i+1,1} & a_{u_i+1,2} & cdots & a_{u_i+1,v_i}& cdots & a_{u_i+1,n}\ vdots & vdots & ddots & vdots & ddots & vdots\ a_{n,1} & a_{n,2} & cdots & a_{n,v_i} & cdots & a_{n,n} end{matrix} ight| + left|egin{matrix} a_{1,1} & a_{1,2} & cdots & a_{1,u_i}& cdots & a_{1,n}\ vdots & vdots & ddots & vdots & ddots & vdots\ a_{u_i-1,1} & a_{u_i-1,2} & cdots & a_{u_i-1,u_i}& cdots & a_{u_i-1,n}\ 0 & 0 & cdots & -1& cdots & 0\ a_{u_i+1,1} & a_{u_i+1,2} & cdots & a_{u_i+1,u_i}& cdots & a_{u_i+1,n}\ vdots & vdots & ddots & vdots & ddots & vdots\ a_{n,1} & a_{n,2} & cdots & a_{n,u_i} & cdots & a_{n,n} end{matrix} ight| ]

于是我们相当于要对后面的两个特殊的矩阵求行列式。

余子矩阵

余子式：把一个矩阵的第(i)行、第(j)列去掉以后得到的矩阵的行列式叫(a_{i,j})的余子式，记做(M_{i,j})。

(c_{i,j}=(-1)^{i+j}M_{i,j})是(a_{i,j})的代数余子式。

(A)的每个位置的代数余子式构成的矩阵叫(A)的余子矩阵。

(A)的行列式，除高斯消元外还有一种求法。任取(A)的一行(k)，则(|A|=sum_{i=1}^{n}a_{k,i}c_{k,i})。

发现我们要求行列式的两个矩阵，除(k=u_i)这行外，其它行都与(A)相同，所以它的余子矩阵的这一行，和(A)的余子矩阵这一行是完全相同的。又因为第(k=u_i)行除了一个位置以外其他位置都是(0)。也就是说，我们只要预处理出(A)的余子矩阵(C)，就可以(O(1))算出后面两个新矩阵的行列式。

按定义预处理余子矩阵的复杂度是(O(n^5))，无法接受。但余子矩阵还有别的求法。

定义矩阵(A)的伴随矩阵(operatorname{adj}(A)=displaystyleleft(frac{A}{|A|} ight)^{-1})。

则(A)的余子矩阵(C)是(A)的伴随矩阵的转置矩阵，即(C=operatorname{adj}(A)^{T})。

模拟一下这个过程：先求伴随矩阵，然后求出余子矩阵，再对每条边算贡献。复杂度(O(n^3+m))。

参考代码：

//problem:zr142
#include <bits/stdc++.h>
using namespace std;

#define pb push_back
#define mk make_pair
#define lob lower_bound
#define upb upper_bound
#define fst first
#define scd second

typedef long long ll;
typedef unsigned long long ull;
typedef pair<int,int> pii;
typedef pair<ll,ll> pll;
typedef pair<int,ll> pil;
typedef pair<ll,int> pli;

namespace Fread{
	const int MAXN=1<<20;
	char buffer[MAXN],*S,*T;
	inline char getchar(){
		if(S==T){
			T=(S=buffer)+fread(buffer,1,MAXN,stdin);
			if(S==T) return EOF;
		}
		return *S++;
	}
}
#ifdef ONLINE_JUDGE
	#define getchar Fread::getchar
#endif
inline int read(){
	int f=1,x=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
inline ll readll(){
	ll f=1,x=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}

const int MAXN=1e5+5,MOD=1e9+7;
inline int mod1(int x){return x<MOD?x:x-MOD;}
inline int mod2(int x){return x<0?x+MOD:x;}
inline int pow_mod(int x,int i){int y=1;while(i){if(i&1)y=(ll)y*x%MOD;x=(ll)x*x%MOD;i>>=1;}return y;}
struct E{int u,v,w;}e[MAXN];
struct Matrix{
	int a[305][305];
	Matrix(){memset(a,0,sizeof(a));}
};
//void print(Matrix A,int n){for(int i=1;i<=n;++i){for(int j=1;j<=n;++j)cout<<A.a[i][j]<<" ";cout<<endl;}}
int det,n,m;
int get_det(Matrix A,int n){
	//行列式
	int res=1;
	for(int i=1;i<=n;++i){
		for(int j=i+1;j<=n;++j){
			while(A.a[j][i]){
				int t=A.a[i][i]/A.a[j][i];
				for(int k=1;k<=n;++k)A.a[i][k]=mod2(A.a[i][k]-(ll)A.a[j][k]*t%MOD);
				swap(A.a[i],A.a[j]);res=mod2(-res);
			}
		}
		res=(ll)res*A.a[i][i]%MOD;
	}
	return res;
}
Matrix get_inv(Matrix A,int n){
	//"逆"矩阵
	int b[n+5][n*2+5];memset(b,0,sizeof(b));
	for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)b[i][j]=A.a[i][j];
	for(int i=1;i<=n;++i)b[i][i+n]=1;
	for(int i=1;i<=n;++i){
		int r=-1;
		for(int j=i;j<=n;++j)if(b[i][j]){r=j;break;}
		if(r==-1){puts("0");exit(0);}
		if(r!=i)for(int j=1;j<=n*2;++j)swap(b[i][j],b[r][j]);
		for(int j=1;j<=n;++j)if(i!=j){
			int t=(ll)b[j][i]*pow_mod(b[i][i],MOD-2)%MOD;
			for(int k=1;k<=n*2;++k)b[j][k]=mod2(b[j][k]-(ll)b[i][k]*t%MOD);
		}
		int t=pow_mod(b[i][i],MOD-2);
		for(int j=1;j<=n*2;++j)b[i][j]=(ll)b[i][j]*t%MOD;
	}
	Matrix C;
	for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)C.a[i][j]=b[i][j+n];
	return C;
}
Matrix get_gay(Matrix A,int det,int n){
	//"伴随"矩阵
	int idt=pow_mod(det,MOD-2);
	for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)A.a[i][j]=(ll)A.a[i][j]*idt%MOD;
	return get_inv(A,n);
}
Matrix get_rev(Matrix A,int n){
	//"转置"矩阵
	Matrix B;
	for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)B.a[j][i]=A.a[i][j];
	return B;
}
Matrix get_FishEgg(Matrix A,int n){
	//"鱼子"矩阵
	return get_rev(get_gay(A,det,n),n);
}
int main() {
	n=read(),m=read();
	Matrix A;
	for(int i=1;i<=m;++i)e[i].u=read(),e[i].v=read(),e[i].w=read(),A.a[e[i].u][e[i].v]--,A.a[e[i].u][e[i].u]++;
	for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)A.a[i][j]=mod2(A.a[i][j]);
	det=get_det(A,n-1);
	Matrix B=get_FishEgg(A,n-1);
	int ans=0;
	for(int i=1;i<=m;++i){
		int new_det=0;
		new_det=B.a[e[i].u][e[i].v];
		new_det=mod1(new_det+det);
		new_det=mod2(new_det-B.a[e[i].u][e[i].u]);
		ans=mod1(ans+(ll)e[i].w*mod2(det-new_det)%MOD);
	}
	cout<<ans<<endl;
	return 0;
}