【BZOJ2655】calc，dp+拉格朗日插值法

传送门
思路：
好题
这个题目做法非常多，我知道的起码还有倍增法和预处理伯努利数
这里仅说一下我的想法
从原始dp的思路出发
f[i][j]表示i个元素中选取j个的分数
写出dp方程
$f[i][j]=f[i-1][j-1]×i×j+f[i-1][j]$
初始状态$f[0][0]=1$
手玩j比较小的情况，容易发现，$f[i][j]$实际上是一个最高次项为$2j$的多项式
那也就是说我们最终的答案$f[A][n]$就是一个$2n$次的多项式
这是个很好的性质，因为我们只用求出$0$到$2n$次项的系数就可以直接求答案了
那这个系数怎么求呢？
有一种比较容易想到的方法是求出较小的$f[i][n]$，然后高斯消元
但这样的复杂度是$O(n^3)$，对这道题来说很不优美
所以我们尝试引入拉格朗日插值法
设$f(x)$是一个$n$次多项式，则有

$$f(x)=\sum_{i=0}^n a_ix^i=\sum_{i=0}^nf(x_i)\prod_{j=0}^n\frac{x-x_j}{x_i-x_j} [i≠j]$$
其中$x1,x2,..xn$互不相等
同样是求出$2n$个$f[i][n]$然后求解$f[A][n]$，但是显然这个复杂度比高斯消元低多了
理论复杂度为$O(n^2\log A)$，因为还要求解逆元
但实际上是可以做到$O(n^2)$的，因为选取的点都是$0$,$1$,$2$,$3$,$4$..$2n$，随便求和预处理一下就可以$O(n\log A)$了，主要复杂度在初始DP那里，但是常数催人泪下……
UPD
2017.3.25
刚从Shallwe大爷那里了解到，原来这里的插值是可以做到$O(n)$的= =，处理一下阶乘的逆元以及$x-x_i$的前缀后缀积就可以了
不过速度好像并没有提升多少？应该是我常数写挫了

#include<cstdio>
#include<iostream>
#define LL long long
using namespace std;
int A,n,p;
int f[1005][505],inv[1005],pre[1005],sub[1005];
main()
{
    f[0][0]=1;
    scanf("%d%d%d",&A,&n,&p);
    for (int i=1;i<=min(n<<1,A);++i)
        for (int j=0;j<=n;++j)
        if (j)
            f[i][j]=((LL)f[i-1][j-1]*i%p*j%p+f[i-1][j])%p;
        else
            f[i][j]=f[i-1][j];
    if (A<=n*2)
        return printf("%d",f[A][n]),0;
    LL ans=0,Inv,fac;
    inv[0]=inv[1]=1;
    pre[0]=A%p;
    for (int i=1;i<=n<<1;++i) pre[i]=1LL*pre[i-1]*(A-i)%p;
    sub[n<<1]=(A-n*2);
    for (int i=n*2-1;i>=0;--i) sub[i]=1LL*sub[i+1]*(A-i)%p;
    for (int i=2;i<=n<<1;++i) inv[i]=1LL*(p-p/i)*inv[p%i]%p;
    for (int i=2;i<=n<<1;++i) inv[i]=1LL*inv[i]*inv[i-1]%p;
    for (int i=0;i<=n<<1;++i)
    {
        if ((n*2-i)&1) Inv=(-1LL)*inv[i]*inv[n*2-i]%p;
        else Inv=(LL)inv[i]*inv[n*2-i]%p;
        fac=1;
        if (i>0) fac=fac*pre[i-1]%p;
        if (i<n*2) fac=fac*sub[i+1]%p;
        ans=(ans+fac*f[i][n]%p*Inv%p)%p;
    }
    printf("%d\n",(ans+p)%p);
}