1. 程式人生 > >FFT的一個小技巧

FFT的一個小技巧

蒟蒻的FFT學習筆記IN

對於一般的,我們將兩個多項式相乘,都是開兩個虛數陣列, A , B A,B ,先將每個多項式的每個係數分別存入兩個陣列的實部,然後先將 A

A 做FFT變換,再將 B B 做FFT變換,然後兩個點值表示式相乘,然後再將其 F F T FFT
變換回來,最後答案就是陣列的實部。

這樣是做了3次 n l o g n nlogn 的FFT變換。


但是有一個小技巧,我們對於兩個多項式,其實只用開一個虛數陣列,將第一個多項式 的係數存入這個陣列的實部,將第二個多項式的係數存入這個陣列的虛部,然後先對其進行FFT變換,然後自己乘自己後再逆變換回來,虛部除以 2

n 2n 便是最後的答案。

此時,我們可以發現,只做了2次 n l o g n nlogn 的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;
}