FFT 【JSOI2012】bzoj4332 分零食

出处：http://blog.csdn.net/todobe
题目大意：

有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.

注意：

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

#include <cstdio>
#include <cmath>
#include <algorithm>
#include <iostream>
#define N 1200000
using 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_Transform
void 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;
}