luogu3803 多项式乘法 (FFT)

Posted 2020-12-25 ressed

FFT讲解传送门

简单记一下做法：

　　1.算法流程：两式的系数表达转化为点值表达（O(nlogn)）->利用点值表达使两式相乘（O(n)）->将结果的点值表达转化回系数表达（O(nlogn)）

　　2.做法：

　　　　$$目标：把一个n项多项式F(x)=sum_{i=0}^{n-1}a_ix^i转化为{(w^k_n,y_k)}的点值表达，其中w^k_n为n次单位根的k次方$$

　　　　不妨设n为2的幂次，如果不是，则可以补上系数为0的高次项

　　　　$$将a_ix^i按照幂次奇偶性分组，得到F(x)=(a_0x^0+a_2x^2+...+a_{n-2}x^{n-2})+(a_1x^1+a_3x^3..a_{n-1}x^{n-1})$$

　　　　$$设F_0(x)=a_0x^0+a_2x^1+a_4x^2+...+a_{n-2}x^{frac{n}{2}-1} , F_1(x)=a_1x^0+a_3x^1+a_5x^2+...+a_{n-1}x^{frac{n}{2}-1}$$

　　　　$$则F(x)=F_0(x^2)+xF_1(x^2)$$

　　　　$$带入n次单位根，F(w^k_n)=F_0(w^{2k}_n)+w^k_nF_1(w^{2k}_n)$$

　　　　$$w^{xk}_{xn}=w^k_n  ==> F(w^k_n)=F_0(w^k_{frac{n}{2}})+w^k_nF_1(w^k_{frac{n}{2}}) $$

　　　　$$w^{k}_{n}=-w^{k+frac{n}{2}}_n ==> F(w^{k+frac{n}{2}}_n)=F_0(w^k_{frac{n}{2}})-w^k_nF_1(w^k_{frac{n}{2}}) $$

剩下明天写（2018-07-25）

 

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #include<cmath>
 5 #define LL long long int
 6 using namespace std;
 7 const int maxn=2097152;
 8 const double Pi=acos(-1);
 9 
10 LL rd(){
11     LL x=0;int neg=1;char c=getchar();
12     while(c<‘0‘||c>‘9‘){if(c==‘-‘) neg=-1;c=getchar();}
13     while(c>=‘0‘&&c<=‘9‘) x=x*10+c-‘0‘,c=getchar();
14     return x*neg;
15 }
16 
17 struct Complex{
18     double x,y;
19     Complex(double x0=0,double y0=0){x=x0,y=y0;}
20 }A[maxn],B[maxn];
21 Complex operator * (Complex a,Complex b){return Complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
22 Complex operator + (Complex a,Complex b){return Complex(a.x+b.x,a.y+b.y);}
23 Complex operator - (Complex a,Complex b){return Complex(a.x-b.x,a.y-b.y);}
24 int rev[maxn],N,M,len=1;
25 
26 void fft(Complex *a,int opt){
27     for(int i=0;i<len;i++){
28         if(i<rev[i]) swap(a[i],a[rev[i]]);
29     }
30     for(int i=1;i<len;i<<=1){
31         Complex wn=Complex(cos(Pi/i),opt*sin(Pi/i));
32         int step=i<<1;
33         for(int j=0;j<len;j+=step){
34             Complex w=Complex(1,0);
35             for(int k=0;k<i;k++,w=w*wn){
36                 Complex x=a[j+k];
37                 Complex y=w*a[j+k+i];
38                 a[j+k]=x+y;
39                 a[j+k+i]=x-y;
40             }
41         }
42     }
43 }
44 
45 int main(){
46     int i,j,k;
47     N=rd()+1,M=rd()+1;
48     for(i=0;i<N;i++) A[i].x=rd();
49     for(i=0;i<M;i++) B[i].x=rd();
50     k=0;while(len<N+M-1) len<<=1,k++;
51     for(i=0;i<len;i++) rev[i]=(rev[i>>1]>>1)|((i&1)<<(k-1));
52     fft(A,1);fft(B,1);
53     for(i=0;i<len;i++) A[i]=A[i]*B[i]; 
54     fft(A,-1);
55     for(i=0;i<N+M-1;i++){
56         printf("%d ",(int)(A[i].x/len+0.5));
57     }
58     return 0;
59 }