题目传送门

分析

两点之间的最小距离其实是由两点高度最大值决定的，

求出长度为 \(n\) 的排列所需距离的方案数，剩下还能放的距离可以用插板法放进去。

也就是 \(\sum_{i=1}^{n^2}f_i*\binom{m-i+n}{n}\)

设 \(dp[i][j][k]\) 表示 \(1\sim i\) 被分成 \(j\) 段所需距离为 \(k\) 的方案数。

新开一段就是 \(dp[i][j+1][k+1]+=dp[i-1][j][k]*(j+1)\)（有 \(j+1\) 个位置可以选）

合并到一段开头或结尾就是 \(dp[i][j][k+i]+=dp[i-1][j][k]*(2*j)\)

合并两段就是 \(dp[i][j-1][k+(i*2-1)]+=dp[i-1][j][k]*(j-1)\)（有 \(j-1\) 个位置可以选）

最后的 \(f_k\) 就是 \(dp[n][1][k]\)，时间复杂度 \(O(n^4)\)

考虑一下排列组合怎么求，由于模数不一定是质数，

我一开始以为要质因数分解，写完之后发现其实可以把 \(x+n\) 加进去之后再将 \(x\) 删掉，

删除的过程实际上直接用最大公约数删除就可以了。

这样的时间复杂度是 \(O(n^3\log m)\) 的。

其实还有一种方法是利用杨辉三角的递推公式矩阵加速递推，

然后再用杨辉三角把所有的值都求出来，虽然复杂度相同，不过常数有点大。

代码（矩阵乘法）

#include 
using namespace std;
const int N=111;
struct maix{int p[N][N];}A,ANS;
int dp[N][N*N],n,m,mod,L,R,mid,ans;
void Mo(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
maix mul(maix A,maix B,int t){
	maix C;
	for (int i=0;i<=t;++i)
	for (int j=0;j<=n;++j){
		C.p[i][j]=0;
		for (int k=0;k<=n;++k)
		    Mo(C.p[i][j],1ll*A.p[i][k]*B.p[k][j]%mod);
	}
	return C;
}
int main(){
	ios::sync_with_stdio(0);
	cin>>n>>m>>mod,dp[1][1]=1,mid=(n+1)>>1;
	L=n*(n+1)/2,R=n*(n+3)-mid*(mid+2);
	if (L>m){
		cout<<0;
		return 0;
	}
	if (R>m) R=m;
	for (int i=2;i<=n;++i){
		for (int k=R;k;--k)
		for (int j=1;j1&&k+(i*2-1)<=R) Mo(dp[j-1][k+(i*2-1)],t*(j-1ll)%mod);
			if (k+i<=R) Mo(dp[j][k+i],2ll*t*j%mod);
			if (k>=1,A=mul(A,A,n))
	    if (t&1) ANS=mul(ANS,A,0);
	for (int i=R;i>=L;--i){
		if (dp[1][i]) Mo(ans,1ll*dp[1][i]*ANS.p[0][n]%mod);
		for (int j=n;j;--j) Mo(ANS.p[0][j],ANS.p[0][j-1]);
	}
	cout<

#include 
using namespace std;
const int N=111;
int dp[N][N*N],n,m,mod,L,R,mid,ans,a[N],o=1;
void Mo(int &x,int y){x=x+y>=mod?x+y-mod:x+y;}
int gcd(int x,int y){return y?gcd(y,x%y):x;} 
int main(){
	ios::sync_with_stdio(0);
	cin>>n>>m>>mod,dp[1][1]=1,mid=(n+1)>>1;
	L=n*(n+1)/2,R=n*(n+3)-mid*(mid+2);
	if (L>m){
		cout<<0;
		return 0;
	}
	if (R>m) R=m;
	for (int i=2;i<=n;++i){
		for (int k=R;k;--k)
		for (int j=1;j1&&k+(i*2-1)<=R) Mo(dp[j-1][k+(i*2-1)],t*(j-1ll)%mod);
			if (k+i<=R) Mo(dp[j][k+i],2ll*t*j%mod);
			if (k1;++j){
		int GCD=gcd(a[j],x);
		x/=GCD,a[j]/=GCD;
	}
	for (int i=R;i>=L;--i){
		int now=1;
		for (int j=1;j<=n;++j) now=1ll*now*a[j]%mod;
		if (dp[1][i]) Mo(ans,1ll*dp[1][i]*now%mod);
		a[o]=m-i+n+1,o=o%n+1;
		for (int j=1,x=m-i+1;j<=n&&x>1;++j){
			int GCD=gcd(a[j],x);
			x/=GCD,a[j]/=GCD;
		}
	}
	cout<