快速傅里叶变换模板题
算法理解请看《算法导论》第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 }
迭代版：（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 }
位逆序置换：（orz trz爷）(27~30行)
TRZ：比如我们要把一个二进制数加一，不是从最低位开始找到一个0，把它变成1，它之后的都变0嘛？那么现在我们有了i的位逆序表示，要求i+1的，不就是从最高位找0嘛，那个循环就是干这个事的。