五颜六色的幻想乡 - 矩阵树定理 - 拉格朗日插值

题目大意：给一张图，边有三种颜色，对于每一种可能的a+b+c=n-1的(a,b,c)问恰好a条红色边b条黄色c条蓝色边的方案数，n<=50,5s。
题解：朴素做法，红色视为x黄色视作y蓝色视作1跑矩阵树，直接做是$O(n^6)$，代入x插值算y的多项式可以变成$O(n^5)$，或者更直接的直接将状态压起来，即令$y=x^n$，然后再代x插值即可。本题实现的时候要注意快速幂常数实际上极大，不能乱用。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<vector>
#define lint long long
#define mod 1000000007
#define gc getchar()
#define rep(i,a,b) for(int i=a;i<=b;i++)
#define Rep(i,v) rep(i,0,(int)v.size()-1)
#define N 55
#define M (N*N*3)
#define K (N*N)
#define debug(x) cerr<<#x<<"="<<x
#define sp <<" "
#define ln <<endl
using namespace std;
inline int inn()
{
	int x,ch;while((ch=gc)<'0'||ch>'9');
	x=ch^'0';while((ch=gc)>='0'&&ch<='9')
		x=(x<<1)+(x<<3)+(ch^'0');return x;
}
inline int fast_pow(int x,int k,int ans=1) { for(;k;k>>=1,x=(lint)x*x%mod) (k&1)?ans=(lint)ans*x%mod:0;return ans; }
inline int inv(int x) { return fast_pow(x,mod-2); }
inline int solve(int a[N][N],int n)
{
	int det=1;
	rep(i,1,n)
	{
		int x=0;rep(j,i,n) if(a[j][i]) { x=j;break; }
		if(!x) return 0;if(x^i) det=(mod-det)%mod;swap(a[x],a[i]);
		int v=inv(a[i][i]);det=(lint)det*v%mod;
		rep(j,i,n) a[i][j]=(lint)a[i][j]*v%mod;
		rep(j,1,n) if(i^j) for(int k=n;k>=i;k--)
			a[j][k]-=(lint)a[j][i]*a[i][k]%mod,(a[j][k]<0?a[j][k]+=mod:0);
	}
	return det?inv(det):0;
}
int g[N][N],u[M],v[M],c[M],y[K],fac[K];
inline int getans(int x,int n,int m)
{
	rep(i,1,n) memset(g[i],0,sizeof(int)*(n+1));
	rep(i,1,m)
	{
		int &a=g[u[i]][v[i]],&b=g[v[i]][u[i]],val;
		if(c[i]==1) val=fast_pow(x,n);
		else if(c[i]==2) val=x;else val=1;
		a+=val,b+=val,(a>=mod?a-=mod:0),(b>=mod?b-=mod:0);
	}
	rep(i,1,n) rep(j,1,n) if(i^j) g[i][i]+=g[i][j],(g[i][i]>=mod?g[i][i]-=mod:0);
	rep(i,1,n) rep(j,1,n) if(i^j) g[i][j]=(g[i][j]?mod-g[i][j]:0);
	return solve(g,n-1);
}
vector<int> tot,tmp,ans,res;
inline vector<int> tms(vector<int> x,int v)
{
	res.resize((int)x.size()+1),res[(int)x.size()]=0;
	Rep(i,x) res[i]=(lint)x[i]*v%mod;
	Rep(i,x) res[i+1]+=x[i],(res[i+1]>=mod?res[i+1]-=mod:0);
	return res;
}
inline vector<int> divd(vector<int> x,int v)
{
	int n=(int)x.size();res.resize(n-1);
	for(int i=n-1;i;i--)
		res[i-1]=x[i],x[i-1]-=(lint)v*x[i]%mod,(x[i-1]<0?x[i-1]+=mod:0);
	return res;
}
int main()
{
	int n=inn(),m=inn(),k=n*(n-1);
	rep(i,1,m) u[i]=inn(),v[i]=inn(),c[i]=inn();
	rep(i,1,k+1) y[i]=getans(i,n,m);
	tot.resize(1),tot[0]=1,ans.resize(k+1);
	rep(i,1,k+1) tot=tms(tot,mod-i);
	fac[0]=1;rep(i,1,k+1) fac[i]=(lint)fac[i-1]*i%mod;
	rep(i,1,k+1)
	{
		y[i]=(lint)y[i]*inv((lint)fac[i-1]*fac[k+1-i]%mod)%mod;//dont use bf!That is too slow!
		if((k+1-i)%2==1) y[i]=(mod-y[i])%mod;tmp=divd(tot,mod-i);
		Rep(j,tmp) ans[j]+=(lint)tmp[j]*y[i]%mod,(ans[j]>=mod?ans[j]-=mod:0);
	}
	rep(i,0,n-1) rep(j,0,n-1-i) printf("%d\n",ans[i*n+j]);return 0;
}