bzoj 2194: 快速傅立叶之二 FFT

题目就是题解系列

既然题目告诉了我们用FFT,我们就用去想FFT.
我们看这个式子,好像和平时的式子差不多
那么我们直接裸上FFT就好了
我们知道,FFT处理的是卷积的问题。
我们在多项式乘法的时候,得到的下标为\(k\)上面
实际上是\(\sum(i*j)*[i+j=k]\)
所以我们发现,使用FFT处理问题的时候,只会统计所有和为k的项,而原题中的

\[C_k = \sum_i(a_i*b_{i-k}) \]

\(i + (i-k)\)明显不是定值,所以不能用FFT直接统计
那么我们把a数组倒过来
式子就变成了

\[C_k = \sum_i(a_{n-i-1}*b_{i-k}) \]

所以我们发现\((n-i-1) + (i-k) = n-k-1\)
是一个只与n,k有关的定值
所以我们使用FFT加速即可

#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;
}
inline int cat_max(const int &a,const int &b){return a>b ? a:b;}
inline int cat_min(const int &a,const int &b){return a<b ? a:b;}
const int maxn = 410010;
const double pi = acos(-1);
struct complex{
	double x,y;
	complex(){}
	complex(double a,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 double &r){return complex(x/r,y/r);}
};
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){
		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;
			int k = m>>1;
			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 a[maxn],b[maxn],c[maxn];
int main(){
	int n;read(n);
	for(int i=0,x;i<n;++i){
		read(x);
		a[n-i-1].x = x;
		read(x);
		b[i].x = x;
	}
	int len;
	for(int i=1;(i>>2)<n;i<<=1) len = i;
	FFT(a,len,1);FFT(b,len,1);
	for(int i=0;i<len;++i) c[i] = a[i]*b[i];
	FFT(c,len,-1);
	for(int i = n-1;i>=0;--i){
		printf("%lld\n",(ll)(c[i].x+0.5));
	}
	getchar();getchar();
	return 0;
}
  
posted @ 2017-01-28 16:32  Sky_miner  阅读(411)  评论(0编辑  收藏  举报