【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(n3)，对这道题来说很不优美
所以我们尝试引入拉格朗日插值法
设f(x)是一个n次多项式，则有

f(x)=∑i=0naixi=∑i=0nf(xi)∏j=0nx−xjxi−xj[i≠j]

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

#include<cstdio>#include<iostream>#define LL long longusing 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);}