count （类插头DP+矩阵快速幂）

题目大意：有n个点，编号为1~n。第i个点和第j个点之间有一条无向边当且仅当|i-j|<=k。求这个图的生成树个数。$k\leq 5,n\leq 10^{15}$。

---

题目分析：Coming在他初二时的资料里找到的一道题，是我校上古大神cdc给的。我不得不吐槽：难道前几届的dalao初二就能做这种题了吗？而且题面还很恶意地给出了怎么用矩阵树定理算无向图的生成树个数……

言归正传，这题很明显就是状态压缩，然后扔进矩乘。首先认为所有边都是从编号大的点连向编号小的点。对于第i个点，用状压记录包括它在内前面k个点的连通块归属情况。由于用最小表示法记录（最小表示法就是将4 2 1 0 0变成0 1 2 3 3这种，尽量让越前面的数越小，相同的数字替换后依旧相同），状态数在k=5的时候也只有52。然后暴力建出矩阵再快速幂即可。

关于这题的具体实现以及转移矩阵的构造，我个人YY出一种比较方便的写法：首先对所有小于等于k位的（不足补0）k+1进制下的数s，都算出它替换成最小表示法后的结果Min[s]。如果Min[s]=s，这就是合法状态，于是在矩阵中为它开出一位并记录下来。容易发现合法状态中一定不会出现某一位的值为k。然后扫一遍所有合法状态，再用$2^k$枚举第i+1个点是否向i~i-k+1连边。如果连了两个已经属于相同连通块的点则不合法；如果没有点和第i-k+1个点属于同一个连通块，而i+1又没有向它连边则不合法。确定连边情况合法后，将i+1涉及到的所有连通块的点的值全部变为k，然后调用Min数组即可得到其最小表示。

最后说一下初始状态的权值。例如初始状态为0 1 2 1 2，则它的权值为$f[1]*f[2]*f[2]$，其中f[i]表示i个点的完全图生成树个数。由于本题只需要算到f[5]，可以手玩可以暴力。不过根据我高一时Semiwaker讲树的prufer序列之类的知识，$f[n]=n^{n-2}$。

---

CODE：

#include<iostream>
#include<string>
#include<cstring>
#include<cmath>
#include<cstdio>
#include<cstdlib>
#include<stdio.h>
#include<algorithm>
using namespace std;

const int maxn=60;
const int maxs=20000;
const long long M=65521;
typedef long long LL;

struct mat
{
    LL val[maxn][maxn];
} e,a;
int num[maxn];
LL pre[maxn];
int N=0;

int Min[maxs];
int id[maxs];

int temp[10];
int cnt[10];
LL f[10];

int ms,k;
LL n;

mat Times(mat x,mat y)
{
    mat z;
    for (int i=0; i<N; i++)
        for (int j=0; j<N; j++) z.val[i][j]=0;
    for (int i=0; i<N; i++)
        for (int j=0; j<N; j++)
            for (int w=0; w<N; w++)
                z.val[i][j]=(z.val[i][j]+ x.val[i][w]*y.val[w][j] )%M;
    return z;
}

int main()
{
    freopen("count.in","r",stdin);
    freopen("count.out","w",stdout);

    scanf("%d%I64d",&k,&n);
    if (k==1)
    {
        printf("1\n");
        return 0;
    }

    ms=1;
    for (int i=1; i<=k; i++) ms*=(k+1);
    for (int i=0; i<ms; i++)
    {
        int s=i;
        for (int j=1; j<=k; j++) temp[j]=s%(k+1),s/=(k+1);
        int x=-1;
        for (int j=0; j<=k; j++) cnt[j]=-1;
        for (int j=1; j<=k; j++)
            if (cnt[ temp[j] ]!=-1) temp[j]=cnt[ temp[j] ];
            else cnt[ temp[j] ]=++x,temp[j]=x;
        s=0;
        for (int j=k; j>=1; j--) s=s*(k+1)+temp[j];
        Min[i]=s;
    }

    for (int i=0; i<ms; i++) id[i]=-1;
    for (int i=0; i<ms; i++)
        if (Min[i]==i)
        {
            id[i]=N;
            num[N]=i;
            N++;
        }

    int mt=(1<<k);
    for (int i=0; i<N; i++)
        for (int t=0; t<mt; t++)
        {
            int s=num[i];
            for (int j=1; j<=k; j++) temp[j]=s%(k+1),s/=(k+1);

            bool sol=true;
            for (int j=1; j<k; j++) if (t&(1<<(j-1)))
                for (int w=j+1; w<=k; w++) if (t&(1<<(w-1)))
                    if (temp[j]==temp[w]) sol=false;
            if (!sol) continue;

            bool zero=false;
            for (int j=2; j<=k; j++) if (!temp[j]) zero=true;
            if ( !zero && !(t&1) ) continue;

            for (int j=0; j<=k; j++) cnt[j]=0;
            for (int j=1; j<=k; j++) if (t&(1<<(j-1))) cnt[ temp[j] ]=-1;
            for (int j=1; j<=k; j++) if (cnt[ temp[j] ]==-1) temp[j]=k;

            s=0;
            for (int j=k; j>=1; j--) s=s*(k+1)+temp[j];
            s/=(k+1);
            s+=(k*ms/(k+1));

            s=id[ Min[s] ];
            a.val[i][s]++;
        }

    f[0]=1;
    f[1]=1;
    f[2]=1;
    f[3]=3;
    f[4]=16;
    f[5]=125; //f[n]=n^(n-2)

    for (int i=0; i<N; i++)
    {
        int s=num[i];
        for (int j=1; j<=k; j++) temp[j]=s%(k+1),s/=(k+1);

        for (int j=0; j<=k; j++) cnt[j]=0;
        for (int j=1; j<=k; j++) cnt[ temp[j] ]++;

        pre[i]=1;
        for (int j=0; j<=k; j++) pre[i]*=f[ cnt[j] ];
    }

    for (int i=0; i<N; i++) e.val[i][i]=1;
    n-=k;
    mat b=e;
    for (int i=61; i>=0; i--)
    {
        b=Times(b,b);
        if (n&(1LL<<i)) b=Times(b,a);
    }

    LL ans=0;
    for (int i=0; i<N; i++) ans=(ans+ pre[i]*b.val[i][0] )%M;
    printf("%I64d\n",ans);

    return 0;
}