2013年高中数学论文 图形计算器应用能力测试活动学生 用快速傅立叶变换实现多项式乘法
2013年高中数学论文 图形计算器应用能力测试活动学生 用快速傅
立叶变换实现多项式乘法
在CLASSPAD330中，main 模块中Interactive 栏目Transformation 中的expand 一直是广大高中生喜闻乐见的好功能。用法简洁易懂。
但是，当你需要算两个很长很长的多项式的乘积的时候，我相信你一定会被冗长的表达式输入问题伤透脑筋（我*，怎么要输这么多的“x m ”……各种残念）。今天，我就带大家认识一种高贵冷艳的计算多项式乘法的方法。
在进入正式介绍之前，我们先要普及一些前提知识。那就是FFT 。
FFT 。俗称“法法塔”。学名快速傅立叶变换。FFT 有什么用呢？它可以在O(nlogn)的时间复杂度内完成一个函数对n 个单位虚根的求值，比朴素的O(n 2)要厉害得多。下面的几段是关于FFT 具体操作过程的介绍，如果没有兴趣可以跳过，因为对我们的实际操作没有太大影响,CLASSPAD330将会为我们完成它的具体工作的。
单位根定义n 次单位根为满足w n =1的复数w 。恰好有n 个单位根e 2πik/n (k=0, 1, ...,n-1)，
其中i 是虚数单位(i 2=-1)，而e iu = cos u + i sin u
定义w n =e 2πik/n ，称为n 次主单位根(principal nth root of unity)，则所有其他单位
根都是它的幂。这样，n 个单位根可以写成w n 0,w n 1,...,w n n-1 ，它们在乘法下构成群，且结构和加法群Zn 一样。容易看出，w n k =-w n k+n/2。
设多项式为∑-==10)(n j j j
x a x A ，那么我们将奇数项和偶数项分开处理。
A [0](x)=a 0+a 2x+...+a n-2x n/2-1
A [0](x)=a 1+a 3x+...+a n-1x n/2-1
那么A(x) = A [0](x 2) + xA [1](x 2)。前面提到过w n k =-w n k+n/2，因此(w n k )2=(w n k+n/2)2这样，我们可以同时计算A(x)在这两个点的值。这样，计算A(x)在n 个点上值的任务被转化成了计算A [0](x)和A [1](x)在n/2个点上值的任务了。
上式的操作实际上是用y [0]k 和y [1]k 算y k 和y k+n/2，这样的操作称为蝴蝶操作(butterfly operation) ，如下图：
这样，只需要花O(n)的时间就可以把序列y [0]和y [1]
合并成序列y 了。递归方程为
T(n)=2T(n/2)+O(n)，因此解为T(n)=O(nlogn)。远低于O(n 2)。实际操作的时候可以用迭代
来完成，这样就可以避开傲娇卖萌极慢无比的系统栈了。唯一可惜的是CLASSPAD330上的program 无法使用数组，所以暂时无法在CLASSPAD330上手动实现一个FFT （包括后面要讲的IFFT ）。读者可以自己用C++或其他高级语言在电脑上实现一下。难度应该不大。
为了能够更好的解决问题我们还需要IFFT ，即FFT 的逆操作。因为是逆操作，所以我们可以不加思索地推出只需要把上面计算过程中的w 换为-w 即可。 Input:Expand((x+2)(x+3))
Output:x 2+5·x+6
接下来让我们切入正题。也许你会觉得这个高贵冷艳的FFT 和我们要用的多项式乘法八
竿子打不着。但是它的确可以在实践中加速多项式乘法的操作，从O(n 2)到O(nlogn)。
对于一个多项式∑-==10)(n j j j
x a x A ，我们有两种表示法。系数表示法和点值表示法。
系数表示法就是记录系数向量a=(a 0,a 1,...,a n-1)。注意要用列向量。这样若两个多项式的系数向量为a,b 。那么多项式加法就相当于a+b （这应该是非常直观的），而多项式乘法就对应于向量卷积a ?b 。
点值表示法是一个点-值对的集合{(x 0,y 0),(x 1,y 1),...,(x n-1,y n-1)}，其中所有xi 两两不相等且y k =A(x k )。这就意味着一个多项式有无穷多种点值表示法。但是如果限定多项式的次数小于n ，那么点值表示法和多项式就一定是一一对应的了。
那么我们可以这样进行多项式乘法的计算：
1.将两个多项式的系数表示法高位补0，得到两个次数为2n 的多项式。
2.将两个多项式的系数表示法转换为点值表示法。
3.把两个点值对应的对应点的数值即y k 相乘得到C 即结果的点值表示法。
4.再把C 的点值表示法转换为系数表示法。
这也是应用FFT 进行多项式计算的原理。
为了应用FFT ，我们不妨将点都取单位虚根，这样第2步就变成了FFT ，第4步也就编程了IFFT 。
那么我们重新表述一下上述过程 1.将两个多项式的系数表示法高位补0，得到两个次数为2n 的多项式。
2.把2n 次单位根作为求值点，用FFT 计算两个多项式的点值表示法。
3.把两个点值对应的对应点的数值即y k 相乘得到C 即结果的点值表示法。
4.用IFFT 把C 转换为系数表示法。
接下来我们用CLASSPAD330来实现一下上述过程。
首先我们引入一个具体问题。
A(x)=1+2x+6x 2
B(x)=3+5x+9x 2+7x 3
求两个多项式的乘积。
则系数向量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)。则
第二步计算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}
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和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