回来第一次复习FFT

今天就写一下FFT的小结

在此安利我所复习的博客,stO 自为风月马前卒的FFT详解

前置知识

1.复数相关

复数的表示法

$a+bi$

复数加减乘

$$(a+bi)+(c+di)=(a+c)+(b+d)i$$

$$(a+bi)-(c+di)=(a-c)+(b-d)i$$

$$(a+bi)*(c+di)=(ac-bd)+(ad+bc)i$$

复数可在直角坐标系中表示

2.点值表示法

如一个多项式

他的系数表示法为

$$A_n=a_0+a_1x+a_2x^2+...+a_nx^n$$

假定这么一个多项式

$$A_n=2+5x+3x^2$$

将若干个 $x$ 代入多项式中,得到一系列点 $(x_i,y_i)$

比如这个多项式,我们将 $x=0,x=1,x=2$ 代入

得到 $(0,2),(1,10),(2,24)$ 这一系列点

我们也可以用这一系列点来表示这个多项式

为什么要用到这个呢

在正常的多项式乘法中

若是两个多项式相乘

需要将两个多项式中每个项都乘一遍,复杂度 $o(n^2)$

若是采用点值表示法,只需要将每两个对应的点相乘,复杂度 $o(n)$

3.1.单位根

以下内容默认n为2的整数次幂

在复数平面上,单位圆就是以原点为圆心,半径为1的一个圆

将这个圆n等分,按照分割线作n个向量

设幅度为正最小的那个向量为 $\omega_n$

我们称之为n次单位根

可以知道其他n-1个向量分别为 $\omega_n^2,\omega_n^3,...,\omega_n^n$

其中 $\omega_n^n=\omega_n^0=1$

这些单位根的值可以通过欧拉公式计算

$$\omega_n^k=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}$$

在代数中,若 $z^k=1$ ,我们称 $z$ 为 $k$ 次单位根

3.2.单位根的性质

$$\omega_n^k=\cos \frac{2k\pi}{n}+i\sin \frac{2k\pi}{n}$$

$$\omega_{2n}^{2k}=\omega_n^k$$

$$\omega_n^{k+\frac{n}{2}}=-\omega_n^k$$

$$\omega_n^n=\omega_n^0=1$$

$k\not= 0$ 时

$$\sum_{i=1}^n\omega_n^{ki}=0$$

$k=0$ 时

$$\sum_{i=1}^n\omega_n^{ki}=n$$

快速傅里叶变换(FFT)(Fast Fourier Transfromation)

我们设一个多项式

$$A(x)=a_0+a_1x+a_2x^2+...+a_{n-1}^{n-1}$$

将其按照奇偶拆开

$$=(a_0+a_2x^2+....+a_{n-2}x^{n-2})+(a_1x+a_3x^3+...+a_{n-1}x^{n-1})$$

奇数项提出一个x

$$=(a_0+a_2x^2+....+a_{n-2}x^{n-2})$$

$$+x(a_1+a_3x^2+...+a_{n-1}x^{n-2})$$

$$A_1(x)=a_0+a_2x+....+a_{n-2}x^{\frac{n}{2}-1}$$

$$A_2(x)=a_1+a_3x+...+a_{n-1}x^{\frac{n}{2}-1}$$

可知

$$A(x)=A_1(x^2)+xA_2(x^2)$$

将 $\omega_n^k$ 代入

$$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+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}A_2(\omega_n^{2k+n})$$

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

我们将两个式子比较

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

$$A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})$$

发现枚举其中一个时,可以在 $o(1)$ 的时间将另一个枚举出来

每次我们可以将问题缩小一半

这样就类似线段树的时间复杂度 $o(nlogn)$

本身需要 $o(n^2)$ 的时间,我们拿 $o(n+2nlogn)$ 就完成了,得到了优化

快速傅里叶逆变换(IFFT)(Fast Fourier Inverse Transformation)

但是,我们常用的表示法并不是点值表示法

我们还应该将点值表示法转化为系数表示法

这个过程叫做傅里叶逆变换

首先

我们设 $(y_1,y_2,...,y_n)$ 为 $(a_1,a_2,...,a_n)$ 的点值表示

设另一个向量 $(c_1,c_2,...,c_n)$ 满足

$$c_k=\sum_{i=0}^{n-1}y_i(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^{i})^j)(\omega_n^{-k})^i$$

$$=\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^{i})^j)(\omega_n^{-k})^i$$

$$=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}$$

于是我们知道

$j=k$ 时

$$\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}=n$$

而 $j\not=k$ 时

$$\sum_{i=0}^{n-1}(\omega_n^i)^{(j-k)}=0$$

于是

$$c_k=na_k$$

$$a_k=\frac{c_k}{n}$$

所以我们对乘起来的点值表示再做一次FFT

然后除去n就是所得到的多项式了

所用时间复杂度为 $o(2nlogn+n+nlogn)$ ,即 $o(nlogn)$

迭代实现

我们也可以尝试用递归的方法来完成FFT,但这样并不是最优的办法

如果我们事先就将需要转换的位置弄明白,这样可以直接转换

可以发现

反转后的序列的二进制数也反转过来

于是我们可以预处理我们需要得到的数列

就可以在 $o(n)$ 的时间复杂度内完成了

板子题

【模板】多项式乘法

#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<iostream>

using namespace std;

#define fr(i,a,b) for(register long long i=a;i<=b;i++)

const double Pi=acos(-1.0);

struct complex{
    double x,y;
    complex(double xx=0,double yy=0){x=xx,y=yy;}
}a[5000010],b[5000010];
complex operator + (complex a,complex b){return complex(a.x+b.x,a.y+b.y);}
complex operator - (complex a,complex b){return complex(a.x-b.x,a.y-b.y);}
complex operator * (complex a,complex b){return complex(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}

long long n,m,limit,l;
long long r[5000010];

inline long long v_in(){
    char ch=getchar();long long sum=0,f=1;
    while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
    while(isdigit(ch)) sum=(sum<<3)+(sum<<1)+(ch^48),ch=getchar();
    return sum*f;
}

void FFT(complex *A,long long type){
    fr(i,0,limit-1) if(i<r[i]) swap(A[i],A[r[i]]);
    for(register long long mid=1;mid<limit;mid<<=1){
        complex Wn(cos(Pi/mid),type*sin(Pi/mid));
        for(register long long j=0;j<limit;j+=(mid<<1)){
            complex W(1,0);
            for(register long long k=0;k<mid;k++,W=W*Wn){
                complex x=A[j+k],y=W*A[j+mid+k];
                A[j+k]=x+y;
                A[j+mid+k]=x-y;
            }
        }
    }
}

void inpt(){
    n=v_in(),m=v_in();
    fr(i,0,n) a[i].x=v_in();
    fr(i,0,m) b[i].x=v_in();
    limit=1,l=0;
    while(limit<n+m+1) limit<<=1,l++;
    fr(i,0,limit-1) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
}

void work(){
    FFT(a,1);
    FFT(b,1);
    fr(i,0,limit) a[i]=a[i]*b[i];
    FFT(a,-1);
}

int main(){
    inpt();
    work();
    fr(i,0,n+m) printf("%lld ",(int)(a[i].x/limit+0.5));
    printf("\n");
    return 0;
}