【UOJ】【34】多项式乘法

快速傅里叶变换模板题

算法理解请看《算法导论》第30章《多项式与快速傅里叶变换》,至于证明插值唯一性什么的看不懂也没关系啦~只要明白这个过程是怎么算的就ok。

递归版:(4252ms  23468kb)

 1 //UOJ 34 递归版
 2 #include<cmath>
 3 #include<vector>
 4 #include<cstdio>
 5 #include<cstring>
 6 #include<cstdlib>
 7 #include<complex>
 8 #include<iostream>
 9 #include<algorithm>
10 #define rep(i,n) for(int i=0;i<n;++i)
11 #define F(i,j,n) for(int i=j;i<=n;++i)
12 #define D(i,j,n) for(int i=j;i>=n;--i)
13 using namespace std;
14 const int N=256256;
15 void read(int &v){
16     v=0;int sign=1; char ch=getchar();
17     while(ch<'0' || ch>'9') {if (ch=='-') sign=-1; ch=getchar();}
18     while(ch>='0'&&ch<='9'){v=v*10+ch-'0'; ch=getchar();}
19     v*=sign;
20 }
21 /****************tamplate***********************/
22 const double PI=acos(-1);
23 typedef complex<double> comp;
24 comp a[N],b[N],c[N];
25 void FFT(comp a[],int n,int type){
26     if (n==1) return;
27     int i;
28     comp y0[n>>1],y1[n>>1];
29     for(int i=0;i<n;i+=2)
30         y0[i>>1]=a[i],y1[i>>1]=a[i+1];
31     FFT(y0,n>>1,type); FFT(y1,n>>1,type);
32     comp w0( cos(type*2*PI/n) , sin(type*2*PI/n) ), w(1,0);
33     for(i=0;i<(n>>1);i++,w*=w0)
34         a[i]=y0[i]+w*y1[i],a[i+(n>>1)]=y0[i]-w*y1[i];
35 }
36 int main(){
37     int i,temp,n,m;
38     read(n); read(m);
39     int x;
40     F(i,0,n) read(x),a[i].real()=x;
41     F(i,0,m) read(x),b[i].real()=x;
42     for(temp=1;temp<=m+n;temp<<=1);
43     FFT(a,temp,1); FFT(b,temp,1);
44     rep(i,temp) c[i]=a[i]*b[i];
45     FFT(c,temp,-1);
46     for(int i=0;i<=m+n;++i)
47         printf("%d ",int(c[i].real()/temp+0.5));
48     puts("");
49     return 0;
50 }
View Code

 迭代版:(1228ms 48088kb)(空间浪费了,其实用不了这么大)

 1 //UOJ 34 迭代版
 2 #include<cmath>
 3 #include<vector>
 4 #include<cstdio>
 5 #include<cstring>
 6 #include<cstdlib>
 7 #include<complex>
 8 #include<iostream>
 9 #include<algorithm>
10 #define rep(i,n) for(int i=0;i<n;++i)
11 #define F(i,j,n) for(int i=j;i<=n;++i)
12 #define D(i,j,n) for(int i=j;i>=n;--i)
13 using namespace std;
14 const int N=1000086;
15 void read(int &v){
16     v=0;int sign=1; char ch=getchar();
17     while(ch<'0' || ch>'9') {if (ch=='-') sign=-1; ch=getchar();}
18     while(ch>='0'&&ch<='9'){v=v*10+ch-'0'; ch=getchar();}
19     v*=sign;
20 }
21 /****************tamplate***********************/
22 const double pi=acos(-1);
23 typedef complex<double> comp;
24 comp a[N],b[N],c[N];
25 int n,m;
26 void FFT(comp *a,int n,int type){
27     for(int i=1,j=0;i<n-1;i++){//j初始为i-1的位逆序置换?……每次算出i的位逆序置换j,下一次从j开始找(二进制加法)
28         //只改变1~n-2,0和n-1位置不变
29         for(int s=n;j^=s>>=1,~j&s;);//
30         if (i<j) swap(a[i],a[j]);
31     }
32     for(int m=1;m<n;m<<=1){
33         double u=pi/m*type; comp wm(cos(u),sin(u));
34         for(int k=0;k<n;k+=(m<<1) ){
35             comp w(1,0);
36             for(int j=0;j<m;j++){
37                 comp &A=a[k+j+m],&B=a[k+j],t=w*A;
38                 A=B-t; B=B+t; w=w*wm;
39             }
40         }
41     }
42     if (type==-1) rep(i,n) a[i].real()/=n;
43 }
44 
45 int main(){
46 #ifndef ONLINE_JUDGE
47     freopen("FFT.in","r",stdin);
48 #endif
49     read(n); read(m);
50     int k=0;
51     for(k=1;k<=n||k<=m;k<<=1); k<<=1;
52     int x;
53     F(i,0,n) read(x),a[i].real()=x;
54     F(i,0,m) read(x),b[i].real()=x;
55     FFT(a,k,1); FFT(b,k,1);
56     rep(i,k) c[i]=a[i]*b[i];
57     FFT(c,k,-1);
58     F(i,0,n+m) printf("%d ",int(c[i].real()+0.5));
59     return 0;
60 }
View Code

 位逆序置换:(orz trz爷)(27~30行)

TRZ:比如我们要把一个二进制数加一,不是从最低位开始找到一个0,把它变成1,它之后的都变0嘛?那么现在我们有了i的位逆序表示,要求i+1的,不就是从最高位找0嘛,那个循环就是干这个事的。

posted @ 2015-01-18 19:22  Tunix  阅读(258)  评论(0编辑  收藏  举报