快速傅里叶变换（FFT）（上）

FFT这东西觉得很玄学啊，绝大部分人看了也是一知半解，似懂非懂，或者根本不懂.....

我今天想尝试一下解释这FFT的很玄学的理解（因为不一定正确）.....

傅里叶级数:
在讲FFT之前，要说清楚一下傅里叶变换。
傅里叶变换最早要追朔到傅里叶级数,傅里叶级数其实和幂级数其实是同一个玩意。
幂级数是说,我们使用x^0,x^1,x^2,x^3....加上不同的系数,可以表达世界上任何一个函数。
而傅里叶级数则说,我使用sin0x,cos0x,sinx,cosx,sin2x,cos2x,sin3x,cos3x.....加上不同的系数,可以表达世界上任何一个周期函数(甚至非周期)。
于是,傅里叶级数的系数,和某个周期函数式,就形成了对应关系。
知道傅里叶级数的所有系数,就能求出对应的函数式。
知道某个函数式,我可以求出傅里叶级数的所有系数。
这个高等数学下册就会有的。

离散傅里叶变换:
什么是离散傅里叶变换(DFT)?

我用一个非常简单的例子来说明一下。

题目1：

已知y=f(x)=2+3x+x^2，求x=0,1,2时y的值。

题目2:

已知一个二次函数f(x)经过点(0,2),(1,6)(2,12)，求f(x)的函数表达式。

相信这两个题目，人都会做。但是，其中却是蕴含着奥秘。嘿嘿嘿....

在题目1，我们知道f(x)的三个系数,求不同的x时y的值。

而问题2则反过来，我们知道三个不同的x时y的值，反过来求系数。

显然f(x)=2+3x+x^2和二次函数f(x)经过(0,2),(1,6)(2,12),是在表达同一个函数,只是用了不同的形式!

离散傅里叶变换就是这样的,函数式和傅里叶系数是周期函数的两种不同表达形式!

而这里,多项式系数和点坐标是多项式函数的两种不同表达形式！
我们知道多项式系数,可以求点坐标.我们知道点坐标,可以反过来求多项式系数.我们称这种变换为离散傅里叶变换。（点值表示法和系数表示法）

有一点要注意的是,要想通过点坐标求出所有系数是有要求的.

如果是一个二次函数,它有3个系数,那么我们需要三个点坐标的信息,并且这3个点不能重合,即x坐标不能一样.

推广的话,n次函数,有n+1个系数,需要n+1个点的信息.

这在线性代数里能得到有关的叙述和严格的证明。

时域和频域:
傅里叶变换有两个很重要的术语,时域和频域.

其实,在DFT里,简单地说,时域就是点坐标,频域就是系数.

这么说好像不容易理解,毕竟时域和频域是在信号与系统那方面命名的,用在这里有点奇怪.

不过,我们就这样记住好了.时域是点坐标,频域是系数.

要想理解时域和频域的命名含义,可以看有关傅里叶变换的书籍.

快速傅里叶变换:
快速傅里叶变换FFT又是什么呢?我们先看离散傅里叶变换DFT.

既然说FFT是个算法,谈到算法我们当然要谈时间复杂度.要实现DFT,时间复杂度是很高的.

比如从频域到时域(知道系数,给出不同的x坐标,求y坐标),我们有O(n)个坐标要算,每个坐标是O(n)个项,复杂度就O(n2)了.

反过来,从时域到频域(知道坐标,求系数),我们要解一个n元一次方程组.用线性代数的高斯消元法,复杂度是O(n3).

那么,FFT就是来降低这个复杂度的,它能把DFT的这两个转换的复杂度降到O(nlogn)。

多项式乘法与高精度乘法:
实现高精度乘法,其实和实现多项式乘法是几乎一样的.

这话怎说?其实,我们用某个进制去表达一个数字,恰恰和多项式的表达一模一样.

比如15386,它其实是1*10^4+5*10^3+3*10^2+8*10^1+6*10^0.

这和函数f(x)=x^4+5*x^3+3*x^2+8*x^1+6是一个套路.

那么,两个数字的乘法,其实就是两个多项式的乘法.

例如12*43,我们可以理解为(x+2)(4x+3).

(x+2)(4x+3)=4x^2+11x+6,所以12*43=4*10^2+11*10+6.

稍有不同的是,数字乘法要进位,所以(4,11,6)要变成(5,1,6),也就是12*43=516.

所以说，多项式乘法和高精度乘法其实是差不多的。

接下来,看看多项式乘法的算法复杂度.分析很简单,两个for循环,O(n2)。

而FFT告诉你,O(nlogn)就可以实现了。

FFT实现多项式乘法（卷积）:
FFT如何实现多项式乘法,我们先探讨多项式乘法的本质.

多项式乘法,其实就是给出你两个多项式的频域,求它们的积的频域!

我们知道两个多项式的表达式,其实就是知道频域(系数),而我们要求的是它们的积的表达式,也就是求频域.

我们把上面的结论收集起来,说明如何用FFT实现O(nlogn)多项式乘法.
1.FFT能O(nlogn)实现频域到时域
2.FFT能O(nlogn)实现时域到频域
3.我们知道两个多项式的频域,要求它们的积的频域.

只有这三个结论是不够的,我们还需要一个很白痴但你可能想不到的结论:

4.我们知道两个多项式的时域,能O(n)时间求出它们的积的时域.

这个结论看似很难,但讲清楚之后你会发现很简单.

