【BZOJ1494】生成树计数（NOI2007）-连通性状压DP+矩阵加速

测试地址：生成树计数
做法：本题需要用到连通性状压DP+矩阵加速。
假设图是从左到右不断添加点的，那么每添加一个点，影响到的只有最后k$k$个点的连通性，所以不难想到以最后k$k$个点的连通性作为状态进行DP。所谓的连通性，就是将这k$k$个点分为若干个集合，相同集合内的点就是连通的，不同集合内的点就不连通，不难发现合法的状态数只有B(k)$B(k)$，即第k$k$个Bell数，而Bell数就是k$k$个元素划分集合的方案数，打表可知B(5)=52$B(5)=52$。
那么我们怎么把这些状态表示出来呢？一般我们为了确保相同的连通情况不重复出现，我们用最小表示法表示一个集合划分方案，即从1$1$到n$n$遍历每个点，如果它没有被编过号，那么就将它和所有和它相连的点编上同一个号，就这样继续下去。容易发现这种状态可以用一个k$k$进制数表示，那么给状态编码和解码的方法就都有了。
接下来我们来想怎么DP。注意到，每加入一个点，我们要做出的决策是：在这个点向左连接的k$k$条边中选择若干条边（可以不选）连接，每种决策都会转移到新的最后k$k$个点的连通状态。我们显然可以暴力模拟加入边的过程，存储一个连通状态可以转移到哪些连通状态。
进行决策的时候有以下几个限制：
一是如果两个点已经连通，那么它们不能通过新的点再相连，这是为了维护树的无环性。
二是添加新点之前，最后k$k$个点中的最前面的点要么跟其它的k−1$k-1$个点连通，要么就必须跟新点连通，这是为了维护树的连通性。
在注意上面两点的情况下，直接DP就好了。注意到DP中每一步的转移都是一样的，显然可以用矩阵快速幂加速这个过程，那么总的时间复杂度就是O(B(k)3logn)$O(B(k{)}^3\log n)$，可以通过此题。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=65521;
int k,cnt[10]={0},st[50010],tot=0;
int pos[50010],nowst[10],nxt[10],a[10];
ll n,full[10],first[50010];
bool vis[10],visst[10];
struct matrix
{
    ll mat[60][60];
}M;

void find_state(int step,int last,int num)
{
    if (step==k)
    {
        st[++tot]=num;
        pos[num]=tot;
        first[tot]=1;
        for(int i=1;i<=k;i++)
            first[tot]*=full[cnt[i]];
        return;
    }
    for(int i=1;i<=last;i++)
    {
        cnt[i]++;
        if (i==last) find_state(step+1,last+1,num*6+i-1);
        else find_state(step+1,last,num*6+i-1);
        cnt[i]--;
    }
}

void init()
{
    scanf("%d%lld",&k,&n);
    full[0]=1;
    for(int i=1;i<=k;i++)
    {
        full[i]=1;
        for(int j=1;j<=i-2;j++)
            full[i]*=i;
    }
    find_state(0,1,0);
}

int encode()
{
    int num=0;
    for(int i=1;i<=k;i++)
        num=num*6+nxt[i]-1;
    return num;
}

void decode(int num)
{
    for(int i=k;i>=1;i--)
    {
        nowst[i]=num%6+1;
        num/=6;
    }
}

void pre_dp()
{
    memset(M.mat,0,sizeof(M.mat));
    for(int i=1;i<=tot;i++)
    {
        decode(st[i]);

        bool flag=0;
        int cnt=0;
        for(int j=1;j<=k;j++)
        {
            if (j>1&&nowst[j]==1) flag=1;
            cnt=max(cnt,nowst[j]);
        }

        for(int j=0;j<(1<<k);j++)
        {
            if (!flag&&!(j&1)) continue;

            int mn=cnt+1;
            memset(visst,0,sizeof(visst));
            memset(vis,0,sizeof(vis));
            for(int p=0;p<k;p++)
            {
                vis[p]=(j&(1<<p));
                visst[nowst[p+1]]=visst[nowst[p+1]]|(j&(1<<p));
                if (vis[p]) mn=min(mn,nowst[p+1]);
            }

            bool nowflag=0;
            for(int p=1;p<=k;p++)
            {
                for(int q=1;q<=k;q++)
                    if (p!=q&&vis[p-1]&&vis[q-1]&&nowst[p]==nowst[q])
                        {nowflag=1;break;}
                if (nowflag) break;
            }
            if (nowflag) continue;

            nxt[k]=mn;
            for(int p=1;p<k;p++)
            {
                nxt[p]=nowst[p+1];
                if (visst[nowst[p+1]]) nxt[p]=mn;
            }
            int nxtcnt=0;
            memset(a,0,sizeof(a));
            for(int p=1;p<=k;p++)
            {
                if (!a[nxt[p]]) a[nxt[p]]=++nxtcnt;
                nxt[p]=a[nxt[p]];
            }

            int nx=encode();
            M.mat[pos[nx]][i]=(M.mat[pos[nx]][i]+1)%mod;
        }
    }
}

void mult(matrix A,matrix B,matrix &C)
{
    memset(C.mat,0,sizeof(C.mat));
    for(int i=1;i<=tot;i++)
        for(int j=1;j<=tot;j++)
            for(int k=1;k<=tot;k++)
                C.mat[i][j]=(C.mat[i][j]+A.mat[i][k]*B.mat[k][j])%mod;
}

void dp(ll n)
{
    matrix S;
    memset(S.mat,0,sizeof(S.mat));
    for(int i=1;i<=tot;i++)
        S.mat[i][i]=1;
    while(n)
    {
        if (n&1) mult(M,S,S);
        mult(M,M,M);n>>=1;
    }
    ll ans=0;
    for(int i=1;i<=tot;i++)
        ans=(ans+S.mat[1][i]*first[i])%mod;
    printf("%lld",ans);
}

int main()
{
    init();
    pre_dp();
    dp(n-(ll)k);

    return 0;
}