FFT（快速傅里叶变换）

前置知识

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{2\pi}{n}),sin(\frac{2\pi}{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_{2\lfloor\frac{n}{2}\rfloor}x^{2\lfloor\frac{n}{2}\rfloor})+(a_1x+a_3x^3...a_{2\lfloor\frac{n}{2}\rfloor+1}x^{2\lfloor\frac{n}{2}\rfloor+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式只有系数不同，所以可以在k\in[0,\frac{n}{2}-1]时顺便把k\in[\frac{n}{2},n-1]计算出来，就可以实现以下时间\)

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

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

细节

为使序列长度为\(2^m\),把不足位补0
把下标奇偶分类后的位置其实就是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;

}