首页 > 代码库 > FFT模板(From MG)

FFT模板(From MG)

 1 #include<cstdio>
 2 #include<cmath>
 3 #include<algorithm>
 4 using namespace std;
 5 struct cp{double x,y;};
 6 int n1,n2,n,m;
 7 double pi=acos(-1);
 8 cp a[500010],b[500010],cur[500010];
 9 cp operator *(cp x,cp y){return (cp){x.x*y.x-x.y*y.y,x.x*y.y+x.y*y.x};}
10 cp operator +(cp x,cp y){return (cp){x.x+y.x,x.y+y.y};}
11 cp operator -(cp x,cp y){return (cp){x.x-y.x,x.y-y.y};}
12 void fft(cp *a,int n,int fl)
13 {
14     for (int i=n>>1,j=1;j<n;j++)
15     {
16         if (i<j) swap(a[i],a[j]);
17         int k=n>>1;
18         for (;k&i;i^=k,k>>=1);
19         i^=k;
20     }
21     for (int m=2;m<=n;m<<=1)
22     {
23         cp w=(cp){cos(2*pi*fl/m),sin(2*pi*fl/m)};
24         cur[0]=(cp){1,0};
25         for (int i=1;i<m;i++) cur[i]=cur[i-1]*w;
26         for (int i=0;i<n;i+=m)
27             for (int j=i;j<i+(m>>1);j++)
28             {
29                 cp u=a[j],v=a[j+(m>>1)]*cur[j-i];
30                 a[j]=u+v;
31                 a[j+(m>>1)]=u-v;
32             }
33     }
34 }
35 int main()
36 {
37     scanf("%d%d",&n1,&n2);n1++;n2++;
38     for (int i=0;i<n1;i++) scanf("%lf",&a[i].x);
39     for (int i=0;i<n2;i++) scanf("%lf",&b[i].x);
40     n=max(n1,n2);
41     m=1;while (m<=n*2) m*=2;
42     fft(a,m,1);fft(b,m,1);
43     for (int i=0;i<=m;i++) a[i]=a[i]*b[i];
44     fft(a,m,-1);
45     for (int i=0;i<n1+n2-1;i++) printf("%d ",(int)(a[i].x/m+0.5));
46 }

 

FFT模板(From MG)