首先令F(i, j)表示选了i个数字，最大的数字不超过j的方案数，转移显然，即$F(i, j)=F(i,j-1)+i*j*F(i-1,j-1)$。后面那个乘以j就是把这个数字插进去。
考虑把F(i, j)看作是一个函数$f_i(j)$，也就是$f_n(x)$，那么：
$f_n(x)=f_n(x-1)+n(x-1)f_{n-1}(x-1)$
归纳这个问题，当n=1的时候$f_1(x)=\frac{x(x+1)}2$是个二次多项式。
由于是数学归纳，$f_{n-1}(x)$是个2n-2次多项式，乘以x后就是个2n-1次多项式，观察$f_n(x)$其实就是在对$xf_{n-1}(x)$这个2n-1次多项式求前缀和后乘以相当是常数的n，因此求完前缀和的$f_n(x)$就要是个2n次多项式。
本题中n不超过千级别，暴力dp比较小的情况然后插值即可。

#include<iostream>
#include<cstring>
#include<cstdio>
#include<algorithm>
#include<assert.h>
#define lint long long
#define N 1010
using namespace std;
int fac[N],facinv[N],dp[N][N],pre[N],suf[N],p;lint xs[3];
inline int mol(lint x) { return x%=p,(int)(x<0?x+p:x); }
inline int fast_pow(int x,int k,int ans=1)
{
    for(;k;k>>=1,x=mol((lint)x*x))
        if(k&1) ans=mol((lint)ans*x);
    return ans;
}
inline int prelude(int n)
{
    for(int i=fac[0]=1;i<=n;i++) fac[i]=mol((lint)fac[i-1]*i);
    facinv[n]=fast_pow(fac[n],p-2);
    for(int i=n-1;i>=0;i--) facinv[i]=mol(facinv[i+1]*(i+1ll));
    return 0;
}
int main()
{
    int m,n,ans=0;scanf("%d%d%d",&m,&n,&p);
    for(int i=0;i<=2*n+1;i++) dp[0][i]=1;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=2*n+1;j++)
            dp[i][j]=mol(dp[i][j-1]+(lint)i*j*dp[i-1][j-1]);
    if(m<=2*n+1) return !printf("%d\n",dp[n][m]);
    prelude(2*n+1),pre[0]=suf[2*n+2]=1,xs[0]=1ll,xs[1]=-1ll;
    for(int i=1;i<=2*n+1;i++) pre[i]=mol((lint)pre[i-1]*(m-i))%p;
    for(int i=2*n+1;i>=1;i--) suf[i]=mol((lint)suf[i+1]*(m-i))%p;
    for(int i=1;i<=2*n+1;i++)
        ans=mol(ans+xs[(i-1)&1]*dp[n][i]*pre[i-1]%p*suf[i+1]%p*facinv[i-1]%p*facinv[2*n-i+1]);
    return !printf("%d\n",ans);
}