FFT
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<ctime> 5 #include<iostream> 6 #include<algorithm> 7 #include<queue> 8 #include<stack> 9 #include<set> 10 #define esp (1e-8) 11 #define maxint (2147483647) 12 #define pi (3.141592653589793) 13 #define l(a) ((a)<<1) 14 #define r(a) ((a)<<1|1) 15 #define b(a) (1<<(a)) 16 #define rep(i,a,b) for(int i=a;i<=(b);i++) 17 #define clr(a) memset(a,0,sizeof(a)) 18 typedef long long ll; 19 using namespace std; 20 int readint(){ 21 int t=0,f=1;char c=getchar(); 22 while(!isdigit(c)){ 23 if(c=='-') f=-1; 24 c=getchar(); 25 } 26 while(isdigit(c)){ 27 t=(t<<3)+(t<<1)+c-'0'; 28 c=getchar(); 29 } 30 return t*f; 31 } 32 ll readll(){ 33 ll t=0ll,f=1ll;char c=getchar(); 34 while(!isdigit(c)){ 35 if(c=='-') f=-1ll; 36 c=getchar(); 37 } 38 while(isdigit(c)){ 39 t=(t<<3ll)+(t<<1ll)+ll(c-'0'); 40 c=getchar(); 41 } 42 return t*f; 43 } 44 const int maxk=100,maxn=500009; 45 int n,m,cnt=0; 46 struct complex{ 47 double r,i; 48 inline complex(double R,double I){ 49 r=R;i=I; 50 } 51 inline complex(){} 52 inline void out(){ 53 if(abs(i)>esp) printf("%.2lf %.2lfi\n",r,i); 54 else printf("%.2lf\n",r); 55 } 56 inline complex operator +(const complex A)const{ 57 return complex(r+A.r,i+A.i); 58 } 59 inline complex operator -(const complex A)const{ 60 return complex(r-A.r,i-A.i); 61 } 62 inline complex operator *(const complex A)const{ 63 return complex(r*A.r-i*A.i,i*A.r+r*A.i); 64 } 65 }pool[maxk][maxn],*S[maxk]; 66 complex*newarray(){ 67 return S[--cnt]; 68 } 69 int rev(int t){//将t化为2^k 70 int T=0; 71 for(T=0;b(T)<t;T++); 72 return b(T); 73 } 74 complex*FFT(complex a[],int t,int opt){//t为2的整数次幂 opt==1->FFT opt==-1->逆FFT 75 if(t==1) return a; 76 complex w=complex(1,0),w1=complex(cos(2*pi/t),sin(2*pi/t*opt));//只需在sin中乘上opt 77 complex*y=newarray(),*a0=newarray(),*a1=newarray(); 78 rep(i,0,t-1){//奇偶分组 79 if(i&1) a1[i>>1]=a[i]; 80 else a0[i>>1]=a[i]; 81 } 82 complex*y0=FFT(a0,t>>1,opt),*y1=FFT(a1,t>>1,opt);//递归 83 rep(i,0,(t>>1)-1){//利用t>>1的情况计算n的情况，注意公式的推导:A(x）=A0(x^2)+xA1(x^2) 84 y[i]=y0[i]+w*y1[i]; 85 y[i+(t>>1)]=y0[i]-w*y1[i]; 86 w=w*w1; 87 } 88 S[cnt++]=a0;S[cnt++]=a1;if(y0!=a0) S[cnt++]=y0;if(y1!=a1) S[cnt++]=y1; 89 return y; 90 } 91 complex*DFT(complex a[],int t,int opt){//t为2的整数次幂 opt==1->FFT opt==-1->逆FFT 92 t=rev(t); 93 complex*ans=FFT(a,t,opt); 94 if(opt==-1) rep(i,0,t-1){ 95 (ans+i)->r/=t;(ans+i)->i/=t; 96 } 97 return ans; 98 } 99 void init(){ 100 rep(i,0,maxk-1) S[cnt++]=pool[i]; 101 } 102 int main(){ 103 //freopen("#intput.txt","r",stdin); 104 //freopen("#output.txt","w",stdout); 105 init(); 106 n=readint();m=readint();n++;m++; 107 complex*a=newarray(),*b=newarray(); 108 rep(i,0,n-1) a[i].r=readint(); 109 rep(i,0,m-1) b[i].r=readint(); 110 int N=max(n,m);N=rev(N)<<1; 111 complex*A=DFT(a,N,1),*B=DFT(b,N,1),*Ans; 112 rep(i,0,N-1) A[i]=A[i]*B[i]; 113 Ans=DFT(A,N,-1);N--; 114 rep(i,0,n+m-2){ 115 printf("%.0lf",(Ans+i)->r+esp); 116 putchar(i==N?'\n':' '); 117 } 118 //fclose(stdin); 119 //fclose(stdout); 120 return 0; 121 }
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<ctime> 5 #include<iostream> 6 #include<algorithm> 7 #include<queue> 8 #include<set> 9 #define pi (3.14159265358979323846) 10 #define maxint (2147483647) 11 #define l(a) ((a)<<1) 12 #define r(a) ((a)<<1|1) 13 #define b(a) (1<<(a)) 14 #define esp (1e-4) 15 #define rep(i,a,b) for(int i=a;i<=(b);i++) 16 #define clr(a) memset(a,0,sizeof(a)) 17 typedef long long ll; 18 using namespace std; 19 int readint(){ 20 int t=0,f=1;char c=getchar(); 21 while(!isdigit(c)){ 22 if(c=='-') f=-1; 23 c=getchar(); 24 } 25 while(isdigit(c)){ 26 t=(t<<3)+(t<<1)+c-'0'; 27 c=getchar(); 28 } 29 return t*f; 30 } 31 void readstr(char s[]){ 32 int p=0;char c=getchar(); 33 while(!isdigit(c)) c=getchar(); 34 while(isdigit(c)){ 35 s[p++]=c;c=getchar(); 36 } 37 } 38 const int maxn=250000; 39 struct complex{ 40 double r,i; 41 inline complex(double R,double I){ 42 r=R;i=I; 43 } 44 inline complex(){ 45 r=i=0; 46 } 47 inline complex operator+(complex A)const{ 48 return complex(A.r+r,A.i+i); 49 } 50 inline complex operator-(complex A)const{ 51 return complex(r-A.r,i-A.i); 52 } 53 inline complex operator*(complex A)const{ 54 return complex(A.r*r-A.i*i,A.r*i+A.i*r); 55 } 56 inline void out(){ 57 if(abs(i)>esp) printf("%.2lf+%.2lfi\n",r,i); 58 else printf("%.2lf\n",r+esp); 59 } 60 }; 61 int n,n1,n2,ans[maxn]; 62 complex A[maxn],B[maxn],_A[maxn],_B[maxn]; 63 char a[maxn],b[maxn]; 64 int init(int n){ 65 int tmp=1;while(tmp<n) tmp<<=1; 66 return tmp; 67 } 68 int rev(int t,int l){ 69 int _t=0; 70 for(int i=0;b(i)<l;i++){ 71 _t<<=1;if(t&b(i)) _t|=1; 72 } 73 return _t; 74 } 75 complex*FFT(complex X[],complex Y[],int t,int opt){ 76 rep(i,0,t-1) Y[rev(i,t)]=X[i]; 77 for(int s=2;s<=t;s<<=1){ 78 complex w1=complex(cos(2*pi/s),sin(2*pi*opt/s)); 79 for(int k=0;k<t;k+=s){ 80 complex w=complex(1,0); 81 rep(i,k,k+(s>>1)-1){ 82 complex t1=Y[i],t2=w*Y[i+(s>>1)]; 83 Y[i]=t1+t2;Y[i+(s>>1)]=t1-t2; 84 w=w*w1; 85 } 86 } 87 } 88 if(opt==-1) rep(i,0,t-1) Y[i].r/=t,Y[i].i/=t; 89 return Y; 90 } 91 int main(){ 92 //freopen("#input.txt","r",stdin); 93 //freopen("#output.txt","w",stdout); 94 n=readint(); 95 readstr(a);readstr(b); 96 while(a[n1]=='0') n1++;while(b[n2]=='0') n2++; 97 rep(i,n1,n-1) A[i-n1].r=a[i]-'0';rep(i,n2,n-1) B[i-n2].r=b[i]-'0'; 98 n1=n-n1;n2=n-n2; 99 int N=init(n); 100 FFT(A,_A,2*N,1);FFT(B,_B,2*N,1); 101 rep(i,0,2*N) _A[i]=_A[i]*_B[i]; 102 FFT(_A,A,2*N,-1); 103 rep(i,0,2*N) ans[i+1]=floor(A[i].r+esp); 104 for(int i=2*N;i>=0;i--) if(ans[i]>9){ 105 ans[i-1]+=ans[i]/10;ans[i]%=10; 106 } 107 int cnt=0;while(!ans[cnt]) cnt++; 108 rep(i,cnt,n1+n2-1) printf("%d",ans[i]); 109 //fclose(stdin); 110 //fclose(stdout); 111 return 0; 112 }
1 #include<cstdio> 2 #include<cstring> 3 #include<cmath> 4 #include<ctime> 5 #include<iostream> 6 #include<algorithm> 7 #include<queue> 8 #include<stack> 9 #include<set> 10 #define esp (1e-8) 11 #define maxint (2147483647) 12 #define pi (3.141592653589793) 13 #define l(a) ((a)<<1) 14 #define r(a) ((a)<<1|1) 15 #define b(a) (1<<(a)) 16 #define rep(i,a,b) for(int i=a;i<=(b);i++) 17 #define clr(a) memset(a,0,sizeof(a)) 18 typedef long long ll; 19 using namespace std; 20 int readint(){//只读入正整数 21 int t=0;char c=getchar(); 22 while(!isdigit(c)) c=getchar(); 23 while(isdigit(c)){ 24 t=(t<<3)+(t<<1)+c-'0'; 25 c=getchar(); 26 } 27 return t; 28 } 29 const int maxn=3000009;//注意maxn的取值 30 int n,m; 31 struct complex{ 32 double r,i; 33 inline complex(double R,double I){ 34 r=R;i=I; 35 } 36 inline complex(){} 37 inline void out(){ 38 if(abs(i)>esp) printf("%.2lf %.2lfi\n",r,i); 39 else printf("%.2lf\n",r+esp); 40 } 41 inline complex operator +(const complex A)const{ 42 return complex(r+A.r,i+A.i); 43 } 44 inline complex operator -(const complex A)const{ 45 return complex(r-A.r,i-A.i); 46 } 47 inline complex operator *(const complex A)const{ 48 return complex(r*A.r-i*A.i,i*A.r+r*A.i); 49 } 50 }; 51 int rev(int t,int len){ 52 int ret=0; 53 for(int i=0;b(i)<len;i++){//求逆序->利用性质 54 ret<<=1;if(t&b(i)) ret|=1; 55 } 56 return ret; 57 } 58 complex x1[maxn],x2[maxn],A[maxn],B[maxn]; 59 void DFT(complex Ans[],complex a[],int t,int opt){ 60 rep(i,0,t-1) Ans[rev(i,t)]=a[i]; 61 for(int s=1;b(s)<=t;s++){//当前迭代计算的子列长度为b(s) 62 int m=b(s); 63 complex w1=complex(cos(2*pi/m),sin(2*pi/m*opt)); 64 for(int k=0;k<t;k+=m){//与递归一样的计算过程 65 complex w=complex(1,0); 66 rep(j,0,(m>>1)-1){ 67 complex t1=Ans[k+j],t2=w*Ans[k+j+(m>>1)]; 68 Ans[k+j]=t1+t2;Ans[k+j+(m>>1)]=t1-t2; 69 w=w*w1; 70 } 71 } 72 } 73 if(opt==-1) rep(i,0,t-1) Ans[i].r/=t,Ans[i].i/=t; 74 } 75 int main(){ 76 //freopen("#intput.txt","r",stdin); 77 //freopen("#output.txt","w",stdout); 78 n=readint();m=readint();n++;m++; 79 rep(i,0,n-1) x1[i].r=readint(); 80 rep(i,0,m-1) x2[i].r=readint(); 81 int t=0,N=max(n,m); 82 while(b(t)<N) t++;N=b(t+1); 83 DFT(A,x1,N,1);DFT(B,x2,N,1); 84 rep(i,0,N-1) A[i]=A[i]*B[i]; 85 DFT(B,A,N,-1); 86 rep(i,0,m+n-2){ 87 printf("%.0lf",B[i].r+esp); 88 putchar(i==m+n-2?'\n':' '); 89 } 90 //fclose(stdin); 91 //fclose(stdout); 92 return 0; 93 }
uoj与洛谷都有题
注意几个点：
1、求rev
2、w从（1，0）开始
3、将n化为2^k
cnblogs.com/rvalue/p/7351400.html
