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.

注意：

亘古不变的精度问题！！！！

在乘的时候不要连乘啊，要不然会变成循环卷积啊什么的云云，会错，就是说每次在做乘法的时候先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
;
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
,p
,f
,tmp
;
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;
}