先考慮一個樸素的DP：設 \(f_{i,j}\) 表示第 \(i\) 個位置填 \(1\sim j\) 的所有升序序列對答案的貢獻。轉移方程：

\[f_{i,j}=f_{i-1,j-1}\times j+f_{i,j-1} \]

這樣時間複雜度是 \(O(nk)\) 的，無法接受。

我們先假設 \(f_{i}(j)=f_{i,j}\) 是關於 \(j\) 的一個多項式，次數設為 \(g(i)\)。我們發現狀態轉移方程中一個差分的形式：

\[f_{i,j}-f_{i,j-1}=f_{i-1,j-1}\times j \]

那麼次數的方程：

\[g(i)-1=g(i-1)+1 \]

所以得出：

\[g(i)=g(i-1)+2 \]

顯然 \(g(0)=0\)，那麼以上假設成立（嚴謹證明可使用數學歸納法）：\(f_i(j)\) 是關於 \(j\) 的 \(2\times i\) 次多項式。所以我們只需要計算出 \(f_{i,1}\sim f_{i,2n+1}\) 的值，然後拉格朗日插值即可。

程式碼：

#include<cstdio>
#include<iostream>
#include<algorithm>
#include<cstring>
#define int long long

using namespace std;

const int N=1009;
int k,n,p,f[N][N],pre[N],suf[N],fac[N],inv_fac[N];

void init()
{
	scanf("%lld %lld %lld",&k,&n,&p);
}

int ksm(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1)
			res=res*a%p;
		b>>=1,a=a*a%p;
	}
	return res;
}

int calc(int x,int k)
{
	pre[0]=suf[k+2]=1;
	for (int i=1;i<=k+1;i++)
		pre[i]=pre[i-1]*(x-i)%p;
	for (int i=k+1;i>=1;i--)
		suf[i]=suf[i+1]*(x-i)%p;
	fac[0]=1;
	for (int i=1;i<=k+1;i++)
		fac[i]=fac[i-1]*i%p;
	inv_fac[k+1]=ksm(fac[k+1],p-2);
	for (int i=k;i>=0;i--)
		inv_fac[i]=inv_fac[i+1]*(i+1)%p;
	int ans=0;
	for (int i=1;i<=k+1;i++)
		ans=(ans+f[n][i]*pre[i-1]%p*suf[i+1]%p*inv_fac[i-1]%p*inv_fac[k+1-i]%p*((k+1-i&1)?-1:1))%p;
	return ans;
}

void work()
{
	if(k<=2*n+1)
	{
		for (int i=1;i<=k;i++)
			f[1][i]=i+f[1][i-1];
		for (int i=2;i<=n;i++)
			for (int j=1;j<=k;j++)
				f[i][j]=(f[i][j-1]+f[i-1][j-1]*j%p)%p;
		int tmp=1;
		for (int i=1;i<=n;i++)
			tmp=tmp*i%p;
		printf("%lld\n",tmp*f[n][k]%p);
	}
	else
	{
		int K=2*n+1;
		for (int i=1;i<=K;i++)
			f[1][i]=i+f[1][i-1];
		for (int i=2;i<=n;i++)
			for (int j=1;j<=K;j++)
				f[i][j]=(f[i][j-1]+f[i-1][j-1]*j%p)%p;
		printf("%lld\n",(calc(k,2*n)+p)%p*fac[n]%p);
	}
}

signed main()
{
	init();
	work();
	return 0;
}