FFT（快速傅里叶变换）

Posted 2021-03-22 sydevil

FFT（快速傅里叶变换）

前置知识

(1.复数)

(2.单位根)

(3.循环结构)

(4.C++)

1.复数

• (定义:形如a+bi的数，其中i^2=-1)
• (计算:1.(a+bi)+(c+di)=(a+c)+(b+d)i)
  
  ( 2.(a+bi)-(c+di)=(a-c)+(b-d)i)
  
  ( 3.(a+bi)(c+di)=(ac-bd)+(ad+bc)i)

• (向量表示法)
  
  (用x表示实部(没有i的部分)，用y表示虚部(i前的系数))

  (性质1.复数乘法满足向量模角相加)

  (证明:arctan(frac{a}{b})+arctan(frac{c}{d})=arctan(frac{frac{a}{b}+frac{c}{d}}{1-frac{ac}{bd}})=arctan(frac{ad+bc}{bd-ac}))

  (是不是巨简单！)

2.单位根

• (定义:n次单位根是n次幂为1的复数。它们位于复平面的单位圆上，)
  (构成正n边形的顶点，其中一个顶点是1。)

• (显然,omega_n^0=omega_n^n=1,omega_n=(cos(frac{2pi}{n}),sin(frac{2pi}{n})),n个复数为omega_n^1,omega_n^2...omega_n^n)

• (性质：omega_{n}^k=omega_{2n}^{2k})

正题

(令A(x)=a_0+a_1x^1...a_nx^n)

(则A(x)=(a_0+a_2x^2...a_{2lfloorfrac{n}{2} floor}x^{2lfloorfrac{n}{2} floor})+(a_1x+a_3x^3...a_{2lfloorfrac{n}{2} floor+1}x^{2lfloorfrac{n}{2} floor+1}),令左边的为A_1(x^2),右边的为xA_2(x^2),A(x)=A_1(x^2)+xA_2(x^2))

(将omega_n^k(k<frac{n}{2})带入得)

(A(omega_n^k)=A_1(omega_n^{2k})+omega_n^kA_2(omega_n^{2k}))

(将omega_n^{k+frac{n}{2}}带入得)

(A(omega_n^k)=A_1(omega_n^{2k})-omega_n^kA_2(omega_n^{2k}))

(我们发现2式只有系数不同，所以可以在kin[0,frac{n}{2}-1]时顺便把kin[frac{n}{2},n-1]计算出来，就可以实现以下时间)

[ T(n)=left{ egin{array}{rcl} 2T(frac{n}{2})+n&& {1<n}1&& {n=1} end{array} ight. ]

(时间复杂度为O(nlog_2n))

---

细节

1. 为使序列长度为(2^m),把不足位补0

2. 把下标奇偶分类后的位置其实就是2进制下的反转，所以可以通过迭代（非递归）实现，实际效率也是递归版的(1.5)~(2)倍

(mathfrak{Talk is cheap,show you the code.})

#include<cstdio>
#include<cmath>
#include<vector>
#include<algorithm>
using namespace std;
# define read read1<int>()
# define Type template<typename T>
Type inline T read1(){
    T n=0;
    char k;
    bool fl=0;
    do (k=getchar())=='-'&&(fl=1);while('9'<k||k<'0');
    while(47<k&&k<58)n=(n<<3)+(n<<1)+(k^48),k=getchar();
    return fl?-n:n;
}
# define f(i,l,r) for(int i=(l);i<=(r);++i)
# define fre(k) freopen(k".in","r",stdin);freopen(k".ans","w",stdout)
const double PI=acos(-1);
class complex{
    public:
        double x,y;
        complex(){x=y=0;}
        complex(double _x,double _y):x(_x),y(_y){}
        complex operator * (complex b){return complex(x*b.x-y*b.y,x*b.y+y*b.x);} 
        complex operator + (complex b){return complex(x+b.x,y+b.y);}
        complex operator - (complex b){return complex(x-b.x,y-b.y);}
        complex operator * (double u){return complex(x*u,y*u);} 
        complex& operator *= (complex b){return *this=*this*b;} 
        complex& operator += (complex b){return *this=*this+b;} 
        complex& operator -= (complex b){return *this=*this-b;} 
};
class Array{
    private:
        vector<int>a;
    public:
        Array(){}
        void push(int n){a.push_back(n);}
        Array(int* l,int* r){while(l!=r)push(*l),++l;}
        int size(){return a.size();}
        int& operator [] (const int x){return a[x];}
};
void FFT(const int len,vector<complex>&a,const int Ty,int *r=NULL){
    if(!r){
        r=new int[len];
        r[0]=0;int L=log2(len);
        f(i,0,len-1){
            r[i]=(r[i>>1]>>1)|((i&1)<<L-1);
            if(i<r[i])swap(a[i],a[r[i]]);
        }
    }
    for(int i=1;i<len;i<<=1){
        complex T(cos(PI/i),Ty*sin(PI/i));
        for(int W=i<<1,j=0;j<len;j+=W){
            complex n(1,0);
            for(int k=0;k<i;++k,n*=T){
                complex x(a[j+k]),y(n*a[i+j+k]);
                a[j+k]=x+y;
                a[i+j+k]=x-y;
            }
        }
    }
}
Array operator * (Array x,Array y){
    int n=x.size()-1,m=y.size()-1;
    int limit=1;
    while(limit<=n+m)limit<<=1;
    vector<complex>_x(limit+1),_y(limit+1);
    Array ans;
    f(i,0,n)_x[i]=complex(x[i],0);
    f(i,0,m)_y[i]=complex(y[i],0);
    FFT(limit,_x,1);
    FFT(limit,_y,1);
    f(i,0,limit)_x[i]*=_y[i];
    FFT(limit,_x,-1);
    f(i,0,n+m)ans.push((int)(_x[i].x/limit+0.5));
    return ans;
}
void into(int n,Array &x){
    f(i,0,n)x.push(read);
}
Array x,y,ans;
int main(){
    int n=read,m=read;
    into(n,x),into(m,y);
    ans=x*y;
    f(i,0,n+m)printf("%d ",ans[i]);
    return 0;
}