lp3803 【模板】多项式乘法（FFT）（非递归版本）

快速傅里叶变换（FFT）的一种较高效的实现

*由于快速傅里叶变换（FFT）常常被用在通信领域，故而FFT的实现效率成为了一个需要被考察的对象。

事实上，观察FFT的实现过程，我们发现它没有非递归不可的理由。

具体来说，我们发现，如果我们想要把FFT改造为递推模式，最大的障碍就是每一次合并操作。
在快速离散傅里叶变换中，每一次合并操作都会导致两个距离为\(\frac{n}{2}\)的点的值合并在一起。因此，如果将原序列传入，那么得到的就是一个乱序的队列。
我们考虑将序列以什么顺序传入。
我们发现，每一次这样的合并操作本质上都是将两个部分「翻转」了。
以0~7为例：
0 1 2 3 4 5 6 7
0 4 2 6 1 5 3 7
将它们写成二进制：
000 001 010 011 100 101 110 111
000 100 010 110 001 101 011 111
发现变换之后每个数的二进制位都倒转了。
故而我们可以通过倒转二进制位的方式来预处理原序列，使得FFT之后的序列顺序是正常的。
但是不知道是什么原因，这个代码跑出来的速度竟然也是4500ms…这就很让人挠头了。
于是向PinkRabbit巨神请教了几个经验：
1.不用中间数组，直接在原数组上操作。
2.调整赋值操作，使得修改尽量少。
3.将蝴蝶变换预处理为一个数组。
综合几点，就成功地卡到了2500ms
当然还是有待卡常，毕竟这还是一个比较危险的时间。

#include<iostream>
#include<cstdio>
#include<cmath>
#define PI 3.141592653589793238462

struct Complex{
	double r;
	double i;
	Complex(double ri=0.0,double ii=0.0):r(ri),i(ii){}
	inline Complex operator+(const Complex &B)const{
		return (Complex){r+B.r,i+B.i}; 
	}
	inline Complex operator-(const Complex &B)const{
		return (Complex){r-B.r,i-B.i};
	}
	inline Complex operator*(const Complex &B)const{
		return (Complex){r*B.r-i*B.i,r*B.i+i*B.r};
	}
}a[1<<21|1],b[1<<21|1],s[1<<21|1];
int r[1<<21|1];
inline void FFT(Complex *A,int N,int typ){
	//ButterflyTransform
	for(int i=0;i<N;++i){
		(r[i]>i)?std::swap(A[i],A[r[i]]),0:0;
	}
	Complex bs,nw,x,y;
	for(int i=2,mid;i<=N;i<<=1){
		bs=Complex(cos(2.0*PI/i),typ*sin(2.0*PI/i));
		mid=i>>1;
		for(int j=0;j<N;j+=i){
			nw=Complex(1.0,0.0);
			for(int k=0;k<mid;++k,nw=nw*bs){
				x=A[j+k],y=(A[mid+j+k]*nw);
				A[j+k]=x+y;
				A[mid+j+k]=x-y;
			}
		}
	}
}
int n,m,cnt=1,byt=0;;
void init(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;++i){
		scanf("%lf",&a[i].r);
	}
	for(int i=0;i<=m;++i){
		scanf("%lf",&b[i].r);
	}
	while(cnt<=n+m){
		cnt<<=1;
		++byt;
	}
	for(int i=0;i<cnt;++i){
		r[i]=(r[i>>1]>>1)|((i&1)<<(byt-1));
	}
	FFT(a,cnt,1);
	FFT(b,cnt,1);
	for(int i=0;i<cnt;++i){
		a[i]=a[i]*b[i];
	}
	FFT(a,cnt,-1);
	//注意这里typ应该是-1 
	for(int i=0;i<=n+m;++i){
		printf("%d ",(int)(a[i].r/cnt+0.5));
		//注意输出应该强转格式。 
	}
	puts("");
}

int main(){
	init();
	return 0; 
}