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)的函数表达式.

相信这两个题目,99%的人都会做.但是,我们不是要做这两个题目,而是要寻找它们的奥妙.

在问题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(nlong)实现时域到频域

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能实现nlogn多项式乘法.
现在,我们剩下的问题,也就是最核心的问题是,如何实现FFT算法?

FFT算法包括频域到时域,以及时域到频域.
我们在这里着重讨论频域到时域,也就是上面那个图的第一个箭头.
现在我们知道的是两个函数的频域(n个多项式系数),我们要求它们的时域(n个点坐标).
还有一点是,这个n一定要是2的幂.如果不是,你可以强制改成是.
比如n=6,你可以加上0*n^7+0*n^6,强行让它有8个项.

首先有一个问题,如何选择时域的一组x坐标呢?
我们知道,时域的n个点并不是确定的,它们只需要保证x坐标互不相同就能构成一个时域.
那么,我们是不是可以随便找n个x坐标就可以了呢?像x=1,2,3...n这样,行吗?
答案是不行!前文说道FFT转换有一个限制,这个限制就在这里:我们不能任意选择x坐标!

你要得到x=1,2...对应的y坐标,抱歉,FFT告诉你我做不到,你要求只能老老实实O(n2).
那么,什么样的一组x坐标,FFT能O(nlogn)求出对应的y坐标呢?
答案是x^n=1这个方程的所有根.

