LINK：calc
容易得到一个nk的dp做法 同时发现走不通了 此时可以考虑暴力生成函数。
不过化简那套不太熟 且最后需要求多项式幂级数及多项式exp等难写的东西。
这里考虑观察优化dp的做法。
不容易看出 f(n,k)是关于k的2n+1次多项式。
证明可以用数学归纳法证明 且还可以从非常规律的转移中看出这应该是一个形似多项式的东西。
可以直接O(n)拉格朗日插值 不过这里懒得写因为 外面dp是(n^2)求点值的所以这里没必要O(n).
注意初始化.
const ll MAXN=1010;
ll n,mod,k;
ll f[MAXN][MAXN];
ll x[MAXN],y[MAXN];
inline ll ksm(ll b,ll p)
{
ll cnt=1;
while(p)
{
if(p&1)cnt=cnt*b%mod;
p=p>>1;b=b*b%mod;
}
return cnt;
}
inline ll lagrange(ll m,ll z)
{
rep(1,m,i)x[i]=i,y[i]=f[n][i];
ll ans=0;
rep(1,m,i)
{
ll cnt1=1,cnt2=1;
rep(1,m,j)if(i!=j)cnt1=(cnt1*(z-x[j]))%mod,cnt2=(cnt2*(x[i]-x[j]))%mod;
cnt2=ksm(cnt2,mod-2);
cnt1=cnt1*cnt2%mod*y[i]%mod;
ans=(ans+cnt1)%mod;
}
return (ans+mod)%mod;
}
signed main()
{
freopen("1.in","r",stdin);
get(k);get(n);get(mod);ll ans=1;
rep(0,(n<<1|1),j)f[0][j]=1;
rep(1,n,i)
{
ans=ans*i%mod;
rep(1,min(k,(n<<1|1)),j)f[i][j]=(j*f[i-1][j-1]+f[i][j-1])%mod;
}
if(k<=(n<<1|1))putl(ans*f[n][k]%mod);
else putl(lagrange(n<<1|1,k)*ans%mod);
return 0;
}