[NOI2007]生成树计数（状压dp+矩阵加速）

【题解】

考虑动态规划的状态转移：

从左往右，加入第n个结点时，只考虑它向前连的边，那么答案就与向前连哪些边，以及之前的n-1个点构成的答案有关 
然而，"向前连哪些边"并不是有2^k种情况 
前n个点构成了许多连通块，而非一棵树 
因此，结点n引出的边既要使结点1~n连通，也不能构成环 
而n只能向 n-k ~ n-1 连边 
所以n的答案只与结点n-k~n-1的连通性有关 

用最小表示法表示出 点n-k~n-1 属于哪个连通块。如k=5时，00121表示不加点n的情况下 n-5与n-4连通，n-3与n-1连通。
情况12114用最小表示法表示为：01002
由于k<=5，可以预处理出所有不超过52种情况，并把它压入十进制数S中作为后k个点的状态集合
预处理的方法为：dfs，规定每个数x之前必须出现不小于 x-1的数 

所以，可以用f[n][p]来表示状态，S[p]为一个不超过5位的十进制数 
则f[n][i]=sigma(f[n-1][j]*g[j][i])，其中g[j][i]代表：点n-k~n-1状态为S[j] 转移到 点n-k+1~n状态为S[i] 的连边方案数 
可以预处理g[j][i]：通过 枚举j、枚举连边方案t(一个k为二进制数)来推出i，若方案t可行，g[j][i]++
*方案t可行 ，需满足：
 1.使结点1~n连通：可能点n-1与n通过n后面的点相连通，但点n-k若与 n-k+1~n-1均不联通，则 n-k 必须与n相连 
 2.没有构成环：t没有将 前面k个点中已连通的点 同时连接 

预处理 f[k][p]：枚举所有状态p，乘法原理即可。用到的结论：n个结点的完全图的生成树个数为n^(n-2)

如此可以O(n*p*p)得出答案 
对于最后两个测试点(n<=10^15)：将g数组写成矩阵，加速之即可 

注意：
1.递归形式的快速幂会爆空间，应该为迭代形式 

2.别忘记取模 

【代码】

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#define MOD 65521
typedef long long LL;
LL f[60],count[10];
int S[60],num[10],link[10],change[10],shi[10];
int k,p=0;
struct juzhen
{
	LL s[60][60];
	juzhen()
	{
		memset(s,0,sizeof(s));
	}
};
juzhen DW;
juzhen cheng(juzhen a,juzhen b)
{
	juzhen res;
	int i,j,l;
	for(i=1;i<=p;i++)
		for(j=1;j<=p;j++)
		{
			for(l=1;l<=p;l++)
				res.s[i][j]+=a.s[i][l]*b.s[l][j];
			res.s[i][j]%=MOD;
		}
	return res;
}
juzhen ksm(juzhen a,LL n)
{
	juzhen ans=DW,t=a;
	for(;n>0;n=n>>1)
	{
		if(n&1) ans=cheng(ans,t);
		t=cheng(t,t);
	}
	return ans;
}
int max(int a,int b)
{
	if(a>b) return a;
	return b;
}
void init()
{
	int i;
	for(i=1;i<60;i++)
		DW.s[i][i]=1;
	shi[0]=1;
	shi[1]=10;
	shi[2]=100;
	shi[3]=1000;
	shi[4]=10000;
	shi[5]=100000;
	count[0]=count[1]=count[2]=1;
	count[3]=3;
	count[4]=16;
	count[5]=125;
}
void getS(int i,int Max)
{
	int j;
	if(i>k)
	{
		p++;
		for(j=1;j<=k;j++)
			S[p]=S[p]*10+num[j];
		return;
	}
	for(j=0;j<=Max;j++)
	{
		num[i]=j;
		getS(i+1,max(Max,j+1));
	}
}
int main()
{
	juzhen g;
	LL n,I,ans=0;
	int i,j,l,t,t2;
	scanf("%d%lld",&k,&n);
	init();
	getS(1,0);
	for(i=1;i<=p;i++)//枚举连续k位的所有状态 
	{
		memset(num,0,sizeof(num));
		for(j=0;j<k;j++)
			num[ (S[i]/shi[j]) % 10 ]++;
		for(j=0;j<=(1<<k)-1;j++)//枚举连接方式 
		{
			memset(link,0,sizeof(link));
			t=S[i];
			for(l=j;l>0;l=l>>1)
			{
				link[t%10]+=l&1;
				t/=10;
			}
			for(l=0;l<k;l++)//排除不合法情况 
				if(num[l]>1&&link[l]>1) break;
			if(l<k) continue;
			if(num[0]==1&&link[0]==0) continue;
			for(t=0;link[t]==0&&t<k;t++);//计算新状态 
			t2=0;
			for(l=k-1;l>=0;l--)
			{
				if( link[ (S[i]/shi[l])%10 ] > 0 ) t2=t2*10+t;
				else t2=t2*10+(S[i]/shi[l])%10;
			}
			t2=t2*10+t;
			for(l=0;l<=k;l++)//用"最小表示法"表示新状态 
				change[l]=-1;
			t=-1;
			for(l=k-1;l>=0;l--)//change:对应法则 
				if( change[ (t2/shi[l]) % 10 ]==-1 ) change[ (t2/shi[l]) % 10 ]=++t;
			t=0;
			for(l=k-1;l>=0;l--)
				t=t*10+change[(t2/shi[l])%10];
			for(l=1;S[l]!=t;l++);
			g.s[i][l]++;
		}
	}
	for(i=1;i<=p;i++)
	{
		memset(num,0,sizeof(num));
		for(j=0;j<k;j++)
			num[ (S[i]/shi[j]) % 10 ]++;
		f[i]=1;
		for(j=0;j<k;j++)
			f[i]*=count[num[j]];
	}
	g=ksm(g,n-k);
	for(i=1;i<=p;i++)
		ans+=f[i]*g.s[i][1];
	printf("%lld",ans%MOD);
	return 0;
}