知道两个多项式的时域,就是知道了f(x)经过(x1,f1),(x2,f2)...(xn,fn),g(x)经过(x1,g1),(x2,g2)....(xn,gn)

设k(x)=f(x)*g(x),那么k(x)在x1,x2...xn上的y坐标是多少?一想就知道是f1*g1,f2*g2....fn*gn!!!因为k(xi)=f(xi)*g(xi)啊!!

也就是说k(x)会经过(x1,f1*g1),(x2,f2*g2)....(xn,fn*gn).这些点,就是k(x)的时域啊!

这里的复杂度不用说都知道是O(n),原来"时域相乘"是这么简单的事情,比"频域相乘"简单多了.

现在我们可以用FFT实现O(nlogn)多项式乘法了!

1.用FFT在O(nlogn)时间,把f(x)和g(x)的频域转变为时域。

2.用O(n)时间,把f(x)和g(x)的时域变成k(x)的时域。

3.用FFT在O(nlogn)时间,把k(x)的时域转变为频域。

接下来的FFT的实现..涉及大量的数学知识..搞不来..更玄学了.....

反正FFT就是能够在O(nlogb)时间内将系数表示法转化为点值表示法，相乘后再将点值表示法转为系数表示法，然后实现卷积。

FFT过程：
两个次数界为n的多项式A（x）和B（x）相乘，输入输出均采用系数表示法。（假定n为2的幂）
1）使次数界增加一倍：A（x）和B（x）扩充为次数界为2n的多项式，并构造起系数表示。
2）求值：两次应用2n阶FFT，计算出A（x）和B（x）的长度为2n的点值表示。
3）点乘：计算多项式C(x)=A（x）*B（x）的点值表示。
4）插值：对2n个点值对应用一次FFT计算出其逆DFT，就可以构造出多项式C（x）的系数表示。
第1步和第3步需要时间为O(n)，第2步和第4步运用FFT需要时间为O(nlogn)。

建议去看算法导论里的FFT....或者去看ACdreamer写的两篇博客....

ACdreamers 

ACdreamers

贴个模板吧....

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int N = 500005;
const double pi = acos(-1.0);
char s1[N],s2[N];
int len,res[N];

struct Complex
{
    double r,i;
    Complex(double r=0,double i=0):r(r),i(i) {};
    Complex operator+(const Complex &rhs)
    {
        return Complex(r + rhs.r,i + rhs.i);
    }
    Complex operator-(const Complex &rhs)
    {
        return Complex(r - rhs.r,i - rhs.i);
    }
    Complex operator*(const Complex &rhs)
    {
        return Complex(r*rhs.r - i*rhs.i,i*rhs.r + r*rhs.i);
    }
} va[N],vb[N];

//雷德算法--倒位序  
void rader(Complex F[],int len) //len = 2^M,reverse F[i] with  F[j] j为i二进制反转
{
    int j = len >> 1;
    for(int i = 1;i < len - 1;++i)
    {
        if(i < j) swap(F[i],F[j]);  // reverse
        int k = len>>1;
        while(j>=k)
        {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}
//FFT实现
void FFT(Complex F[],int len,int t)
{
    rader(F,len);
    for(int h=2;h<=len;h<<=1) //分治后计算长度为h的DFT 
    {
        Complex wn(cos(-t*2*pi/h),sin(-t*2*pi/h)); //单位复根e^(2*PI/m)用欧拉公式展开
        for(int j=0;j<len;j+=h)
        {
            Complex E(1,0); //旋转因子
            for(int k=j;k<j+h/2;++k)
            {
                Complex u = F[k];
                Complex v = E*F[k+h/2];
                F[k] = u+v; //蝴蝶合并操作
                F[k+h/2] = u-v;
                E=E*wn; //更新旋转因子
            }
        }
    }
    if(t==-1)   //IDFT
        for(int i=0;i<len;++i)
            F[i].r/=len;
}
//求卷积 
void Conv(Complex a[],Complex b[],int len) //求卷积
{
    FFT(a,len,1);
    FFT(b,len,1);
    for(int i=0;i<len;++i) a[i] = a[i]*b[i];
    FFT(a,len,-1);
}
/*
void init(char *s1,char *s2)
{
    int n1 = strlen(s1),n2 = strlen(s2);
    len = 1;
    while(len < 2*n1 || len < 2*n2) len <<= 1;
    int i;
    for(i=0;i<n1;++i)
    {
        va[i].r = s1[n1-i-1]-'0';
        va[i].i = 0;
    }
    while(i<len)
    {
        va[i].r = va[i].i = 0;
        ++i;
    }
    for(i=0;i<n2;++i)
    {
        vb[i].r = s2[n2-i-1]-'0';
        vb[i].i = 0;
    }
    while(i<len)
    {
        vb[i].r = vb[i].i = 0;
        ++i;
    }
}

void gao()
{
    Conv(va,vb,len);
    memset(res,0,sizeof res);
    for(int i=0;i<len;++i)
    {
        res[i]=va[i].r + 0.5;
    }
    for(int i=0;i<len;++i)
    {
        res[i+1]+=res[i]/10;
        res[i]%=10;
    }
    int high = 0;
    for(int i=len-1;i>=0;--i)
    {
        if(res[i])
        {
            high = i;
            break;
        }
    }
    for(int i=high;i>=0;--i) putchar('0'+res[i]);
    puts("");
}

*/
int main()
{   /*
    while(scanf("%s %s",s1,s2)==2)
    {
        init(s1,s2);
        gao();
    }  */
    return 0;
}