关于快速傅里叶变换的一些总结

以前打FFT只想到背板就行了，现在看来还是有必要了解一下其原理，故写出此篇总结。

概述

FFT是一类用来计算多项式乘积的算法，如果用普通的方法进行，复杂度将会是$O(n^2)$,而FFT可以做到$O(n\log n)$时间复杂度内解决。

多项式及其乘法

多项式的点值表示
通过拉格朗日插值的唯一性可知，一个$n$次多项式由$n+1$个函数上的点确定。

已知多项式$A(x)=\sum_{i=0}^na_ix^i$,$B(x)=\sum_{i=0}^mb_ix^i$,要求得两多项式的乘积$C(x)=A(x)B(x)$,首先知道$C(x)$是一个$n+m$次多项式，那么选取$n+m+1$个点值即可确定，而FFT就是加速求取点值的一种方法。

Cooley-Tukey算法

首先将这个多项式的次数补充至$2^n-1$次（高位补充0），现在来介绍这一算法的过程。
现在设n为这一多项式的次数。
先将函数域扩展到虚数域上，问题转化为求取单位根$\omega_n^0,\omega_n^1...\omega_n^{n-1}$（后面会发现，选取这些数的原因是因为他们彼此重复计算了很多部分，而且可以快速整合）。

有：

$$F(\omega_n^k)=\sum_{i=0}^n{a_i\omega_n^{ik}}$$

考虑将系数奇偶分类，有

$$F(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ik}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ik}$$

同时有

$$F(\omega_n^{k+\frac{n}{2}})=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ik}+\omega_n^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ik}$$
容易发现如果一开始系数按照二进制分组的话，两边都是一个求解的子问题，而且最多往下 $\log$层，所以可以在 $n\log n$时间内解决求取点值问题。

接下来考虑求出原函数(以下证明来自:传送门)：

应用

1.多项式乘法(uoj34)
这个是一个裸题了，先放个板。

#include<bits/stdc++.h>
typedef double db;
using namespace std;

inline int read(){
    char ch=getchar();int i=0,f=1;
    while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
    while(isdigit(ch)){i=(i<<1)+(i<<3)+ch-'0';ch=getchar();}
    return i*f;
} 
inline void W(int x){
    static int buf[50];
    if(!x){putchar('0');return;}
    if(x<0){putchar('-');x=-x;}
    while(x){buf[++buf[0]]=x%10;x/=10;}
    while(buf[0])putchar(buf[buf[0]--]+'0');
}
const int Maxn=4e5+50;
const double PI=acos(-1.0);
int n,m;
struct cp{
    db r,i;
    cp(){}
    cp(db r,db i):r(r),i(i){}
    friend inline cp operator +(const cp &a,const cp &b){return cp(a.r+b.r,a.i+b.i);}
    friend inline cp operator -(const cp &a,const cp &b){return cp(a.r-b.r,a.i-b.i);}
    friend inline cp operator *(const cp &a,const cp &b){return cp(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}
    friend inline cp operator /(const cp &a,db b){return cp(a.r/b,a.i/b);}
}A[Maxn],B[Maxn];
struct FFT{
    int pos[Maxn],k;cp w[Maxn];
    inline void init(int len){
        for(k=1;k<len;k<<=1);
        cp wn(cos(PI*2/k),sin(PI*2/k));
        for(int i=1;i<k;i++)pos[i]=((i&1)?((pos[i>>1]>>1)^(k>>1)):(pos[i>>1]>>1)); 
        w[0]=cp(1,0);
        for(int i=1;i<=k;i++)w[i]=w[i-1]*wn; 
    }
    inline void dft(cp *a){
        for(int i=1;i<k;i++)if(pos[i]>i)swap(a[pos[i]],a[i]);
        for(int i=1;i<k;i<<=1){
            int bl=(i<<1),ct=k/bl;
            for(int j=0;j<k;j+=bl){
                for(int z=0;z<i;++z){
                    cp &t1=a[j+z],&t2=a[j+z+i],t=t2;
                    t2=t1+w[(z+i)*ct]*t,t1=t1+w[z*ct]*t;
                }
            } 
        }
    }
    inline void mul(cp *a,cp *b,cp *c){
        dft(a),dft(b);
        for(int i=0;i<k;i++)c[i]=a[i]*b[i];
        dft(c);reverse(c+1,c+k);
        for(int i=0;i<k;i++)c[i]=c[i]/(double)k; 
    }
}fft;
int main(){
    n=read(),m=read();
    for(int i=0;i<=n;i++)A[i].r=read();
    for(int j=0;j<=m;j++)B[j].r=read();
    fft.init(n+m+1);fft.mul(A,B,A);
    for(int i=0;i<=n+m;i++)W((int)(A[i].r+0.5)),putchar(' ');
}

2.循环卷积。

详见:http://blog.csdn.net/qq_35649707/article/details/78660144