bzoj 4332: JSOI2012 分零食 快速傅立叶变换

题目:

Description
同学们依次排成了一列，其中有A位小朋友，有三个共同的欢乐系数O，S和U。如果有一位小朋友得到了x个糖果，那么她的欢乐程度就是(f(x)=O*x^2+S*x+U)
现在校长开始分糖果了，一共有M个糖果。有些小朋友可能得不到糖果，对于那些得不到糖果的小朋友来说，欢乐程度就是1。如果一位小朋友得不到糖果，那么在她身后的小朋友们也都得不到糖果。（即这一列得不到糖果的小朋友一定是最后的连续若干位）
所有分糖果的方案都是等概率的。现在问题是：期望情况下，所有小朋友的欢乐程度的乘积是多少？呆呆同学很快就有了一个思路，只要知道总的方案个数T和所有方案下欢乐程度乘积的总和S，就可以得到答案Ans=S/T。现在他已经求出来了T的答案，但是S怎么求呢？他就不知道了。你能告诉他么？
因为答案很大，你只需要告诉他S对P取模后的结果。

题解:

首先这道题我们可以考虑枚举一下有多少人得到了零食
设(g[i][j])表示(i)个人里分下去了(j)个零食得到的值,n为人数,m为零食数
这样我们有(ans = sum_{i=1}^ng[i][m])
(g[i][j])的递推我们有

[g[i][j] = sum_{k=1}^{j-1}g[i-1][j-k]*F(k) ]

其中(F(k))表示将(k)个零食分给一个人得到的权
然后我们惊奇地发现后面的式子是一个卷积的形式
所以我们可以得到(g_i = g_{i-1}*F)
由于卷积满足结合律,所以我们有(g_i = g_0*F^i)
这样的话我们能够完成(nlogn)的单点求值,但是我们要求的是(g_i)的一个和
我们把答案(ans)拓展为一个多项式(f(x)),答案储存在第f(n)的第(m)位
则有(f_n = sum_{i=1}^{n}g_i)
然后我们发现:

[f_n = f_{frac{n}{2}} + sum_{i=frac{n}{2}+1}^{n}g_i ]

[f_n = f_{frac{n}{2}} + sum_{i=1}^{frac{n}{2}}g_{frac{n}{2}+i} ]

在继续推导之前首先我们需要证明: (g_{i+j} = g_i*g_j)
由(g_i = g_0*F^i)可得:(g_i*g_j = g_0*g_0*F^i*F^j)
因为:(g_0*g_0 = g_0)所以有(g_0*g_0*F^i*F^j = g_0*F^i*F^j = g_{i+j})得证

所以继续上式的推导我们有

[f_n = f_{frac{n}{2}} + sum_{i=1}^{frac{n}{2}}g_{frac{n}{2}}*g_i ]

我们将卷积的形式再拆解开来:

[f_n = f_{frac{n}{2}} + sum_{i=1}^{frac{n}{2}}sum_{j=1}^{m-1}g_{frac{n}{2},m-j}*g_{i,j} ]

[f_n = f_{frac{n}{2}} + sum_{j=1}^{m-1}g_{frac{n}{2},m-j}*sum_{i=1}^{frac{n}{2}}g_{i,j} ]

[f_n = f_{frac{n}{2}} + sum_{j=1}^{m-1}g_{frac{n}{2},m-j}*f_{frac{n}{2},j} ]

[f_n = f_{frac{n}{2}} + g_{frac{n}{2}}*f_{frac{n}{2}} ]

我们又知道:

[g_n = g_{frac{n}{2}}*g_{frac{n}{2}} ]

所以迭代倍增即可求解
复杂度(O(nlog^2n))

#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
inline void read(int &x){
	x=0;char ch;bool flag = false;
	while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
	while(x=10*x+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
const int maxn = 41000;
const double pi = acos(-1);
int mod;
struct complex{
	long double x,y;
	complex(){x=y=0;}
	complex(long double a,long double b){x=a;y=b;}
	complex operator + (const complex &r){return complex(x+r.x,y+r.y);}
	complex operator - (const complex &r){return complex(x-r.x,y-r.y);}
	complex operator * (const complex &r){return complex(x*r.x-y*r.y,x*r.y+y*r.x);}
	complex operator / (const long double &r){return complex(x/r,y/r);}
};
inline void FFT(complex *x,int n,int p){
	for(int i=0,t=0;i<n;++i){
		if(i > t) swap(x[i],x[t]);
		for(int j=n>>1;(t^=j) < j;j>>=1);
	}
	for(int m=2;m<=n;m<<=1){
		int k = m>>1;
		complex wn(cos(p*2*pi/m),sin(p*2*pi/m));
		for(int i=0;i<n;i+=m){
			complex w(1,0),u;
			for(int j=0;j<k;++j,w=w*wn){
				u = x[i+j+k]*w;
				x[i+j+k] = x[i+j] - u;
				x[i+j] = x[i+j] + u;
			}
		}
	}
	if(p == -1 ) for(int i=0;i<n;++i) x[i] = x[i]/n;
}
complex ca[maxn],cb[maxn],cc[maxn];int len,m;
inline int mul(int *a,int *b,int *c){
	for(int i=0;i<len;++i){
		ca[i] = complex((long double)a[i],0);
		cb[i] = complex((long double)b[i],0);
	}
	FFT(ca,len,1);FFT(cb,len,1);
	for(int i=0;i<len;++i) cc[i] = ca[i]*cb[i];
	FFT(cc,len,-1);
	for(int i=0;i<=m;++i){
		c[i] = ((int)floor(cc[i].x + 0.5)) % mod;
	}
}
int f[maxn],g[maxn],arr[maxn],tmp[maxn];
inline void qpow(int k){
	if(k == 1){
		for(int i=0;i<=m;++i) f[i] = g[i] = arr[i];
		return ;
	}qpow(k>>1);
	mul(f,g,tmp);mul(g,g,g);
	for(int i=0;i<=m;++i){
		f[i] += tmp[i];
		if(f[i] >= mod) f[i] -= mod;
	}
	if(k&1){
		mul(g,arr,g);
		for(int i=0;i<=m;++i){
			f[i] += g[i];
			if(f[i] >= mod) f[i] -= mod;
		}
	}
}
int main(){
	read(m);read(mod);
	for(len = 1;(len) <= (m<<1);len<<=1);
	int n,a,b,c;read(n);read(a);read(b);read(c);
	a %= mod;b %= mod;c %= mod;
	for(int i=1;i<=m;++i) arr[i] = ((a*i*i % mod) + (b*i % mod) + c) % mod;
	qpow(n);printf("%d
",f[m]);
	getchar();getchar();
	return 0;
}