传送门
给定一张 N 个点，M 条边的无向图，边有红、绿、蓝三种颜色，分别用 1,2,3 表示。
求这张图有多少生成树，满足绿色边数量不超过 x，蓝色边数量不超过 y，答案对10^9 + 7 取模。
1 ≤ N ≤ 40,1 ≤ M ≤ 10^5

矩阵树定理

• 专门用来处理无向连通图生成树有关的计数问题。
• 首先定义基尔霍夫（Kirchhoff）矩阵为$度数矩阵A-邻接矩阵E$（都是n阶的）。$A[i][i]$为点i的度数。$E[i][j]=1$表示$i,j$有一条边。
• 定理内容：该矩阵的任意一个代数余子式相等，且生成树的个数即为该矩阵的任意一个代数余子式。（代数余子式即挖掉某一行某一列后剩下的$n-1$阶的矩阵的行列式）。
• 运用高斯消元即可做到$n^3$的时间复杂度。

矩阵树定理的推广

1. 广义上来说求得的应该是所有生成树的边权乘积之和。对于简单的无向图来说边权就是1，就可以套用上面的结论。如果$(i,j)$间有很多条平行边，那么对应的邻接矩阵的位置也就不是1了，而是边的条数，相当于是套了一个权值。
2. 由于求行列式是一系列相乘的计算的，所以我们矩阵的元素还可以是多项式（具体参照这道题目），而求行列式的过程涉及到了多项式乘法。应该注意到的是，度数也要根据组成的不同分成多项式。

Solution

• 首先这道题目显然与矩阵树定理密切相关。
• 而多种颜色的边我们可以用多项式表示，每一对$(i,j)$的边权可以用(a+bx+cy)表示，常数项，x的系数，y的系数分别表示红绿蓝。那么求行列式的多项式乘法最终的答案也是一个多项式，$x^iy^j$的系数自然就是$i$条绿，$j$条蓝边的方案数。
• 感性理解。
• 那么问题在于模数和时间复杂度，使得我们不能直接用多项式乘法，求逆去做高斯消元求行列式。
• 注意到最后次数在n以内，所以我们拿出我们处理多项式问题的又一利器——拉格朗日插值。我们直接暴力$n^2$枚举n种i，n种j，做二维的拉格朗日插值即可。二维的拉格朗日插值与一维的大同小异。
  
• 至此我们需要$n^2$枚举$i,j$，里面一个$n^3$的行列式，以及一个$n^2$的多项式暴力乘法.
• $O(n^5)$

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define maxn 41
#define ll long long 
#define mo 1000000007
using namespace std;

int n,m,i,j,k,cx,cy,du[maxn],x,y,z,a[maxn][maxn][3];
ll inv[maxn],F[maxn][maxn],A[maxn][maxn],X[maxn],Y[maxn],Z[maxn];

void read(int &x){
	x=0; char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar());
	for(;ch>='0'&&ch<='9';ch=getchar()) x=x*10+ch-'0';
}

ll ksm(ll x,ll y){
	ll s=1;
	for(;y;y/=2,x=x*x%mo) if (y&1)
		s=s*x%mo;
	return s;
}

ll Gauss(int x,int y){
	for(i=1;i<n;i++) for(j=1;j<n;j++) 
		A[i][j]=-a[i][j][0]-a[i][j][1]*x-a[i][j][2]*y;
	for(i=1;i<n;i++) {
		A[i][i]=0;
		for(j=1;j<=n;j++) 
			A[i][i]+=a[i][j][0]+a[i][j][1]*x+a[i][j][2]*y;
	}
	ll ans=1;
	for(i=1;i<n;i++) {
		if (!A[i][i]) {
			for(j=i+1;j<n;j++) if (A[j][i]){
				ans=-ans;
				for(k=1;k<n;k++) swap(A[i][k],A[j][k]);
				break;
			}
			if (!A[i][i]) return 0;
		}
		int inv0=ksm(A[i][i],mo-2); ans=ans*A[i][i]%mo;
		for(j=i;j<n;j++) A[i][j]=A[i][j]*inv0%mo;
		for(j=i+1;j<n;j++) for(k=n-1;k>=i;k--)
			(A[j][k]-=A[j][i]*A[i][k])%=mo;
	}
	for(i=1;i<n;i++) ans=ans*A[i][i]%mo;
	return ans;
}

int main(){
	read(n),read(m),read(cx),read(cy);
	for(i=1;i<=n;i++) inv[i]=ksm(i,mo-2);
	for(i=1;i<=m;i++) {
		read(x),read(y),read(z);
		a[x][y][z-1]++,a[y][x][z-1]++;
		du[x]++,du[y]++;
	}
	for(x=1;x<=n;x++) for(y=1;y<=n;y++){
		ll tmp=Gauss(x,y);
		memset(X,0,sizeof(X)),memset(Y,0,sizeof(Y));
		X[0]=Y[0]=1;
		for(k=0,i=1;i<=n;i++) if (i!=x){
			tmp=tmp*inv[abs(x-i)]%mo*((x>i)?1:-1);
			for(j=++k;j>=1;j--) X[j]=(-X[j]*i+X[j-1])%mo;
			X[0]=-X[0]*i%mo;
		}
		for(k=0,i=1;i<=n;i++) if (i!=y){
			tmp=tmp*inv[abs(y-i)]%mo*((y>i)?1:-1);
			for(j=++k;j>=1;j--) Y[j]=(-Y[j]*i+Y[j-1])%mo;
			Y[0]=-Y[0]*i%mo;
		}
		for(i=0;i<n;i++) for(j=0;j<n;j++)
			F[i][j]+=X[i]*Y[j]%mo*tmp%mo;
	}
	for(i=0;i<n;i++) for(j=0;j<n;j++)
		F[i][j]=(F[i][j]%mo+mo)%mo;
	ll ans=0;
	for(i=0;i<=cx;i++) for(j=0;j<=cy;j++)
		ans+=F[i][j]%mo;
	printf("%lld",(ans%mo+mo)%mo);
}