FFT的一个小技巧

蒟蒻的FFT学习笔记IN

对于一般的,我们将两个多项式相乘,都是开两个虚数数组,A,BA,B,先将每个多项式的每个系数分别存入两个数组的实部,然后先将AA做FFT变换,再将BB做FFT变换,然后两个点值表达式相乘,然后再将其FFTFFT变换回来,最后答案就是数组的实部。

这样是做了3次nlognnlogn的FFT变换。


但是有一个小技巧,我们对于两个多项式,其实只用开一个虚数数组,将第一个多项式 的系数存入这个数组的实部,将第二个多项式的系数存入这个数组的虚部,然后先对其进行FFT变换,然后自己乘自己后再逆变换回来,虚部除以2n2n便是最后的答案。

此时,我们可以发现,只做了2次nlognnlogn的FFT变换。

目前原理待博主研究后再写上。

上份代码:
模板

Luogu FFT 模板题
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define db double
using namespace std;
const int M=3e6+10;
const db pi=acos(-1);
struct complex{
	db x,y;
	complex(){}
	complex(db a,db b):x(a),y(b){}
}A[M];
complex operator +(const complex &a,complex &b){return complex(a.x+b.x,a.y+b.y);}
complex operator -(const complex &a,complex &b){return complex(a.x-b.x,a.y-b.y);}
complex operator *(const complex &a,complex &b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
int R[M];
void DFT(complex *a,int n,int f){
	for(int i=0;i<=n;i++)if(i<R[i])swap(a[i],a[R[i]]);
	for(int i=2;i<=n;i<<=1){
		int now=i>>1;
		complex wn=complex(cos(2*pi/i),f*sin(2*pi/i));
		for(int j=0;j<n;j+=i){
			complex w=complex(1,0),x,y;
			for(int k=j;k<j+now;k++,w=w*wn){
				x=a[k],y=w*a[k+now];
				a[k]=x+y;a[k+now]=x-y;
			}
		}
	}
	if(f==-1)for(int i=0;i<=n;i++)a[i].y/=(2*n);
}
void FFT(complex *a,int n,int m){
	int K,lg;
	for(K=1,lg=0;K<=n+m;K<<=1,++lg);--lg;
	for(int i=0;i<=K;i++)R[i]=(R[i>>1]>>1)|((i&1)<<lg);
	DFT(a,K,1);
	for(int i=0;i<=K;i++)a[i]=a[i]*a[i];
	DFT(a,K,-1);
}
int n,m,x;
int main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++)scanf("%d",&x),A[i].x=x;
	for(int i=0;i<=m;i++)scanf("%d",&x),A[i].y=x;
	FFT(A,n,m);
	for(int i=0;i<=n+m;i++){
		printf("%d ",int(A[i].y+0.5));
	}
	return 0;
}
posted @ 2018-12-01 16:49  VictoryCzt  阅读(123)  评论(0编辑  收藏  举报