FFT快速傅里叶
来源:互联网 发布:刚柔流空手道 知乎 编辑:程序博客网 时间:2024/06/16 02:12
传送门
表示只会抄板子,看了个半懂。。
我因为重载运算符出了点儿问题,调了好久好久……
贴代码:
#include<cstdio>#include<cstring>#include<iostream>#include<cmath>#include<algorithm>#include<cstdlib>using namespace std;inline int read(){ int x=0;char ch=' ';int f=1; while(ch!='-'&&(ch<'0'||ch>'9'))ch=getchar(); if(ch=='-')f=-1,ch=getchar(); while(ch>='0'&&ch<='9')x=(x<<3)+(x<<1)+(ch^48),ch=getchar(); return x*f;}const int N=2e5+5;const double pi=acos(-1);struct cp{ double r,i; cp(){} cp(double _r,double _i):r(_r),i(_i){} friend cp operator + (const cp& a,const cp& b){return cp(a.r+b.r,a.i+b.i);} friend cp operator - (const cp& a,const cp& b){return cp(a.r-b.r,a.i-b.i);} friend cp operator * (const cp& a,const cp& b){return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}}a[N],b[N];int n,L,R[N],c[N];char ch[N];inline void FFT(cp *a,int f){ for(int i=0;i<n;++i)if(i<R[i])swap(a[i],a[R[i]]); for(int i=1;i<n;i<<=1){ cp wn(cos(pi/i),f*sin(pi/i)); for(int j=0;j<n;j+=(i<<1)){ cp w(1,0); for(int k=0;k<i;++k,w=w*wn){ cp x=a[j+k],y=w*a[j+k+i]; a[j+k]=x+y;a[j+k+i]=x-y; } } } if(f==-1)for(int i=0;i<n;++i)a[i].r/=n;}int main(){ n=read()-1; scanf("%s",ch);for(int i=0;i<=n;++i)a[i].r=ch[n-i]-'0'; scanf("%s",ch);for(int i=0;i<=n;++i)b[i].r=ch[n-i]-'0'; int m=(n<<1);for(n=1;n<=m;n<<=1)++L; for(int i=0;i<n;++i)R[i]=(R[i>>1]>>1)|((i&1)<<(L-1)); FFT(a,1);FFT(b,1); for(int i=0;i<=n;++i)a[i]=a[i]*b[i]; FFT(a,-1); for(int i=0;i<=m;++i)c[i]=(int)(a[i].r+0.5); for(int i=0;i<=m;++i)if(c[i]>=10){c[i+1]+=c[i]/10;c[i]%=10;if(i==m)++m;} while(m){if(c[m])break;else m--;} while(~m)putchar(c[m--]+'0'); return 0;}
阅读全文