bzoj1494 生成树计数 (dp+矩阵快速幂)

题面欺诈系列...

因为一个点最多只能连到前k个点，所以只有当前的连续k个点的连通情况是对接下来的求解有用的

那么就可以计算k个点的所有连通情况，dfs以下发现k=5的时候有52种。

我们把它们用类似于并查集的方式表达（比如12132代表点1和点3连通，2和5连通，3自己），然后再压缩一下。

但要注意的是，12132和23213这两种实际对应的是一种连通情况，我们只要把它都化成字典序最小的那种就可以了

然后考虑增加一个新点以后状态的转移，可以枚举这个点与前k个点（始状态S）的连边情况，其中有一些是不合法的：

　　1.连到了两个本来就连通的点上（导致成环）

　　2.在1号点不与其他点连通的情况下，没有连到1号点（导致不连通）

然后再根据连边情况得到终状态E，将trans[S][E]++。最后trans[i][j]表示的就是加入一个点后由i状态到j状态的方案数

那么就可以得到递推式$f[i][j]=sum^{总状态数}_{k=1}{f[i-1][k]*trans[k][j]}$，其中f[i][j]表示以i为结尾的k个点状态是j的方案数

那么答案就是f[N][1]（假设1是都连通的状态）

然后初值就应该是f[K][i]=i状态的方案数，其中i状态的方案数为它其中每个联通块方案数的乘积。

那么联通块的方案数怎么算呢？其实题面已经说了...n个点的方案数就是$n^{n-2)}$

然后就可以愉快地dp一脸啦

容易发现递推的形式其实和矩阵乘法是相同的，把f[i]和trans看作矩阵，就是$f[N]=f[K]*trans^{N-K}$

然后就可以倍增做矩阵快速幂了。方法和整数的快速幂是一样的。

复杂度：

　　K=5的状态数接近50，$log(10^{15})$接近50。

　　所以基本上是$50^3*50=6250000$的。

代码写的很蛋疼...

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<map>
#include<vector>
#define LL long long 
#define inf 0x3f3f3f3f
using namespace std;
const LL maxn=1e15+5,maxs=55,p65=50000,mod=65521;

LL rd(){
    LL x=0;char c=getchar();int neg=1;
    while(c<'0'||c>'9'){if(c=='-') neg=-1;c=getchar();}
    while(c>='0'&&c<='9') x=x*10+c-'0',c=getchar();
    return x*neg;
}

LL N;int K,sct,num[6]={0,1,1,3,16,125};
LL f[maxs];
LL trans[2][maxs][maxs],out[2][maxs][maxs];
struct ST{
    int s[6];
    ST(int x=0){memset(s,x,sizeof(s));}
}sta[maxs];
struct Mat{
    LL m[maxs][maxs];
    Mat(){memset(m,0,sizeof(m));}
};
int mp[p65];

void print(ST s){
    for(int i=1;i<=K;i++) printf("%d",s.s[i]);
}
int STtoInt(ST s){
    int x=0;for(int i=1;i<=K;i++) x=x*6+s.s[i];return x;
}

ST toLeast(ST s){
    ST flag=ST(-1),re=ST();int n=0;
    for(int i=1;i<=K;i++){
        if(flag.s[s.s[i]]==-1) flag.s[s.s[i]]=n++;
        re.s[i]=flag.s[s.s[i]];
    }return re;
}

void dfs1(int x,ST s,int m){
    if(x>=K+1){
        sta[++sct]=s;mp[STtoInt(s)]=sct;return;
    }
    for(int i=0;i<=m+1;i++){
        s.s[x]=i;
        dfs1(x+1,s,max(m,i));
    }
}

void dfs2(int x,ST s,ST from){
    if(x>=K+1){
        bool flag[6],bb=s.s[1];ST to=ST();int mi=5;
        memset(flag,0,sizeof(flag));
        //print(from);printf("	");print(s);printf("
");
        for(int i=1;i<=K;i++){
            if(!s.s[i]) continue;
            if(flag[from.s[i]]) return;
            flag[from.s[i]]=1;mi=min(mi,from.s[i]);
        }
        for(int i=2;i<=K;i++){
            to.s[i-1]=flag[from.s[i]]?mi:from.s[i];
            if(from.s[1]==to.s[i-1]) bb=1;
        }to.s[K]=mi;
        if(!bb) return;
        to=toLeast(to);
        //print(from);printf("	");print(s);printf("	");print(to);printf("
");
        trans[0][mp[STtoInt(from)]][mp[STtoInt(to)]]++;
        
    }else{
        dfs2(x+1,s,from);s.s[x]=1;dfs2(x+1,s,from);
    }
}

void sets(){
    for(int i=1;i<=sct;i++){
        dfs2(1,ST(),sta[i]);
        //for(int j=1;j<=sct;j++) printf("%d ",trans[0][i][j]);printf("
");
        ST cnt=ST();
        for(int j=1;j<=K;j++) cnt.s[sta[i].s[j]]++;
        f[i]=1;for(int j=0;j<K;j++) if(num[cnt.s[j]])f[i]*=num[cnt.s[j]];
    }
}

void solve(){
    LL p=N-K;bool b=1,c=1;int i,j,k;
    for(i=1;i<=sct;i++) out[0][i][i]=1;
    while(p){
        if(p&1){
            memset(out[c],0,sizeof(out[c]));
            for(i=1;i<=sct;i++){
                for(j=1;j<=sct;j++){
                    for(k=1;k<=sct;k++){
                        out[c][i][j]=(out[c][i][j]+out[c^1][i][k]*trans[b^1][k][j])%mod;
                    }
                }
            }c^=1;
        }
        memset(trans[b],0,sizeof(trans[b]));
        for(i=1;i<=sct;i++){
            for(j=1;j<=sct;j++){
                for(k=1;k<=sct;k++){
                    trans[b][i][j]=(trans[b][i][j]+trans[b^1][i][k]*trans[b^1][k][j])%mod;
                }
            }
        }p>>=1;b^=1;
    }LL ans=0;
    for(i=1;i<=sct;i++){
        ans=(ans+f[i]*out[c^1][i][1])%mod;
    }printf("%d
",ans);
}

int main(){
    int i,j,k;
    K=rd(),N=rd();
    dfs1(2,ST(),0);sets();
    solve();
    return 0;
}