FFT（快速傅里叶变换）

今天是2018/4/8，DCDCBigBig的第三十四篇博文

啃完算导上的FFT我深深的感觉到自己的数学真是太烂了（因为太菜哇的一声哭了出来）
先膜一发算导QAQ
FFT在信号处理，声音合成等方面有巨大的应用空间，而在OI中主要是用于优化某些式子，用的最广泛的就是快速求两个多项式的积，把时间复杂度从O(n^2)降到O(nlogn)
这里稍微讲一下我的理解（由于本人太菜所以全部没有证明，想要具体证明请自行膜算导）
给出一个次数界为n的多项式A(x)和次数界为m的多项式B(x)，要求它们的积，即多项式C(x)
那么我们很快可以知道C的次数界是n+m（不知道为什么的同学可以按Alt+F4百度一下）
但是直接暴力乘起来的话时间复杂度就会是O(n^2)的，在项数大于10000的时候就GG了，所以要考虑优化，即FFT
FFT的具体做法分三步：
1.对A、B分别进行离散傅里叶变换（即DFT），通俗来讲就是将他们从系数表示转化为用点值表示。以次数界为n的多项式A(x)为例，可以通过n个点(xi,A(xi))(i=0,1,2,…,n-1)来表示，且用此种方法表示的多项式是唯一的
2.通过变换后的A和B求出用点值表示的C，这个可以在线性时间做到，因为对于x的每一个取值都满足C(x)=A(x)*B(x)，所以只需要将离散后的A与B，每一次项分别直接相乘即可。但是注意C的次数界是(n+m)，所以这里实际上需要取A和B的(n+m)个点，可以采取的一个方法是在A和B的高次项补0直到次数为(n+m-1)。（当然计算时间复杂度的时候这些就直接写为n）
3.通过拉格朗日插值法（拉格朗是谁）进行IDFT（逆离散傅里叶变换），通过点值还原回多项式C
但是这样子看复杂度还是O(n^2)的，因为在第一步需要取n个值，对于每个取值需要计算n次（n项），乘起来就是O(n^2)
怎么办呢？此时FFT最精妙的地方就出现了：“通过精妙的取值，使复杂度降为O(nlogn)”
所谓“精妙的取值”，就是应用了单位复数根w
这里的“n次单位复数根”指的是一类满足w^n=1的复数w，且n次单位复数根必然有且仅有n个。因为他们有一些非常妙的性质，从而可以快速计算（请自行百度“消去引理”、“折半引理”和“求和引理”）
这样子取n个点，就可以通过每次折半而变为(logn)个子问题，·对于每一个子问题进行O(n)次计算，就可以达成O(nlogn)的时间复杂度。稍微具体一点就是利用单位复数根在复平面上的平均分布，把每次的子问题分每一项的次数奇偶折半，然后传递到下层子问题（并没有变得具体）
emmm……你以为这就是FFT的全部？并不是，我们只完成了第一步，第二步O(n)不计，那第三步呢？
这时我们需要推一些神秘公式（请自行百度“范德蒙特矩阵”），然后就发现其实IDFT和DFT的公式非常像，只需要把结果和条件颠倒一下，再除以n即可
所以就大功告成了！！！
对了，还有一个非常重要的事情，由于FFT的优化性质，传统意义上的FFT只能解决次数界为2的整数次幂的问题（考虑优化的折半），而对于非2的整数次幂，“策略已存在，但超出了本书范围”%%%算导QAQ
那么就类似于第二步时说的，把A和B高位补0，补到2的整数次幂即可（看起来很暴力，实际上并没有损失多少时间复杂度）
然后在具体实现上，由于C++数组传递的尿性以及递归的巨大常数，我们一般使用非递归方法实现FFT，这时就要引入一个新的东东：rev，翻译一下就是把一个数转化为二进制后每一位取反得到的数值，这个东西对应着FFT递归到最底层时这个数的位置（证明略）
（其实真正的大佬看到这里会发现我什么都没有讲，只是口胡然后浪费了观看者几分钟的时间）
终于讲完辣！！！下面放代码

#include<iostream>
#include<cstring>
#include<cstdio>
#include<cmath>
#define pw(n) (1<<n)
using namespace std;
const double pi=acos(-1);
struct complex{
    double a,b;
    complex(double _a=0,double _b=0){
        a=_a;
        b=_b;
    }
    friend complex operator +(complex x,complex y){return complex(x.a+y.a,x.b+y.b);}
    friend complex operator -(complex x,complex y){return complex(x.a-y.a,x.b-y.b);}
    friend complex operator *(complex x,complex y){return complex(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);}
    friend complex operator *(complex x,double y){return complex(x.a*y,x.b*y);}
    friend complex operator /(complex x,double y){return complex(x.a/y,x.b/y);}
}a[100001],b[100001];
int n,m,bit,bitnum=0,rev[pw(20)];
void getrev(int l){//Reverse
    for(int i=0;i<pw(l);i++){
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(l-1));
    }
}
void FFT(complex *s,int op){
    for(int i=0;i<bit;i++)if(i<rev[i])swap(s[i],s[rev[i]]);
    for(int i=1;i<bit;i<<=1){
        complex w(cos(pi/i),op*sin(pi/i));
        for(int p=i<<1,j=0;j<bit;j+=p){//Butterfly
            complex wk(1,0);
            for(int k=j;k<i+j;k++,wk=wk*w){
                complex x=s[k],y=wk*s[k+i];
                s[k]=x+y;
                s[k+i]=x-y;
            }
        }
    }
    if(op==-1){
        for(int i=0;i<=bit;i++){
            s[i]=s[i]/(double)bit;
        }
    }
}
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++)scanf("%lf",&a[i].a);
    for(int i=0;i<=m;i++)scanf("%lf",&b[i].a);
    m+=n;
    for(bit=1;bit<=m;bit<<=1)bitnum++;
    getrev(bitnum);
    FFT(a,1);
    FFT(b,1);
    for(int i=0;i<=bit;i++)a[i]=a[i]*b[i];
    FFT(a,-1);
    for(int i=m;i>=0;i--)printf("%d ",(int)(a[i].a+0.5));
    return 0;
}

由于现在联赛和省赛中喜(sang)闻(xin)乐(bing)见(kuang)的多了很多数学题，模拟赛也充斥着看不懂的数学名词，本蒟蒻只好被迫来学习数学……而FFT和NTT（下次再讲，这玩意不用大量浮点数计算，所以精度损失小）又是最常用的之一，所以说背好FFT的板子还是很重要的！！！
ps：在这里开一个学莫比乌斯反演、高斯消元、矩阵求逆和欧拉筛的大坑，看看要多久才能达成……