FFT 【JSOI2012】bzoj4332 分零食

题目大意：

有n个小朋友，m块糖。
给小朋友分糖，如果一个小朋友分不到糖，那他后面的小朋友也分不到糖。
每个小朋友有一个喜悦值，有三个参数，O，S，U，设一个小朋友分到糖数为x，则这个小朋友的喜悦值为O*x x+ S x +U，分不到糖的小朋友的喜悦值为1。
求所有分糖方案下 所有小朋友喜悦值乘积 的和。

题目分析：

如果没有小朋友必须是前面的连续的一段这个要求，那就是一个FFT模板，把喜悦值数组自乘n次即可。
但是有这个要求之后就不能这么做了。

这个时候我们考虑DP。
设数组g［ｉ］［ｊ］为给前ｉ个小朋友分ｊ块糖（每个小朋友都能分到糖）的答案。
设数组ｆ［ｉ］为一个小朋友分到ｉ块糖的喜悦值。

那么转移方程就可以写成ｇ［ｉ］［ｊ］=∑（1<=k<=m）g[ｉ－１］［ｊ－ｋ］＊ｆ［ｋ］。
这样就是标准的卷积的形式，就可以用快速傅里叶变换来优化，可以用快速幂迅速求出g[min(n,m)][m]，括弧但是我们并没有用到快速幂括回笑。

我们求出g[min(n,m)][m]并没有什么用,因为不一定会分到第n个小朋友。

所以我们要求的是∑（1<=i<=min(n,m)）g[i][m]。
设数组p［ｉ］［ｊ］＝∑（１＜＝ｋ＜＝ｉ）ｇ［ｋ］［ｊ］。
我们只要求出ｐ［ｍｉｎ（ｎ，ｍ）］［ｍ］即可。

ｐ可以递推来求
当　ｉ　ｍｏｄ　２　为０时　，ｐ［ｉ］＝ｐ［ｉ／２］＋ｐ［ｉ／２］＊ｇ［ｉ／２］；
当　ｉ　ｍｏｄ　２　为１时　，ｐ［ｉ］＝ｐ［ｉ／２］＋ｐ［ｉ／２］＊ｇ［ｉ／２］＋ｇ［ｉ］；

这个可以感性的证明一下：
首先ｇ［ａ＋ｂ］＝ｇ［ａ］＊ｇ［ｂ］。
给ａ＋ｂ个小朋友分糖，就相当于给前ａ个小朋友分糖和给ａ＋１到ａ＋ｂ个小朋友分糖的答案求一个卷积，在点值表达式的意义下计算是合法的。

那么ｐ［ｉ／２］为ｇ［１］到ｇ［ｉ／２］的和，把ｐ［ｉ／２］乘以ｇ［ｉ／２］，那我们就得到了ｇ［ｉ／２＋１］到ｇ［ｉ］的和，再加上ｐ［ｉ／２］，就是所求的ｐ［ｉ］。

当ｉ为奇数时，我们就先求ｐ［ｉ－１］，然后再加上ｇ［ｉ］就可以了，括弧还是挺好理解的吧括回笑。

设ｔｏｐ＝ｍｉｎ（ｎ，ｍ）．
我们要求ｐ［ｔｏｐ］和ｇ［ｔｏｐ］，那就要求出ｐ［ｔｏｐ／２］　和　ｇ［ｔｏｐ／２］。
那就可以递归处理啦。

时间复杂度为Ｏ（ｍｌｏｇｍｌｏgｍｉｎ（ｎ，ｍ））。

此篇题解部分借鉴自sengxian’s blog.

注意：

1. 亘古不变的精度问题！！！！
2. 在乘的时候不要连乘啊，要不然会变成循环卷积啊什么的云云，会错，就是说每次在做乘法的时候先DFT过去，乘完之后立刻IDFT回来，然后把超过m的部分全部去掉，再进行下一次乘法。

代码：

#include <cstdio>#include <cmath>#include <algorithm>#include <iostream>#define N 1200000using namespace std;const double pi=acos(-1);const double DFT=2.0;const double IDFT=-2.0;long long n,m,O,S,U,mod;int len;int pos[N];struct complex{    double a,b;    complex(double a=0,double b=0):a(a),b(b){}    complex operator + (const complex &c) { return complex(a+c.a,b+c.b); }    complex operator - (const complex &c) { return complex(a-c.a,b-c.b); }    complex operator * (const complex &c) { return complex(a*c.a-b*c.b,a*c.b+b*c.a); }}g[N],p[N],f[N],tmp[N];void init(int len){    for(int i=1;i<len;i++)    {        pos[i]=pos[i>>1]>>1;        if(i&1) pos[i]|=len>>1;    }    return;}void Fast_Fourier_Transform(complex x[],int len,double mode){    for(int i=1;i<len;i++)    if(i<pos[i]) swap(x[i],x[pos[i]]);    for(int i=2;i<=len;i<<=1)    {        int step=i>>1;        complex wm=complex(cos(2*pi/i),sin(mode*pi/i));        for(int j=0;j<len;j+=i)        {            complex w=complex(1,0);            for(int k=j;k<j+step;k++)            {                complex a=x[k],b=w*x[k+step];                x[k]=a+b;x[k+step]=a-b;                w=w*wm;            }        }    }    if(mode==IDFT)        for(int i=0;i<len;i++) x[i].a=(x[i].a)/len;    return;}#define FFT Fast_Fourier_Transformvoid modify(complex x[]){    for(int i=0;i<len;i++)    {        if(i<=m) x[i].a=((long long)(x[i].a+0.5)%mod+mod)%mod,x[i].b=0;        else x[i]=complex(0,0);    }}void solve(long long top){    if(top==1)    {        for(int i=0;i<len;i++) p[i]=f[i],g[i]=f[i];        FFT(p,len,IDFT); FFT(g,len,IDFT);        modify(p); modify(g);        return;    }    solve(top/2);    FFT(p,len,DFT); FFT(g,len,DFT);    for(int i=0;i<len;i++) p[i]=p[i]+p[i]*g[i];    for(int i=0;i<len;i++) g[i]=g[i]*g[i];    FFT(p,len,IDFT); FFT(g,len,IDFT);    modify(p); modify(g);    if(top&1)    {        FFT(g,len,DFT); FFT(p,len,DFT);        for(int i=0;i<len;i++) g[i]=g[i]*f[i];        for(int i=0;i<len;i++) p[i]=p[i]+g[i];        FFT(g,len,IDFT); FFT(p,len,IDFT);        modify(p); modify(g);    }    return;}int main(){    scanf("%lld%lld",&m,&mod);    scanf("%lld",&n);    scanf("%lld%lld%lld",&O,&S,&U);    for(int i=1;i<=m;i++) f[i].a=(O*i*i+S*i+U)%mod;    long long top=min(n,m);    len=1;    while(len<=(m<<1)) len<<=1;    init(len);    FFT(f,len,DFT);    solve(top);    printf("%lld\n",((long long)(p[m].a))%mod);    return 0;}