(下面所说的需要比较多的数学知识,请看不懂的同志自行百度恶补,这里有一篇参考文章: http://wenku.baidu.com/link?url=RWHc2uPEZWdbo0TTVM3z320hJidSR5hk-6YcYQbG6n9n2WzrDKcx1E4r_XS0ZtIqkLjgxODMq_Q8GlK-wGJgZCyHsd1aPhFjBC9Ute2yB0G)
哈?!这个方程的解不就只有1,如果n是偶数就是1和-1.我可是要n个x坐标喔!
注意了,我这里是说x^n=1的所有复数根,而不只是实数跟.
可能有些同学复数学得少,在这里就困惑了.x^n=1这方程一定有n个复数根吗?
答案是一定的,一元n次方程有且只有n个复数根.这是代数基本定理.
并且,我还可以告诉你这n个复数根分别是什么,使用复数开方公式就可以了.

对于一个复数z,我们使用其三角函数形式表达,z=r(cosθ+isinθ)
那么z开n次方的所有根就是下图:


啥?k∈Z,那这里不就有无限个复数?
并不是哦,注意k是在cos和sin里面的,是有周期的.
并且很容易知道,周期就是n,也就是说k=0和k=n的结果是一样的.
就看式子可能还是很抽象,我们把这n个根画在复平面上,就一目了然了.


这是 算法导论上的图片,它定义了一个符号 .( 下面我使用w[k,n]去表达这个符号)
w[i,n]其实是表示x^n=1这个方程的所有根中,令k=i的那一个(k就是刚才那个式子的k)
事实上w[k,n]的k是k次方的意思,而w[n]则是e^(2πi/n).
但这两个定义是等价的,前者是我自己的理解,后者是算法导论上的定义.(详情参考 算法导论)
在这里,可以很清楚地看到x^8=1中真的有8个复数根,并且那个周期真的是8.

好,我们得到了n个x坐标,但是是复数的.
有人可能会问,这个x坐标是复数的,怎么对应在xy平面上的某个点啊?
答案当然是不能对应,正确来说是不能在xy平面上对应.
这里的x坐标和y坐标,其实已经是一个四维面上的两条轴了,不能单单用xy平面去理解了.
不过,这里FFT并不关心你是实数还是复数,它就是根据x坐标,算出y坐标而已.

那么,FFT算法要怎么O(nlogn)把这n个x坐标转换成对应的y坐标呢?
答案是分治,FFT算法是一个分治的算法,分而治之.
那么,两个关键的问题就是,如何"分"?如何"治"?

考虑如何分,就要先说清楚原问题和子问题是什么,我们现在说清楚.
FFT算法,是已知一个n-1次的多项式函数,求出n个x坐标对应的y坐标,这n个x坐标是x^n=1的n个复数根.(注意是n-1次,不是n次,因为n-1次多项式有n项)
举个例子,现在我有一个(8-1)次的多项式,我要求出x^8=1的8个根作为x坐标,对应的y坐标.
那么子问题是,我有一个(4-1)次的多项式,我要求出x^4=1的4个根作为x坐标,对应的y坐标.

那么,如何从原问题转变成子问题?
首先,我们需要几个挺容易理解的公式.
公式1:w[k,2n]^2=w[k,n]

公式2:w[dk,dn]=w[k,n]

公式3:w[j+k,n]=w[(j+k)%n,n]

公式4:w[-k,n]=1/w[k,n]

现在,以n=4为例子,y=a0+a1*x+a2*x^2+a3*x^3.
我们的问题是,要求出x^4=1的4个根作为x坐标,对应的y坐标.
我们按照x的指数的奇偶性,把y分成两部分ya和yb.
ya=a0+a2*x^2
yb=a1*x+a3*x^3=x(a1+a3*x^2)
然后我们定义y1,y2
y1=a0+a2*x
y2=a1+a3*x

然后,我们得到了子问题,我们已知y1和y2这两个(2-1)次函数.
我要求出x^2=1的2个根作为x坐标,对应的y1和y2.

为什么这样就是子问题了,或者说为什么要选这样的y1和y2作为子问题的多项式?
因为首先它符合子问题的要求,其次是我们有了这两个子问题的答案,能给出原问题的答案.
也就是,它能帮我们实现分而治之的"治".

怎么治?
我们从子问题得到的答案是y1(w[0,2]),y1(w[1,2]),y2(w[0,2])和y2(w[1,2]).
而我们要求的答案是y(w[0,4]),y(w[1,4]),y(w[2,4]),y(w[3,4]).
在这里,上面的几条公式就有用了.

我们看y和y1y2的关系.
y(x)=y1(x^2)+x*y2(x^2)
那么,我们把w[0,4]代入x里,得到:
y(w[0,4])=y1(w[0,4]^2)+w[0,4]*y2(w[0,4]^2)
根据公式1,w[0,4]^2=w[0,2],所以:
y(w[0,4])=y1(w[0,2])+w[0,4]*y2(w[0,2])
我们也代入w[1,4],得到:
y(w[1,4])=y1(w[1,2])+w[1,4]*y2(w[1,2])
那w[2,4]呢?
y(w[2,4])=y1(w[2,2])+w[2,4]*y2(w[2,2])
w[2,2]是啥,根据公式三的周期性质,w[0,2]=w[2,2],所以:
y(w[2,4])=y1(w[0,2])+w[2,4]*y2(w[0,2])
同样道理写出w[3,4]:
y(w[3,4])=y1(w[1,2])+w[3,4]*y2(w[1,2])

写在一起就是:
y(w[0,4])=y1(w[0,2])+w[0,4]*y2(w[0,2])
y(w[1,4])=y1(w[1,2])+w[1,4]*y2(w[1,2])
y(w[2,4])=y1(w[2,2])+w[2,4]*y2(w[2,2])
y(w[3,4])=y1(w[1,2])+w[3,4]*y2(w[1,2])

有没有发现,求y(w[0,4]),y(w[1,4]),y(w[2,4]),y(w[3,4])的式子中,除了w[0,4],w[1,4],w[2,4],w[3,4],其他项都是子问题返回的答案,y1(w[0,2]),y1(w[1,2]),y2(w[0,2])和y2(w[1,2]).
那剩下的问题是怎么求w[0,4],w[1,4],w[2,4],w[3,4]了.
刚开始,我以为这些根是用某种特殊形式表现的.
结果发现,它就直接用浮点复数表现了!
也就是说,即使你想实现整数多项式相乘,使用FFT得到的结果是一个浮点数多项式,你需要四舍五入去得到正确的整数!

那w[0,4],w[1,4],w[2,4],w[3,4]到底怎么求?你可以直接用开方的那个式子!
那式子中要开方的复数是1,所以r=1,θ=0即可.因此:
w[k,n]=cos(2kπ/n)+isin(2kπ/n)
将k和n都代进去算的话就得到:
w[0,4]=1
w[1,4]=i
w[2,4]=-1
w[3,4]=-i
这里n=4的情况下,这些复数的实部和虚部都是整数.
但从8开始就不是整数的了,甚至大部分都是无理数,因此只能用浮点数去表示.

还有一个小问题是递归的终点,终点就是n=1的时候.
这个时候,我们要求的是w[0,1]其实就是1作为x坐标,对应的y坐标.
在这里我们的多项式只有一项,也就是只有常数,因此我们返回这个常数就行.

到这里,我们真的实现了从频域到时域的转换.
不难知道这个算法复杂度是O(nlogn).(f(n)=2*f(n/2)+O(n),和归并排序类似.)
我们用python去实现一下,这个优美的FFT算法.

代码:

运行结果:

代码解释:

多项式我们使用列表来储存,从前到后分别是0次项,1次项,2次项...
子问题返回的n个y坐标也是用列表储存,从前到后分别是y0,y1,y2...

python中复数能直接表示和运算,无需特殊的类.
并且,虚数i在python中用j表示.

第6,7行是递归终点的特殊处理.
第9到13行是构建y1和y2.
第15,16行是递归调用,得到子问题答案.
第19行到第23行是计算并返回原问题答案.
第20行和第21行的计算公式看起来比较复杂,但根据前文推导就发现也不是太难.

我们的样例是f(x)=2+3x+x^2+2x^3,也就是频域是[2,3,1,2]
从运行结果可以看到,FFT得到的时域是一堆浮点数.
我们整理后得到的时域是[8,1+i,-2,1-i]
如果你将x^4=1的四个根1,i,-1,-i代入多项式,得到的结果也就是这四个值.
这说明,算法正确!

我们已经能实现从频域到时域了,那我们要怎么做才能实现时域到频域?
数学家发现,从频域到时域和从时域到频域的实现几乎一样!
这要使用线性代数的知识去解释.

从频域到时域,其实就是频域的n个系数组成一个列,然后乘一个n*n矩阵,得到另一个列,而这个列就是y坐标组成的列.如下图(出自 算法导论):


为什么中间的n*n矩阵是这个样子?你用矩阵乘法算一算某个行乘列a就懂了,这其实就是把w[k,n]代进了多项式后的式子.
为什么我们要把从频域到时域理解成矩阵乘法呢?因为矩阵乘法有一个特殊的性质.
如果y=M*a,矩阵M是可逆的,那么a=M^(-1)*y.
它的意思就是,如果M是可逆矩阵,那么我们拿M的逆去乘时域y,就能得到频域a.
也就是,它能实现从时域到频域.

那M可逆吗?这个矩阵其实是有名字的,叫 范德蒙矩阵.
嗯,在线性代数里可以证明,它是可逆的.也就是说,有戏!
那M的逆是什么呢?这个线性代数也帮我们算好了,如下图:


有了逆矩阵有什么用呢?矩阵乘法复杂度可是O(n3)的.
谁说用矩阵乘法了?我们还是用FFT.
我们把得到的时域y理解为新的一组系数,把要求的频域a理解为新的一组y坐标,然后新选取的一组x坐标是w[-1,n],w[-2,n]...
你就会发现,这样FFT转换的含义刚好和这个矩阵乘法式子吻合.并且选取的这组x坐标,也能像刚才那组一样分而治之(你可以像上面一样试一下.)
也就是说,我们可以依样画葫芦,使用FFT实现从时域到频域,只是选取的x坐标换了,并且最后得到的结果要除以一个n.
w[-1,n],w[-2,n]怎么求...公式四w[-k,n]=1/w[k,n],完.

代码:

运行结果:

代码解释:

我们在原来的FFT函数里加上inverse参数,inverse为true的话,表示我们的FFT要实现从时域到频域.
然后在第21行加上一句,如果inverse为true则w=1/w,即现在的w表示w[-i,n]而不是w[i,n].

样例还是使用原来的,我们先将频域转为时域得到y,再从时域转为频域得到new_fx.
由于FFT函数里没有帮我们将结果除以n,所以我们在第30行自己实现了.

整理一下运行结果后可以看到,new_fx=[2,3,1,2].
也就是说,new_fx=fx!我们的FFT转换是对的!
我们实现了从频域到时域,也实现了从时域到频域.

FFT算法实现多项式乘法:

最后我们来用FFT算法实现多项式乘法,在这里我们要特殊处理一些东西.
一个是多项式的长度,不是2^k要变成2^k.
另一个是积的时域的问题,假如两个要乘的多项式都是有n项,我们在FFT转换后得到n个点.
这n个点的y坐标相乘得到目标多项式的n个点,但是目标多项式理论上最多是有2n-1项的(如果是整数乘法则可以有2n位).
所以,这n个点会造成信息的丢失,它并不能构成目标的时域.
怎么办?我们一开始求2n个点就好了!

也就是把要乘的那两个多项式的长度增长到2n,增长的部分的系数全部填0,OK~!
那假如两个要乘的多项式的项数不一样呢?简单,把短的拉成和长的一样.
(PS:这里的项数的意思其实是最高次数+1,例如x^2+1我们要说成3项.)

代码:

运行结果:

代码解释:

上图代码调用的fft函数就是前一份代码里的fft函数.

代码里的27行到34行,都是在处理长度问题,包括长度统一,找符合条件的最小的2的幂,以及扩大两倍.

36,37行是第一个箭头,实现频域到时域.

39行则是第二个箭头,求出目标多项式的时域.

41,42行是第三个箭头,实现时域到频域.

我们用(1+2x+3x^2)*(4+5x)作为样例,手算一下答案是4+13x+22x^2+15x^3.

那运行结果呢,一堆很难看的浮点数.

我们整理一下,得出结果就是4+13x+22x^2+15x^3.

虽然结果有8项,但很容易发现后面4项都是0.

FFT实现多项式乘法,在实际的操作上有许多种版本,我这里只是其中的一种,并且是效率不高的一种.

但是,只要思路掌握了,无论什么样的版本都是一个套路.