BZOJ 1494 NOI2007 生成树计数 状压DP+矩阵乘法

题目大意：给定$n(n\leq 10^{15})$个点，编号差不超过$k(k\leq 5)$的点之间有连边，问生成树个数
将$k$个点的连通性用最小表示法压成状态，那么最多有$52$种状态
计算出每个状态的生成树个数，作为初始行向量$A$
对于每种状态考虑新加入一个点并向这$k$个点连边，每种连法可以转移到哪些状态，得到转移矩阵$B$
那么答案就是$A*B^n[$所有点都连通的状态 ] <script type="math/tex" id="MathJax-Element-331">]</script>

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define M 60
#define MOD 65521
using namespace std;
int k;
long long n;
int pos[1<<15],tot;
namespace Matrix_Multiplication{
    struct Matrix{
        unsigned int a[M][M];
        Matrix(bool flag)
        {
            int i;
            memset(a,0,sizeof a);
            for(i=1;i<=tot;i++)
                a[i][i]=flag;
        }
        unsigned int* operator [] (int x)
        {
            return a[x];
        }
        friend Matrix& operator *= (Matrix &x,Matrix y)
        {
            Matrix z(false);
            int i,j,k;
            for(i=1;i<=tot;i++)
                for(j=1;j<=tot;j++)
                    for(k=1;k<=tot;k++)
                        (z[i][j]+=x[i][k]*y[k][j])%=MOD;
            return x=z;
        }
    }a(false);
    Matrix Quick_Power(Matrix x,long long y)
    {
        Matrix re(true);
        while(y)
        {
            if(y&1) re*=x;
            x*=x; y>>=1;
        }
        return re;
    }
}
void Express(int a[],int sta)
{
    int i;
    for(i=0;i<k;i++)
        a[i]=(sta>>i*3)&7;
}
int Impress(int a[])
{
    int i,re=0;
    for(i=0;i<k;i++)
        re|=a[i]<<i*3;
    return re;
}
void Standardize(int a[])
{
    static int b[10];
    int i,cnt=0;
    memset(b,0,sizeof b);
    for(i=0;i<k;i++)
        if(!b[a[i]])
            b[a[i]]=++cnt;
    for(i=0;i<k;i++)
        a[i]=b[a[i]];
}
bool Is_Valid(int sta)
{
    static int a[10];
    int i,now=1;
    Express(a,sta);
    for(i=0;i<k;i++)
    {
        if(a[i]==0)
            return false;
        if(a[i]>now)
            return false;
        if(a[i]==now)
            ++now;
    }
    return true;
}
void Calculate1(int sta)
{
    static int a[10],b[10],cnt[10];
    int i,_sta,limit=0;
    Express(a,sta);
    memset(cnt,0,sizeof cnt);
    for(i=0;i<k;i++)
    {
        cnt[a[i]]++;
        limit=max(limit,a[i]);
    }
    for(_sta=(cnt[1]==1);_sta<1<<limit;_sta+=(cnt[1]==1)+1)
    {
        int temp=1;
        memcpy(b,a,sizeof b);
        b[k]=9;
        for(i=1;i<=limit;i++)
            if(_sta&(1<<i-1))
                temp*=cnt[i];
        for(i=0;i<k;i++)
            if(_sta&(1<<a[i]-1))
                b[i]=9;
        Standardize(b+1);
        Matrix_Multiplication::a[ pos[sta] ][ pos[Impress(b+1)] ]=temp;
    }
}
int Slow_Power(int x,int y)
{
    int re=1;
    while(y>0)
    {
        re*=x;
        y--;
    }
    return re;
}
void Calculate2(unsigned int A[],int sta)
{
    static int a[10],cnt[10];
    int i,re=1;
    Express(a,sta);
    memset(cnt,0,sizeof cnt);
    for(i=0;i<k;i++)
        cnt[a[i]]++;
    for(i=1;cnt[i];i++)
        re*=Slow_Power(cnt[i],cnt[i]-2);
    A[ pos[sta] ]=re;
}
namespace Union_Find_Set{
    int fa[M];
    void Initialize()
    {
        memset(fa,0,sizeof fa);
    }
    int Find(int x)
    {
        if(!fa[x]||fa[x]==x)
            return fa[x]=x;
        return fa[x]=Find(fa[x]);
    }
}
int main()
{
    using namespace Matrix_Multiplication;
    int i;
    cin>>k>>n;
    if(k>n) k=n;
    for(i=0;i<1<<k*3;i++)
        if( Is_Valid(i) )
            pos[i]=++tot;
    static unsigned int A[M];
    for(i=0;i<1<<k*3;i++)
        if(pos[i])
        {
            Calculate1(i);
            Calculate2(A,i);
        }
    Matrix B=Quick_Power(a,n-k);
    unsigned int ans=0;
    for(i=1;i<=tot;i++)
        (ans+=A[i]*B[i][1])%=MOD;
    cout<<ans<<endl;
    return 0;
}