Skip to content

快速傅里叶变换 FFT

~~众所周知，FFT全称是Fast Fast TLE~~

这是一个十分巧妙且优美的算法，包括了分治迭代等重要思想。

前置知识

具体参考这篇文章,里面讲的很详细，包括多项式乘法(卷积)，系数表示法与点值表示法，朴素算法(离散傅里叶变换DFT) ，复数加法乘法，单位根及其性质等诸多知识点。

FFT

这里只讲一些具体的，原文没讲清楚的细节。

  1. 首先，假设原来的系数是\(a_0,a_1,...,a_n\),求出来的点值表示为\((x_0,y_0),(x_1,y_1),...,(x_n,y_n)\),再将\(y_0,y_1,...,y_n\)作为系数逆变换出点值表示\((x_0',c_0),(x_1',c_1),...,(x_n',c_n)\),则此时满足\(a_n=\frac{c_n}{k}\)

~~作为蒟蒻看原文竟然看了三天才看懂，所以解释一下~~

  1. 假设多项式\(a,b\)为\(n,m\)次多项式，则卷积卷出来的多项式为\(n+m\)次，因此在算法中，需要将\(a,b\)两个多项式的系数补到\(n+m\)次。就算系数是\(0\)也没事。

  2. 原文中"问题缩小一半，可以分治"体现为：假设算法进行中已经算出\(a_1,a_2\)两个子多项式的点值表示，且他们的点值分别对应多项式\(a_1,a_2\)取\(W_n^0 ,W_n^1,...,W_n^k\)时的值，那么可以用这些值算出多项式\(a\)取\(W_n^0 ,W_n^1,...,W_n^{2k}\)的值，这样从\(k=1\)一直到\(k=n\)，我们用\(O(n\log n)\)的时间复杂度算出了整个多项式的点值表示。

另外，因为递归的空间开销过高，所以要采用迭代(递推)的方式实现。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
using namespace std;
typedef double db;
const int N=4e6+10;
const db Pi=acos(-1.0);
struct c{
    db x,y;
    c(db xx=0,db yy=0){x=xx,y=yy;}
}a[N],b[N];
c operator +(c a,c b){return c(a.x+b.x,a.y+b.y);}
c operator -(c a,c b){return c(a.x-b.x,a.y-b.y);}
c operator *(c a,c b){return c(a.x*b.x-a.y*b.y,a.y*b.x+b.y*a.x);}
int n,m,k,l;
int r[N];
void init(){
    l=0,k=1;while(k<=n+m)k<<=1,++l;
    for(int i=0;i<k;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));//找规律的东西就背一背好了
}
void FFT(c* a,int t){
    for(int i=0;i<k;++i)if(i<r[i])swap(a[i],a[r[i]]);//得到最底层的多项式系数
    for(int mid=1;mid<k;mid<<=1){//枚举中点，即区间长的一半
        c Wn(cos(Pi/mid),t*sin(Pi/mid));//单位根
        for(int L=mid<<1,j=0;j<k;j+=L){//L为区间长，j为区间起点
            c w(1,0);//p=0的初始状态
            for(int p=0;p<mid;++p,w=w*Wn){//这里的p就是单位根的幂次，为0~mid
                c x=a[j+p],y=a[j+mid+p]*w;//蝴蝶操作
                a[j+p]=x+y;//前半
                a[j+mid+p]=x-y;//后半
            }
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;++i)scanf("%lf",&a[i].x);
    for(int i=0;i<=m;++i)scanf("%lf",&b[i].x);
    init();
    FFT(a,1);//正变换
    FFT(b,1);
    for(int i=0;i<k;++i)a[i]=a[i]*b[i];//点值直接O(n)相乘
    FFT(a,-1);//逆变换
    for(int i=0;i<=n+m;++i)printf("%d ",(int)(a[i].x/k+0.5));
    return 0;
}