FFT食用指北

概述

快速傅里叶变换(Fast Fourier Transform,简称FFT)在OI中主要作用是求多项式乘法

朴素的多项式乘法即为系数相乘,时间复杂度是的。可以通过FFT降为


预备知识

多项式的表示

对于n次多项式常用的表示方法有两种:
系数表示法:即 即用n+1个系数表示一个多项式
点值表示法:选取互不相同,用二元组 表示一个多项式 其中

不难发现这两种方法都是可以表示一个唯一确定的多项式

复数

为实数, ,形如 的数叫做复数,其中称为虚数单位。复数域是已知的最大的域

复平面

在复平面中,轴代表实数,轴代表虚数。每一个复数对应一个从的向量。

该向量的长度叫做模长。表示从 轴正半轴到该向量的转角的有向(以逆时针为正方向)角叫做幅角。

复数相加遵循平行四边形定则。

复数相乘时,模长相乘,幅角相加。

单位复数根

n次单位复数根是满足的复数次单位复数根恰好有个。对于,这些根是

欧拉公式来解释这个表达式

这个公式是可以基于泰勒展开证明,此处略去。

以复平面的原点为圆心的单位半径的圆周上。值 称为主n次单位根(这里参考的是算导的定义),所有其他的n次单位复数根都是的幂次

主n次单位根在复平面上就是轴正半轴的单位向量逆时针旋转 的向量

单位根的性质

单位根有很多有用的性质,其中 有很多联系

比如

其中性质一可以通过几何意义不难看出

性质二的证明如下

考虑其值

考虑几何意义即为旋转 后的向量 即实轴取反,虚轴取反。

性质三更简单了,直接把幂次展开就好


多项式乘法

如果都是次数为的多项式

则他们的乘积

朴素的乘法显然是的,我们发现如果用点值表示法做乘法显然复杂度是的如果能够找到一种方法快速从点值表示法和系数表示法相转换,我们就可以快速地计算多项式的乘法了,快速傅里叶变换就可以做到这一点。

快速傅里叶变换

现在的思路就是找一些特殊的点,求出点值表示法,再用的复杂度得出乘完之后的点值表示法,之后再转化为系数表示法

考虑计算的点值 其中的幂次

的偶数次项系数组成的多项式, 是奇数。那么有:

另外

所以说,如果我们有了的点值表示,就可以在的时间内求的的点值表示,而又是原问题的一半。因为我们保证了的幂次,所以我们可以采用分治算法。

根据主定理 时间复杂度为

傅里叶逆变换(IDFT)

有了求得点值的变换,还需要通过点值恢复出多项式的变换。

考虑变换的本质,等价于把系数列向量左乘一个矩阵得到点值列向量。并且这个矩阵是个非常特殊的范德蒙德矩阵,第 行的公比是

记上面的系数矩阵为现在考虑下面这个矩阵

设它们相乘的结果为
则有:

时,显然有
时 根据等比数列求和公式,有

容易发现其中表示单位矩阵。那么有

那么求系数的方法就是把变成 做一次类似快速傅里叶变换的过程,再将结果每个数除以 ,即为傅里叶逆变换的结果。

高效FFT实现

FFT的递归实现

只是模拟我们推出的过程即可。
还有讲课的老师说STL的提供了complex的模板
但是不建议用,手写的也不多。
被某个毒瘤卡STL岂不是很GG

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cctype>
#include<cmath>
using namespace std;

inline int read(){
    char c=getchar();int x=0,f=1;
    while(!isdigit(c)){
        if(c=='-')f=-1;
        c=getchar();
    }
    while(isdigit(c)){
        x=(x<<3)+(x<<1)+(c^48);
        c=getchar();
    }
    return x*f;
}

const double pi=acos(-1.0);
const int N=2e6+10;
struct complex{
    double x,y;
    complex(double tx=0,double ty=0){
        x=tx;y=ty;
    }
    complex operator +(const complex &o)const{
        return complex(x+o.x,y+o.y);
    }
    complex operator -(const complex &o)const{
        return complex(x-o.x,y-o.y);
    }
    complex operator *(const complex &o)const{
        return complex(x*o.x-y*o.y,x*o.y+y*o.x);
    }
}a[N],b[N];

inline void FFT(int limit,complex *a,int type){
    if(limit==1)return ;
    complex a0[limit>>1],a1[limit>>1];
    for(int i=0;i<=limit;i+=2){
        a0[i>>1]=a[i];
        a1[i>>1]=a[i+1];
    }
    FFT(limit>>1,a0,type);
    FFT(limit>>1,a1,type);
    complex Wn=complex(cos(2.0*pi/limit),type*sin(2.0*pi/limit)),w=complex(1,0);
    for(int i=0;i<(limit)>>1;i++,w=w*Wn){
        a[i]=a0[i]+w*a1[i];
        a[i+(limit>>1)]=a0[i]-w*a1[i];
    }
}

