BZOJ2655 calc（动态规划+拉格朗日插值法）

　　考虑暴力dp：f[i][j]表示i个数值域1~j时的答案。考虑使其值域++，则有f[i][j]=f[i][j-1]+f[i-1][j-1]*i*j，边界f[i][i]=i!*i!。

　　注意到值域很大，考虑能不能在这一维上优化。完全不会证地有f[i][j]是一个关于j的2i次多项式。那么dp出一部分后就可以直接拉格朗日插值求出多项式，代入即可。

#include<iostream> 
#include<cstdio>
#include<cmath>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;
int read()
{
    int x=0,f=1;char c=getchar();
    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}
    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();
    return x*f;
}
#define N 510
int n,m,P,f[N][N<<1],fac[N],ans=0;
int ksm(int a,int k)
{
    if (k==0) return 1;
    int tmp=ksm(a,k>>1);
    if (k&1) return 1ll*tmp*tmp%P*a%P;
    else return 1ll*tmp*tmp%P;
}
int inv(int x){return ksm(x,P-2);}
int main()
{
#ifndef ONLINE_JUDGE
    freopen("bzoj2655.in","r",stdin);
    freopen("bzoj2655.out","w",stdout);
    const char LL[]="%I64d
";
#else
    const char LL[]="%lld
";
#endif
    m=read(),n=read(),P=read();
    fac[0]=1;for (int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%P;
    for (int i=0;i<=n;i++)
    {
        f[i][i]=1ll*fac[i]*fac[i]%P;
        for (int j=i+1;j<=(n<<1);j++)
        f[i][j]=(f[i][j-1]+1ll*f[i-1][j-1]*i%P*j%P)%P;
    }
    for (int i=0;i<=(n<<1);i++)
    {
        int w=f[n][i],v=1;
        for (int j=0;j<=(n<<1);j++)
        if (i!=j) w=(1ll*w*(m-j+P)%P)%P,v=1ll*v*(i-j+P)%P;
        ans=(ans+1ll*w*inv(v))%P;
    }
    cout<<ans;
    return 0;
}