bzoj 1494: [NOI2007]生成树计数

55555555555555555555555555555555被虐哭了

调了两个小时才调出来，简直不要太坑。

论文还是比较良心的，就是初始矩阵好难构造的说。

最后迫不得已写了个很挫的方法构造。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<queue>
using namespace std;
typedef long long ll;
const int p=65521;
const int inf=1e9;
int m;
struct matrix{
	ll a[55][55];
	matrix(){
		memset(a,0,sizeof(a));
	}
	matrix operator * (const matrix &b)const{
		static matrix ans;
		memset(ans.a,0,sizeof(ans.a));
		for(int i=1;i<=m;i++)
		for(int j=1;j<=m;j++)
		for(int k=1;k<=m;k++)
		ans.a[i][j]+=a[i][k]*b.a[k][j];
		for(int i=1;i<=m;i++)
		for(int j=1;j<=m;j++)
		ans.a[i][j]%=p;
		return ans;
	}
	matrix operator ^ (const ll &k)const{
		static matrix ans,b;
		memcpy(b.a,a,sizeof(a));
		for(int i=1;i<=m;i++)
		ans.a[i][i]=1;
		ll h=k;
		for(;h;h>>=1,b=b*b)if(h&1)ans=ans*b;
		return ans;
	}
};
int pa[6],s[55][6];
int k;
void dfs(int i,int tot){
	if(i>k){
		m++;
		memcpy(s[m],pa,sizeof(pa));
	}else{
		for(int j=1;j<=tot;j++){
			pa[i]=j;
			dfs(i+1,tot);
		}
		pa[i]=tot+1;dfs(i+1,tot+1);
	}
}
int size(int i,int j){
	int sz=0;
	for(int t=1;t<=k;t++)
	sz+=s[i][t]==j;
	return sz;
}
bool same(int i,int j,int k){
	return s[k][i]==s[k][j];
}
int mul(int a,int b){
	int ans=1;
	for(int i=1;i<=b;i++)ans*=a; 
	return ans;
}
int size(int i){
	int ans=0;
	for(int j=1;j<=k;j++)
	ans=max(ans,s[i][j]);
	return ans;
}
int calc(int i){
	int sz=size(i),ans=1;
	for(int j=1;j<=sz;j++)
	ans=ans*mul(size(i,j),size(i,j)-2);
	return ans;
}
bool check(int pa1[6],int pa2[6]){
	for(int i=1;i<=k;i++)
	if(pa1[i]!=pa2[i])return false;
	return true;
}
int find(int pa[6]){
	for(int i=1;i<=m;i++)
	if(check(pa,s[i]))return i;
}
bool check(int s,int t){
	for(int i=1;i<=k;i++)
	for(int j=i+1;j<=k;j++)
	if(same(i,j,t)&&(s&(1<<i-1))&&(s&(1<<j-1)))return false;
	return true;
}
matrix trans;
ll f[55];
int main(){
	//freopen("a.in","r",stdin);
	//freopen("a.out","w",stdout);
	ll n;scanf("%d%lld",&k,&n);
	dfs(1,0);
	for(int i=1;i<=m;i++){
		int pa[6];
		for(int state=0;state<(1<<k);state++){
			if(size(i,1)==1&&(!(state&(1<<0))))continue;
			if(!check(state,i))continue;
			memset(pa,0,sizeof(pa));
			int g[7][7];memset(g,0,sizeof(g));
			for(int j=1;j<=k;j++)
			for(int h=j+1;h<=k;h++)
			if(same(j,h,i))
			g[j][h]=g[h][j]=1;
			for(int j=1;j<=k;j++)
			if(state&(1<<j-1))g[k+1][j]=g[j][k+1]=1;
			for(int c=1;c<=k+1;c++)
			for(int a=1;a<=k+1;a++)
			for(int b=1;b<=k+1;b++)
			if(g[a][c]&&g[c][b])g[a][b]=1;
			int tot=0;
			for(int j=2;j<=k+1;j++){
				int h;
				for(h=2;h<j;h++)
				if(g[j][h])break;
				if(j==h)pa[j-1]=++tot;
				else pa[j-1]=pa[h-1];
			}
			int t=find(pa);
			trans.a[i][t]++;
		}
	}
	for(int i=1;i<=m;i++)f[i]=calc(i);
	trans=trans^(n-k);
	ll ans=0;
	for(int i=1;i<=m;i++)
	ans=(ans+f[i]*trans.a[i][1])%p;
	printf("%lld\n",ans);
	return 0;
}