int main(){
    int n=read(),m=read();
    for(int i=0;i<=n;i++)a[i].x=read();
    for(int i=0;i<=m;i++)b[i].x=read();
    int limit=1;
    while(limit<=n+m)limit<<=1;
    FFT(limit,a,1);
    FFT(limit,b,1);
    for(int i=0;i<=limit;i++)
    a[i]=a[i]*b[i];
    FFT(limit,a,-1);
    for(int i=0;i<=n+m;i++)
    printf("%d ",(int)(a[i].x/limit+0.5));
    return 0;
} 

蝴蝶变换

看样子这个东西的名字很NB
但是十分的简单。撕拷后发现复数的运算是很慢的。而我们将的值记录下来就像这样。

    for(int i=0;i<(limit)>>1;i++,w=w*Wn){
        complex tmp=w*a1[i];
        a[i]=a0[i]+tmp;
        a[i+(limit>>1)]=a0[i]-tmp;
    }

FFT的迭代实现

FFT的递归实现常数巨大,会TLE
考虑一下每次递归后的边界
盗图

一个显然的发现
我们求的序列是原序列的二进制反转。
不用按奇偶性分类只需要找到我们需要合并的区间不断向上合并即可。

具体操作是这样的

求数的二进制反转

  for(int i=0;i<limit;i++)
        r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ;

可以这么理解

在原序列中的的关系是可以看成左移一位得到的结果。那么在反转后的序列就要右移一位。
同时要处理一下奇数。

有了上一步将数转换为其二进制反转的操作,并求出其数组,之后的步骤主要是

将数列变换为递归实现的最后一步的序列

for(int i=0;i<limit;i++) 
        if(i<r[i]) swap(A[i],A[r[i]]);

有了这一步就可以将原数列转换为上图递归实现的最后一步。

迭代操作

讲解在代码中了。
可能有锅,请一定要联系我

void FFT(complex *A,int type){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(A[i],A[r[i]]);
    //将数列变换为递归实现的最后一步的序列
    for(int cur=1;cur<limit;cur<<=1){
    //这里枚举的cur是将要计算的区间的一半 
        complex Wn(cos(pi/cur),type*sin(pi/cur));
        //用欧拉公式求出单位根 
        int len=cur<<1;
        //这里的len表示计算的区间长度 
        for(int j=0;j<limit;j+=len){
        //这里枚举每一个即将要操作的区间,其中的j表示这段区间的左端点 
            complex w(1,0);
            for(int k=0;k<cur;k++,w=w*Wn){
            //具体枚举区间的每一个数 
                complex x=A[j+k],y=w*A[j+cur+k];
                //蝴蝶变换 
                A[j+k]=x+y;
                A[j+cur+k]=x-y;
            }
        }
    }
}

板子

#include<cstdio>
#include<cctype>
#include<cmath>
#include<algorithm>
using namespace std;

inline int read(){
    char c=getchar();int x=0,f=1;
    while(!isdigit(c)){
        if(c=='-')f=-1;
        c=getchar();
    } 
    while(isdigit(c)){
        x=(x<<3)+(x<<1)+(c^48);
        c=getchar();
    }
    return x*f;
}

const int N=1e7+29;
const double pi=acos(-1.0);
struct complex{
    double x,y;
    complex(double tx=0,double ty=0){
        x=tx;y=ty;
    }
    complex operator +(const complex &o)const{
        return complex(x+o.x,y+o.y);
    }
    complex operator -(const complex &o)const{
        return complex(x-o.x,y-o.y);
    }
    complex operator *(const complex &o)const{
        return complex(x*o.x-y*o.y,y*o.x+x*o.y);
    }
}a[N],b[N];
int limit=1,l,r[N];

inline void FFT(complex *A,int type){
    for(int i=0;i<limit;i++)
        if(i<r[i])
            swap(A[i],A[r[i]]);
    for(int cur=1;cur<limit;cur<<=1){
        complex Wn(cos(pi/cur),type*sin(pi/cur));
        int len=cur<<1;
        for(int j=0;j<limit;j+=len){
            complex w(1,0);
            for(int k=0;k<cur;k++,w=w*Wn){
                complex x=A[j+k],y=w*A[j+cur+k];
                A[j+k]=x+y;
                A[j+cur+k]=x-y;
            }
        }
    }
}

int main(){
    int n=read(),m=read();
    for(int i=0;i<=n;i++)
        a[i].x=read();
    for(int i=0;i<=m;i++)
        b[i].x=read();
    while(limit<=n+m){
        limit<<=1;
        l++;
    }
    for(int i=0;i<limit;i++)
        r[i]=((r[i>>1]>>1)|((i&1)<<(l-1)));
    FFT(a,1);
    FFT(b,1);
    for(int i=0;i<=limit;i++)
        a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=0;i<=m+n;i++)
        printf("%d ",(int)(a[i].x/limit+0.5));
    return 0;
} 

参考资料

发表评论

电子邮件地址不会被公开。 必填项已用*标注