FFT模板

来源：互联网发布：principle软件编辑：程序博客网时间：2024/06/02 06:49

FFT模板确实不太好写，虽然看思想比较快，但是一看代码依旧一脸茫然

FFT主要是利用了DFT二进制的性质来进行奇技淫巧的转换

还是把写的模板贴下(抄的大神的2333333)

#include <stdio.h>#include <string.h>#include <math.h>#include <algorithm>#define pi acos(-1.0)using namespace std;const int maxn=400007;struct complex{double r,i;complex(double real=0.0,double image=0.0){r=real,i=image;}complex operator + (complex a1){return complex(r+a1.r,i+a1.i);}complex operator - (complex a1){return complex(r-a1.r,i-a1.i);}complex operator * (complex a1){return complex(r*a1.r-i*a1.i,r*a1.i+i*a1.r);}}x1[maxn],x2[maxn]; char a[maxn/2],b[maxn/2];int sum[maxn];int vis[maxn];void brc(complex *a,int l){memset(vis,0,sizeof(vis));for(int i=1;i<l-1;i++){int x=i,y=0;if(vis[x]) continue;int m=(int)log2(l)+0.1;while(m--){y<<=1;y|=(x&1);x>>=1;}vis[i]=vis[y]=1;swap(a[i],a[y]);}}void fft(complex *a,int l,double on){complex u,t;brc(a,l);for(int h=2;h<=l;h<<=1){complex wn(cos(on*2*pi/h),sin(on*2*pi/h));for(int j=0;j<l;j+=h){complex w(1,0);for(int k=j;k<j+h/2;k++){u=a[k];t=w*a[k+h/2];a[k]=u+t;a[k+h/2]=u-t;w=w*wn;}}}if(on==-1) for(int i=0;i<l;i++) a[i].r/=l;}int main(){scanf("%s%s",a,b);int l1=strlen(a),l2=strlen(b);int l=1;while(l<l1*2||l<l2*2) l<<=1;for(int i=0;i<l1;i++){x1[i].r=a[l1-i-1]-'0';x1[i].i=0.0;}for(int i=l1;i<l;i++)x1[i].r=x1[i].i=0.0;for(int i=0;i<l2;i++){x2[i].r=b[l2-i-1]-'0';x2[i].i=0.0;}for(int i=l2;i<l;i++) x2[i].r=x2[i].i=0.0;fft(x1,l,1);fft(x2,l,1);for(int i=0;i<l;i++) x1[i]=x1[i]*x2[i];fft(x1,l,-1);for(int i=0;i<l;i++) sum[i]=x1[i].r+0.5;for(int i=0;i<l;i++){sum[i+1]+=sum[i]/10;sum[i]%=10;}l=l1+l2-1;while(sum[l]<=0&&l) l--;for(int i=l;i>=0;i--) printf("%d",sum[i]);printf("\n");return 0;}

0 0