BZOJ 3527 [ZJOI 3014] 力 （FFT）

题目链接：BZOJ 3527

因为这道题在BZOJ上没有贴，所以我还是附上数据吧。

input：

5
4006373.885184
15375036.435759
1717456.469144
8514941.004912
1410681.345880

output：
-16838672.693
3439.793
7509018.566
4595686.886
10903040.872

做完这道题的第一感受是：我也是会写FFT的人了！！！（虽说只会写模板= =）。这道题还算是比较裸的FFT，我也是来学习写FFT的。题解什么的，由于涉及到许多高数知识，博主蒟蒻知道的少，所以我还是简单说说吧。。。

简单说一下FFT吧。之前一直不知道FFT怎么求卷积，现在才知道。对于函数y[ ],f[ ],g[ ],满足y[n]=f[n]*g[n]，快速求得y[n]，设函数F[f]=DET(f[n])，G[f]=DET(g[n])，那么函数y[n]可以表示为：y[n]=IDFT(DFT(f[n])*DFT(g[n]))。而对于卷积(f*g)[m]=segma(f[n]*g[m-n])就可以快速求出，简化乘法和加法运算。因为FFT是基于复数运算的，对于一次复数的乘法运算相当于四次实数乘法运算和二次实数加法，对于一次复数加法相当于二次实数加法，这样就简化了运算。

代码中有一个小技巧，即观察计算DFT和IDFT的式子，我们发现两者只相差一个符号和一个常数。代码中的v参数就是来改变符号的，至于对于要除掉的那个常数，可以在最后输出答案的时候除掉，这样就只用写一个函数。

#include<cstdio>
#include<cstring>
#include<iostream>
#include<cmath>
using namespace std;
int N,M;
double pi=acos(-1);

struct complex{
	double r,i;
	complex(double r=0,double i=0):r(r),i(i){}
	complex operator + (const complex &A)const{
		return complex(r+A.r,i+A.i);
	}
	complex operator - (const complex &A)const{
		return complex(r-A.r,i-A.i);
	}
	complex operator * (const complex &A)const{
		return complex(r*A.r-i*A.i,r*A.i+i*A.r);
	}
};
complex A[600000+10],B[600000+10];

void DFT(complex A[],int v=0){//v用于辨别当前是进行DFT还是IDFT操作
	for(int i=0,j=0;i<M;i++){
		if(i<j)swap(A[i],A[j]);
		for(int k=M>>1;(j^=k)<k;k>>=1)
			;
	}
	for(int i=1;i<M;i<<=1){
		complex w(1.,0.),wn(cos(pi/i),sin(pi/i));
		if(v)wn.i=-wn.i;
		for(int j=0;j<i;j++){
			for(int k=j;k<M;k+=i*2){
				complex u=A[k],v=A[k+i]*w;
				A[k]=u+v,A[k+i]=u-v;
			}
			w=w*wn;
		}
	}
}

void FFT(){
	for(M=1;M<2*N-1;M<<=1)
		;
	M<<=1;
	for(int i=N;i<M;i++)A[i]=0;
	for(int i=2*N-1;i<M;i++)B[i]=0;
	DFT(A);
	DFT(B);
	for(int i=0;i<M;i++)A[i]=A[i]*B[i];
	DFT(A,1);
	for(int i=0;i<N;i++){
		double t=A[i+N-1].r/M;//因为最后一步是IDFT操作所以要除以一个常数
		printf("%lf\n",t);
	}
}

int main(){
	scanf("%d",&N);
	for(int i=0;i<N;i++)scanf("%lf",&A[i]);
	for(int i=N-2;i>=0;i--){
		double t=i+1;
		B[N-i-2]=-1.0/t/t;
	}
	B[N-1]=0;
	for(int i=1;i<N;i++){
		double t=i;
		B[i+N-1]=1.0/t/t;
	}
	FFT();
	return 0;
}