Description：

$N$ 个点，每个点只能向 $[N-k,N-1]$ 连边，求生成树个数。

$1\le N \le 10^{15}$

Solution：

$DP$ ，记录 $[i-k,i-1]$ 的状态，为了防止重复用最小表示法表示连通性。先搜出来所有状态，在转移时考虑新加进来的一个点和谁连边，如果第一个点没和任何点连边，那他就必须和最后一个点连边，其他点的连边情况 $dfs$ 枚举，用并查集维护每时每刻的连通性，注意不能路径压缩因为回溯时要断。

然后做 $n-k$ 次矩阵乘法，对于刚开始的 $k$ 个点，手算出大小为 $k$ 的联通块生成树个数，最后统计初始情况每种颜色个数乘法原理乘起来再累加到答案里。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<vector>

#include<cstring>

using namespace std;

#define MOD 65521

typedef long long ll;

ll n;

ll cot[6] = {1,1,1,3,16,125};

int d;

int t = 0;

struct matrix

{

	ll d[60][60];

	matrix(){memset(d,0,sizeof(d));}

	void init()

	{

		memset(d,0,sizeof(d));

		for(int i = 1;i <= t;++i)

		{

			d[i][i] = 1;

		}

		return;

	}

	ll *operator [](int x){return d[x];}

	friend matrix operator * (matrix a,matrix b)

	{

		matrix c;

		for(int i = 1;i <= t;++i)

		{

			for(int j = 1;j <= t;++j)

			{

				for(int k = 1;k <= t;++k)

				{

					c[i][j] = (c[i][j] + a[i][k] * b[k][j] % MOD) % MOD;

				}

			}

		}

		return c;

	}

	friend matrix operator ^ (matrix a,ll k)

	{

		matrix res;res.init();

		while(k > 0)

		{

			if(k & 1)res = res * a;

			a = a * a;

			k = k >> 1;

		}

		return res;

	}

}f;

map<vector<int>,int> tr;

struct state

{

	int s[7];

	int& operator [](int x){return s[x];}

	void reset()

	{

		int v[7];memset(v,0,sizeof(v));

		int ptr = 0;

		for(int i = 1;i <= d;++i)

		{

			if(v[s[i]] == 0)

			{

				v[s[i]] = ++ptr;

			}

			s[i] = v[s[i]];

		}

		return;

	}

	vector<int> insert()

	{

		vector<int> l;

		for(int i = 1;i <= d;++i)

		{

			l.push_back(s[i]);

		}

		return l;

	}

}s;

map<int,state> p;

void get(int bit,int cur)

{

	if(bit == d + 1)

	{

		p[++t] = s;

		tr[s.insert()] = t;

		return;

	}

	for(int i = 1;i <= cur + 1;++i)

	{

		s[bit] = i;

		get(bit + 1,max(cur,i));

	}

	return;

}

int fa[7],siz[7];

int find(int x){return (x == fa[x] ? x : find(fa[x]));}

void merge(int p,int q)

{

	p = find(p);q = find(q);

	if(p != q)

	{

		fa[p] = q;

		siz[q] += siz[p];

	}

	return;

}

state nxt;

void dfs(int now,int cur)

{

	if(cur == d + 1)

	{

		for(int i = 2;i <= d + 1;++i)

		{

			nxt[i - 1] = find(i);

		}

		nxt.reset();

		++f[now][tr[nxt.insert()]];

		return;

	}

	dfs(now,cur + 1);

	if(find(cur) == find(d + 1))return;

	int l = find(cur);

	fa[l] = d + 1;

	dfs(now,cur + 1);

	fa[l] = l;

	return;

}

void build()

{

	for(int i = 1;i <= t;++i)

	{

		for(int k = 1;k <= d + 1;++k)

		{

			fa[k] = k;siz[k] = 1;

		}

		for(int k1 = 1;k1 <= d;++k1)

		{

			for(int k2 = 1;k2 <= d;++k2)

			{

				if(p[i][k1] == p[i][k2])

				{

					merge(k1,k2);

				}

			}

		}

		for(int k1 = (siz[find(1)] > 1 ? 0 : 1);k1 <= 1;++k1)

		{

			if(k1 == 1)merge(1,d + 1);

			dfs(i,2);

		}

	}

	return;

}

int main()

{

	cin >> d >> n;

	get(1,0);

	build();

	f = f ^ (n - d);

	ll ans = 0,tmp,cnt;

	for(int i = 1;i <= t;++i)

	{

		tmp = 1;cnt = 0;

		for(int k = 1;k <= d;++k)

		{

			cnt = 0;

			for(int j = 1;j <= d;++j)

			{

				if(p[i][j] == k)

				{

					++cnt;

				}

			}

			tmp = tmp * cot[cnt] % MOD;

		}

		ans = (ans + tmp * f[i][1] % MOD) % MOD;

	}

	cout << ans << endl;

	return 0;

}