題目大意：給定n(n≤1015)個點，編號差不超過k(k≤5)的點之間有連邊，問生成樹個數
將k個點的連通性用最小表示法壓成狀態，那麼最多有52種狀態
計算出每個狀態的生成樹個數，作為初始行向量A
對於每種狀態考慮新加入一個點並向這k個點連邊，每種連法可以轉移到哪些狀態，得到轉移矩陣B
那麼答案就是A∗Bn[所有點都連通的狀態]

#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;
}