辽宁省沈阳市第十五中学2013年高中数学论文 图形计算器应用能
力测试活动学生 用快速傅立叶变换实现多项式乘法
宁波效实中学 1502班 周杨
在CLASSPAD330中,main模块中Interactive栏目Input:Expand((x+2)(x+3)) Transformation中的expand一直是广大高中生喜闻乐见Output:x2+5·x+6 的好功能。用法简洁易懂。
但是,当你需要算两个很长很长的多项式的乘积的时候,我相信你一定会被冗长的表达式输入问题伤透脑筋(我*,怎么要输这么多的“xm”……各种残念)。今天,我就带大家认识一种高贵冷艳的计算多项式乘法的方法。
在进入正式介绍之前,我们先要普及一些前提知识。那就是FFT。 FFT。俗称“法法塔”。学名快速傅立叶变换。FFT有什么用呢?它可以在O(nlogn)的时间复杂度内完成一个函数对n个单位虚根的求值,比朴素的O(n2)要厉害得多。下面的几段是关于FFT具体操作过程的介绍,如果没有兴趣可以跳过,因为对我们的实际操作没有太大影响,CLASSPAD330将会为我们完成它的具体工作的。
π
单位根定义n次单位根为满足wn=1的复数w。恰好有n个单位根e2ik/n(k=0, 1, ...,n-1),其中i是虚数单位(i2=-1),而eiu = cos u + i sin u
π
定义wn=e2ik/n,称为n次主单位根(principal nth root of unity),则所有其他单位
根都是它的幂。这样,n个单位根可以写成wn0,wn1,...,wnn-1 ,它们在乘法下构成群,且结构和加法群Zn一样。容易看出,wnk=-wnk+n/2。
设多项式为A(x)??axjj?0n?1j,那么我们将奇数项和偶数项分开处理。
A[0](x)=a0+a2x+...+an-2xn/2-1 A[0](x)=a1+a3x+...+an-1xn/2-1
那么A(x) = A[0](x2) + xA[1](x2)。前面提到过wnk=-wnk+n/2,因此(wnk)2=(wnk+n/2)2这样,我们可以同时计算A(x)在这两个点的值。这样,计算A(x)在n个点上值的任务被转化成了计算A[0](x)和A[1](x)在n/2个点上值的任务了。
上式的操作实际上是用y[0]k和y[1]k 算yk和yk+n/2,这样的操作称为蝴蝶操作(butterfly operation) ,如下图:
这样,只需要花O(n)的时间就可以把序列y和y合并成序列y了。递归方程为
2
T(n)=2T(n/2)+O(n),因此解为T(n)=O(nlogn)。远低于O(n)。实际操作的时候可以用迭代来完成,这样就可以避开傲娇卖萌极慢无比的系统栈了。唯一可惜的是CLASSPAD330上的program无法使用数组,所以暂时无法在CLASSPAD330上手动实现一个FFT(包括后面要讲的IFFT)。读者可以自己用C++或其他高级语言在电脑上实现一下。难度应该不大。
为了能够更好的解决问题我们还需要IFFT,即FFT的逆操作。因为是逆操作,所以我们可以不加思索地推出只需要把上面计算过程中的w换为-w即可。
1
[0][1]
接下来让我们切入正题。也许你会觉得这个高贵冷艳的FFT和我们要用的多项式乘法八
2
竿子打不着。但是它的确可以在实践中加速多项式乘法的操作,从O(n)到O(nlogn)。
对于一个多项式A(x)??axjj?0n?1j,我们有两种表示法。系数表示法和点值表示法。
系数表示法就是记录系数向量a=(a0,a1,...,an-1)。注意要用列向量。这样若两个多项式的系数向量为a,b。那么多项式加法就相当于a+b(这应该是非常直观的),而多项式乘法就对应于向量卷积a?b。
点值表示法是一个点-值对的集合{(x0,y0),(x1,y1),...,(xn-1,yn-1)},其中所有xi两两不相等且yk=A(xk)。这就意味着一个多项式有无穷多种点值表示法。但是如果限定多项式的次数小于n,那么点值表示法和多项式就一定是一一对应的了。
那么我们可以这样进行多项式乘法的计算:
1.将两个多项式的系数表示法高位补0,得到两个次数为2n的多项式。 2.将两个多项式的系数表示法转换为点值表示法。
3.把两个点值对应的对应点的数值即yk相乘得到C即结果的点值表示法。 4.再把C的点值表示法转换为系数表示法。 这也是应用FFT进行多项式计算的原理。
为了应用FFT,我们不妨将点都取单位虚根,这样第2步就变成了FFT,第4步也就编程了IFFT。
那么我们重新表述一下上述过程
1.将两个多项式的系数表示法高位补0,得到两个次数为2n的多项式。 2.把2n次单位根作为求值点,用FFT计算两个多项式的点值表示法。 3.把两个点值对应的对应点的数值即yk相乘得到C即结果的点值表示法。 4.用IFFT把C转换为系数表示法。
接下来我们用CLASSPAD330来实现一下上述过程。 首先我们引入一个具体问题。 A(x)=1+2x+6x2
B(x)=3+5x+9x2+7x3
求两个多项式的乘积。
则系数向量a=(1,2,6),b=(3,5,9,7)。
第一步添0,则a’=(1,2,6,3,0,0,0,0,0),b’=(3,5,9,7,0,0,0,0,0)。则
Input:a:={1,2,6,0,0,0,0,0} {1,2,6,0,0,0,0,0} Input:b:={3,5,9,7,0,0,0,0} {3,5,9,7,0,0,0,0}
第二步计算FFT。不过在计算前要记得把模式调为Cplx(虚数)模式。 Input:x:=FFT(a) Input:y:=FFT(b) 第三步求对应点乘。 Input:c:=x*y 第四步求IFFT Input:z:=IFFT(c) Output: {3,11,37,55-1.25e-12i,68,42,0,-1.25e-11+1.25e-12i}
2
这个答案看上去似乎很奇怪。但仔细观察你会发现后面跟的那个虚数几乎可以忽略不计。其实这是在FFT和IFFT计算过程中产生的精度问题。一般情况下我们可以直接忽略后面的虚数。那么得到的就是乘积的系数表示,即c=(3,11,37,55,68,42).
但实际上,这个精度的问题也是可以避免的。只需要将第二步到第四步写在一个表达式里面就可以了。如下所示: Input:IFFT(FFT(a)*FFT(b)) Output: {3,11,37,55,68,42,0,0}
最后,附上一张模拟器运行图。
最后我们验证一下。记得先用Edit栏里的Clear All Variables清空一下变量。 Input:expand((1+2x+6x^2)*(3+5x+7x^2+7x^3)) Output:42*x^5+68*x^4+55*x^3+37*x^2+11*x+3
验证成功。
这个方法在高端洋气之中透出一点2b气息,但是我觉得还是不失为一项实用的技术。欢迎读者对此方法提出批评建议。
PS:本文中logn默认为log2n
3