码迷,mamicode.com
首页 > 其他好文 > 详细

FFT多项式乘法加速

时间:2015-05-30 21:11:04      阅读:144      评论:0      收藏:0      [点我收藏+]

标签:

FFT基本操作。。。讲解请自己看大学信号转置系列。。。

15-5-30更新:改成结构体的,跪烂王学长啊啊啊啊机智的static。。。

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cmath>
 4 #include<algorithm>
 5 #include<queue>
 6 #include<cstring>
 7 #define PAU putchar(‘ ‘)
 8 #define ENT putchar(‘\n‘)
 9 #pragma comment(linker,"/STACK:10240000,10240000")
10 using namespace std;
11 const double PI=acos(-1.0);
12 const int maxn=800000+10;
13 struct FFT{
14     struct cox{
15         double r,i;cox(double _r = 0.0,double _i = 0.0){r=_r;i=_i;}
16         cox operator +(const cox &b){return cox(r+b.r,i+b.i);}
17         cox operator -(const cox &b){return cox(r-b.r,i-b.i);}
18         cox operator *(const cox &b){return cox(r*b.r-i*b.i,r*b.i+i*b.r);}
19     }f[maxn];int len,lenx;
20     void init(int*s,int L,int Len){
21         len=L;lenx=Len;
22         for(int i=0;i<L;i++) f[i]=cox(s[L-1-i],0);
23         for(int i=L;i<Len;i++) f[i]=cox(0,0);return;
24     }
25     void change(){
26         for(int i=1,j=lenx>>1;i<lenx-1;i++){
27             if(i<j)swap(f[i],f[j]);
28             int k=lenx>>1;
29             while(j>=k) j-=k,k>>=1;
30             if(j<k) j+=k;
31         } return;
32     }
33     void cal(int tp){
34         change();
35         for(int h=2;h<=lenx;h<<=1){
36             double tr=-tp*2*PI/h;
37             cox wn(cos(tr),sin(tr));
38             for(int j=0;j<lenx;j+=h){
39                 cox w(1,0);
40                 for(int k=j;k<j+(h>>1);k++){
41                     cox u=f[k],t=w*f[k+(h>>1)];
42                     f[k]=u+t;f[k+(h>>1)]=u-t;w=w*wn;
43                 }
44             }
45         } if(tp==-1) for(int i=0;i<lenx;i++) f[i].r/=lenx;return;
46     }
47 };
48 void mul(int*s1,int*s2,int L1,int L2,int&L,int*ans){
49     L=1;while(L<L1<<1||L<L2<<1) L<<=1;
50     static FFT a,b;a.init(s1,L1,L);b.init(s2,L2,L);a.cal(1);b.cal(1);
51     for(int i=0;i<L;i++) a.f[i]=a.f[i]*b.f[i];a.cal(-1);
52     for(int i=0;i<L;i++) ans[i]=(int){a.f[i].r+0.5};return;
53 }
54 int s1[maxn>>1],s2[maxn>>1],ans[maxn],L1,L2,L;
55 void init(){
56     char ch;int tot;
57     do ch=getchar(); while(!isdigit(ch));tot=0;
58     while(isdigit(ch)){s1[tot++]=ch-0;ch=getchar();}L1=tot;
59     do ch=getchar(); while(!isdigit(ch));tot=0;
60     while(isdigit(ch)){s2[tot++]=ch-0;ch=getchar();}L2=tot;
61     mul(s1,s2,L1,L2,L,ans);
62     return;
63 }
64 void work(){
65     for(int i=0;i<L;i++){
66         ans[i+1]+=ans[i]/10;ans[i]%=10;
67     } L=L1+L2-1;return;
68 }
69 void print(){
70     while(ans[L]<=0&&L>0) L--;
71     for(int i=L;i>=0;i--) putchar(ans[i]+0);
72     return;
73 }
74 int main(){init();work();print();return 0;}

原来写的丑哭了TAT

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cmath>
 4 #include<algorithm>
 5 #include<queue>
 6 #include<cstring>
 7 #define PAU putchar(‘ ‘)
 8 #define ENT putchar(‘\n‘)
 9 using namespace std;
10 const double PI=acos(-1.0);
11 struct complex{
12     double r,i;
13     complex(double _r = 0.0,double _i = 0.0){r=_r;i=_i;}
14     complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
15     complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
16     complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
17 };
18 void change(complex y[],int len){
19     int i,j,k;
20     for(i=1,j=len/2;i<len-1;i++){
21         if(i<j)swap(y[i],y[j]);
22         k=len/2;
23         while(j>=k) j-=k,k/=2;
24         if(j<k) j+=k;
25     } return;
26 }
27 void fft(complex y[],int len,int on){
28     change(y,len);
29     for(int h=2;h<=len;h<<=1){
30         complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
31         for(int j=0;j<len;j+=h){
32             complex w(1,0);
33             for(int k=j;k<j+h/2;k++){
34                 complex u=y[k],t=w*y[k+h/2];
35                 y[k]=u+t;
36                 y[k+h/2]=u-t;
37                 w=w*wn;
38             }
39         }
40     }
41     if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
42     return;
43 }
44 const int MAXN=2000+10;
45 complex x1[MAXN],x2[MAXN];
46 char str1[MAXN>>1],str2[MAXN>>1];
47 int sum[MAXN];
48 inline int read(){
49     int x=0,sig=1;char ch=getchar();
50     while(!isdigit(ch)){if(ch==-)sig=-1;ch=getchar();}
51     while(isdigit(ch))x=10*x+ch-0,ch=getchar();
52     return x*=sig;
53 }
54 inline void write(int x){
55     if(x==0){putchar(0);return;}if(x<0)putchar(-),x=-x;
56     int len=0,buf[15];while(x)buf[len++]=x%10,x/=10;
57     for(int i=len-1;i>=0;i--)putchar(buf[i]+0);return;
58 }
59 int len,len1,len2;
60 void init(){
61     scanf("%s%s",str1,str2);
62     len1=strlen(str1),len2=strlen(str2),len=1;
63     return;
64 }
65 void work(){
66     while(len<len1<<1||len<len2<<1) len<<=1;
67     for(int i=0;i<len1;i++) x1[i]=complex(str1[len1-1-i]-0,0);
68     for(int i=len1;i<len;i++) x1[i]=complex(0,0);
69     for(int i=0;i<len2;i++) x2[i]=complex(str2[len2-1-i]-0,0);
70     for(int i=len2;i<len;i++) x2[i]=complex(0,0);
71     fft(x1,len,1);fft(x2,len,1);
72     for(int i=0;i<len;i++) x1[i]=x1[i]*x2[i];
73     fft(x1,len,-1);
74     for(int i=0;i<len;i++) sum[i]=(int)(x1[i].r+0.5);
75     for(int i=0;i<len;i++){
76         sum[i+1]+=sum[i]/10;
77         sum[i]%=10;
78     } len=len1+len2-1;return;
79 }
80 void print(){
81     while(sum[len]<=0&&len>0) len--;
82     for(int i=len;i>=0;i--) putchar(sum[i]+0);
83     return;
84 }
85 int main(){init();work();print();return 0;}

 

FFT多项式乘法加速

标签:

原文地址:http://www.cnblogs.com/chxer/p/4507345.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!