6461. 【GDOI2020模拟02.05】生成树（矩阵树及其扩展、二维拉格朗日插值）

题目描述

给定一张 N 个点，M 条边的无向图，边有红、绿、蓝三种颜色，分别用 1,2,3 表示。

求这张图有多少生成树，满足绿色边数量不超过 x，蓝色边数量不超过 y，答案对10^9 + 7 取模。

1 ≤ N ≤ 40,1 ≤ M ≤ 10^5,1 ≤ ci ≤ 3

行列式

定义矩阵A的行列式det(A)或|A|

\(|A|=\sum_{排列p}{(-1)^{p的逆序对个数}\prod{A_{i,p[i]}}}\)

行列式的性质

（\(A\)的转置矩阵\(A^T\)：把\(A\)的行列互换）

①\(|A|=|A^T|\)

②\(|A||B|=|AB|\)

见 https://www.zhihu.com/question/48497108

③\(|A的某一行乘x|=x|A|\)

④\(|A的某行第j项=aj+bj|=|A的某行第j项=aj|+|A的某行第j项=bj|\)

⑤交换相邻两行行列式变号

因为在每一种排列中符号都会改变

⑥交换任意两行行列式变号

两两依次交换的次数为2k+1，根据⑤可得必然变号

⑦有两行完全相同时行列式为0

否则交换两行后行列式改变

⑧把某一行乘k后加到另一行上行列式不变

相当于左乘了一个对角线为1，某个位置为k的矩阵， 这个矩阵的行列式显然是1

根据②可知不变

⑨上三角矩阵的行列式为对角线的积

其实可以不用化成上三角矩阵，只需要消元后n^2算逆序对即可

有了上面这些定理即可把矩阵高斯消元后求得行列式

（注意某行单独乘x时最终的结果要乘上x）

矩阵树定理

余子式：\(M_{i,j}\)表示矩阵\(A\)去掉第i行第j列后剩余矩阵的行列式

基尔霍夫矩阵：最简单的形式即为 度数矩阵D-邻接矩阵C，其中度数矩阵\(D_{i,i}=\)i的度数

其实度数矩阵的真正形态是这样的：\(D_{i,i}=\sum_{i≠j}{C[i][j]}\)

所以当Ci,j不为1，甚至是一个多项式时也是成立的

矩阵树定理：图的生成树个数=基尔霍夫矩阵的任意余子式\(M_{i,i}\)

我也不会证

二维拉格朗日插值

已知x=xi,y=yj时的点值\(a_{x_iy_j}\)，构造\(\sum_{i,j}{a_{x_iy_j}f_{x_i}(x)}g_{x_i}(x)\)使得\(x=x_i\)时\(f_{x_i}(x)=1\)，否则=0，g同理

可以发现和一维的一模一样，\(f_{x_i}(x)=\prod_{j≠i}{\frac{x-x_j}{x_i-x_j}}\)

题解

知道了上面的这题就是模板题了

设三种边边权为1,x,y，那么答案就是\(x^Xy^Y\)（XY为输入的数）的系数

枚举xy，求出行列式后插值即可

code

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define mod 1000000007
#define Mod 1000000005
#define file
using namespace std;

long long A[41][41];
int a[41][41][3];
long long b[41][41];
long long c[41][41];
long long g[41];
long long ans[41][41];
int f[41];
int F[41];
int n,m,X,Y,i,j,k,l,x,y;
long long s,Ans;

long long qpower(long long a,int b)
{
	long long ans=1;
	
	while (b)
	{
		if (b&1)
		ans=ans*a%mod;
		
		a=a*a%mod;
		b>>=1;
	}
	
	return ans;
}

int main()
{
	freopen("tree.in","r",stdin);
	#ifdef file
	freopen("tree.out","w",stdout);
	#endif
	
	scanf("%d%d%d%d",&n,&m,&X,&Y);
	fo(i,1,m)
	scanf("%d%d%d",&j,&k,&l),--l,++a[j][k][l],++a[k][j][l];
	
	fo(x,0,n-1)
	{
		fo(y,0,n-1)
		{
			A[x][y]=1;
			
			fo(i,1,n)
			{
				b[i][i]=0;
				
				fo(j,1,n)
				if (i!=j)
				b[i][j]=-(a[i][j][0]+a[i][j][1]*x+a[i][j][2]*y),b[i][i]-=b[i][j];
			}
			
			memset(f,0,sizeof(f));
			fo(i,1,n-1)
			{
				fo(j,1,n-1)
				if (b[i][j])
				{
					if (!f[j])
					{
						A[x][y]=A[x][y]*b[i][j]%mod;
						
						b[i][j]=qpower(b[i][j],Mod);
						fo(k,j+1,n-1)
						b[i][k]=b[i][k]*b[i][j]%mod;
						b[i][j]=1;
						
						f[j]=i;
						F[i]=j;
						
						fo(k,1,i-1)
						if (F[k]>F[i])
						A[x][y]=-A[x][y];
						break;
					}
					else
					{
						fo(k,j+1,n-1)
						b[i][k]=(b[i][k]-b[i][j]*b[f[j]][k])%mod;
						b[i][j]=0;
					}
				}
				
				if (j>n-1)
				{
					A[x][y]=0;
					break;
				}
			}
		}
	}
	
	fo(i,0,n-1)
	{
		c[i][0]=s=1;
		
		fo(j,0,n-1)
		if (i!=j)
		{
			memset(g,0,sizeof(g));
			s=s*(i-j)%mod;
			
			fo(k,0,n-2)
			{
				g[k]=(g[k]-c[i][k]*j)%mod;
				g[k+1]=(g[k+1]+c[i][k])%mod;
			}
			fo(k,0,n-1) c[i][k]=g[k];
		}
		
		s=qpower(s,Mod);
		fo(j,0,n-1)
		c[i][j]=c[i][j]*s%mod;
	}
	
	fo(i,0,n-1)
	{
		memset(g,0,sizeof(g));
		fo(j,0,n-1)
		{
			fo(k,0,n-1)
			g[k]=(g[k]+c[j][k]*A[i][j])%mod;
		}
		
		fo(j,0,n-1)
		{
			fo(k,0,n-1)
			ans[j][k]=(ans[j][k]+c[i][j]*g[k])%mod;
		}
	}
	
	fo(i,0,X)
	{
		fo(j,0,Y)
		Ans=(Ans+ans[i][j])%mod;
	}
	
	printf("%lld\n",(Ans+mod)%mod);
	
	fclose(stdin);
	fclose(stdout);
	
	return 0;
}