快速傅里叶变换FFT

FFT是基于分治思想对DFT和IDFT的优化。用于多项式的快速乘法

DFT: (mathbb{C} ^{N} ightarrow mathbb{C} ^{N}, left( x_{1},ldots ,x_{N} ight) mapsto left( y_{1},ldots ,y_{N} ight))
IDFT是DFT的逆映射
记(omega _{N}=e^{dfrac {2pi i}{N}})，则

[DFT: y_{k}=sum ^{N-1}_{n=0}omega ^{nk}_{N}x_{n}\ IDFT: x_{k}=sum ^{N-1}_{n=0}dfrac {omega ^{-nk}_{N}}{N}y_{n}]

于是，DFT与IDFT都是N维复空间上的线性变换，其矩阵为Vandermonde矩阵。由于是矩阵乘法，复杂度为O(n^2)。注意到复单位根的性质使得DFT与IDFT高度相似，因此，只需要考虑降低DFT的复杂度即可。

定义：

[E_{k}=sum ^{N/2-1}_{m=0}omega ^{2mk}_{N}x_{2m}=sum ^{N/2-1}_{m=0}omega ^{mk}_{N/2}x_{2m}\ O_{k}=sum ^{N/2-1}_{m=0}omega ^{2mk}_{N}x_{2m+1}=sum ^{N/2-1}_{m=0}omega ^{mk}_{N/2}x_{2m+1}]

则$$y_{k}=E_{k}+omega ^{k}{N}O{k}
y_{k+N/2}=E_{k}-omega ^{k}{N}0{k}$$
显然，如此递归可以将复杂度降至O(nlogn)。另外，N必须是2的幂次。

应用
由多项式的知识知$$left{ fleft( x ight) |deg f=n-1 ight} cong left{ left( a_{0},ldots ,a_{n-1} ight) |a_{n-1} eq 0 ight}
cong left{ left{ left( x_{i},fleft( x_{i} ight) ight) ight} |x_{i} eq x_{j},forall i,j ight}$$
于是，通过FFT可以快速将多项式函数的“系数表示”转化为“点值表示”。
当计算多项式的乘积时，我们知道用“点值表示”计算是方便的。因此可以先将待乘多项式转化为“点值表示”，再将结果转回“系数表示”。

代码（递归）

void separate (complex<double>* a, int n) {
    complex<double>* b = new complex<double>[n/2];
    for(int i=0; i<n/2; i++) 
        b[i] = a[i*2+1];
    for(int i=0; i<n/2; i++) 
        a[i] = a[i*2];
    for(int i=0; i<n/2; i++) 
        a[i+n/2] = b[i];
    delete[] b;
}
// N must be a power-of-2
//DFT: sig=1
//IDFT: sig=-1
// N input samples in X[] are FFT'd and results left in X[]
void fft2 (complex<double>* X, int N,int sig) {
    if(N < 2) {
        // bottom of recursion.
    } else {
        separate(X,N);
        fft2(X,     N/2,sig);
        fft2(X+N/2, N/2,sig);
        // combine results of two half recursions
        for(int k=0; k<N/2; k++) {
            complex<double> e = X[k    ];   // even
            complex<double> o = X[k+N/2];   // odd
                         // w is the "twiddle-factor"
            complex<double> w = exp( complex<double>(0,sig*2.*M_PI*k/N) );
            X[k    ] = e + w * o;
            X[k+N/2] = e - w * o;
        }
    }
}

模拟递归。
参考网上的博客，这种模拟首先需要做一种置换：对指标二进制表示下对称的元素做对换。（至于如何实现这一置换，可以反向做二进制加法。例如1100的后一位是0010，正如0011的后一位是0100一样。过程见代码。）接着，如同递归函数从栈顶到栈底的过程一样，做模拟的递归过程。

void rearrange(complex<double> x[],int n)
{
    for(int i=1,j=n/2;i<n;++i)
    {
        if(i<j)swap(x[i],x[j]);
        int tmp=n/2;
        while(tmp&&j>=tmp){j-=tmp;tmp/=2;}
        if(j<tmp)j+=tmp;
    }
}
void fft(complex<double> x[],int n,int sig)
{
    rearrange(x,n);
    for(int i=2;i<=n;i*=2)
    {
        for(int j=0;j<n;j+=i)
        {
            for(int k=j;k<j+i/2;++k)
            {
                complex<double> e=x[k],o=x[k+i/2],w=exp(complex<double>(0,sig*2.*pi*(k-j)/i));
                x[k]=e+w*o;
                x[k+i/2]=e-w*o;
            }
        }
    }
    if(sig==-1)
    {
        for(int i=0;i<n;++i)x[i]/=n;
    }
}