FFT

Posted on |

title: 快速傅里叶变换
date: 2020/2/24
tag: 算法
mathjax: true

P1919 A×B Problem
快速傅里叶变换FFT 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
#include<bits/stdc++.h>
using namespace std;
const int maxn=(1<<21)+1;
const double PI=acos(-1);
struct Complex{
	double r,i;
	Complex(){r=0,i=0;}
	Complex(double real,double imag):r(real),i(imag){}
}F[maxn],G[maxn];
Complex operator +(Complex a,Complex b){return Complex(a.r+b.r,a.i+b.i);}
Complex operator -(Complex a,Complex b){return Complex(a.r-b.r,a.i-b.i);}
Complex operator *(Complex a,Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
int rev[maxn],len,lim=1,totF=0,totG=0,ans[maxn];
char n[maxn],m[maxn];
void FFT(Complex *a,double opt){
	for(int i=0;i<lim;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
	for(int dep=1;dep<=log2(lim);++dep){
		int m=(1<<dep);
		Complex wn=Complex(cos(2.0*PI/m),opt*sin(2.0*PI/m));
		for(int k=0;k<lim;k+=m){
			Complex w=Complex(1,0);
			for(int j=0;j<(m>>1);++j){
				Complex t=w*a[k+j+(m>>1)];
				Complex u=a[k+j];
				a[k+j]=u+t;
				a[k+j+(m>>1)]=u-t;
				w=w*wn;
			}
		}
	}
	if(opt==-1) for(int i=0;i<lim;++i) a[i].r/=(double)lim;
}
int main(){
	scanf("%s%s",&n,&m);
	int nlen=strlen(n),mlen=strlen(m);
	for(int i=nlen-1;i>=0;--i) F[totF++].r=(n[i]^48);
	for(int i=mlen-1;i>=0;--i) G[totG++].r=(m[i]^48);
	mlen--;nlen--;
	while(lim<=nlen+mlen) lim<<=1,++len;
	for(int i=0;i<lim;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	FFT(F,1);FFT(G,1);
	for(int i=0;i<=lim;++i) F[i]=F[i]*G[i];
	FFT(F,-1);
	int tot=0;
	for(int i=0;i<=lim;++i){
		ans[i]+=(int)(F[i].r+0.5);
		if(ans[i]>=10) ans[i+1]+=ans[i]/10,ans[i]%=10,lim+=(i==lim);
	}
	while((!ans[lim])&&lim>=1) --lim;
	++lim;
	while(lim--) cout<<ans[lim];
	return 0;
}
0%