基于多项式快速傅立叶变换的大整数乘法

【关键字】

多项式乘法 快速傅立叶变换(FFT) 离散傅立叶变换(DFT) 大整数乘法 分治策略 ACM 高精度算法库
【摘要】

大整数的高精度运算是一个很重要的算法课题。这篇文章讨论了如何实现大整数的快速四则运算问题。目的是为了建立一个大整数算法库。主要讨论了如何实现 乘法。
【正文】

一：引言

大整数的四则运算中，加法和减法没有太多的优化余地；对于除法，则应充分用到乘法的优势，这要通过巧妙的运用牛顿迭代算法。因而乘法运算是大整数四则运算中的核心内容。进行乘法运算时，如果我们模拟平时进行整数乘法的步骤复杂度将达到O(n^2)，应用快速傅立叶变换可以将复杂度降到O(n*lg(n)) ，从而大大提高程序的运行效率并扩大了使用范围。本文先从理论角度介绍了大整数乘法运算的FFT(Fast Fourier Transforms)算法。并从编程角度讨论了一些实现时应注意的边界和精度问题。（针对ACM/ICPC程序设计比赛设计，因而某些地方为了提高效率而牺牲了一部分其它性能）

二：大整数的表示方法

一个N进制的整数可以表示成SUM(ai*(N^i))。也就是说整数运算问题可以转换为多项式的运算问题。例如：对N进制的整数X=P(N)和Y=Q(N)有XY＝R(N)。（其中P(N)和Q(N)是以N为变量的n-1次多项式X和Y的各位是其系数，且P(N)*Q(N)=R(N)）

所以下面的讨论就针对多项式的乘法运算展开。

我们都比较熟悉多项式的这种以系数ai来表示的形式。但我们今天更敢兴趣的是它的另一种表示形式——点值表示法。由n个不同的变量x及其对应的多项式取值y所构成的n点可以唯一确定一个n－1次多项式。这也就形成了在这种联系基础上的另一种多项式表示方法：用这些点(x,y)表示多项式的形式。

三：在多项式系数表示法的基础上做乘法运算

如果我们要把两个多项式乘起来。用系数表示法如右图：我们需要进行O(n^2)次基本乘法运算。这也是我们所手工计算是经常采用的方法。这种方法看上去已经没有什么可以优化的了。因为它看似没有做任何多余的事情。其中的每一组乘法运算都会对最终的结果产生影响。而且我们不能通过调整计算顺序或使用什么技巧降低算法的时间复杂度。最多只能通过一些语法和其它编程技巧使得其常数因子稍微小一些。

四：在点值表示法的基础上做乘法运算

对于两个多项式Y1(x),Y2(x).他们的一组点值表示法形式分别为(x[0],y1[0]),(x[1],y1[1])...(x[n-1], y1[n-1])和(x[0],y2[0]),(x[1],y2[1])...(x[n-1], y2[n-1]) 在这种形式的基础上做乘法运算就可以用O(n*lg(n))的复杂度得到多项式的成绩。因为Y(x)=Y1(x)*Y2(x)在x[i]处的值是y1[i]*y2[i] 。同理可求出其它点的值。

要把n－1次多项式和m－1次多项式做成绩，我们只需要在n＋m－1个不同的的x点分别求处对应的y1，y2值。该点所对应的y值就为y1×y2。只需做 这样的运算就可以得到乘积多项式的点值表示形式。

现在，速度瓶颈已经转移到如何在点值表示和系数表示之间进行快速转换。如果任给n个x值算y值从而得到n个点的话时间复杂度为 。把这n个点值表示法的点对应成为n－1次多项式也得花费相同复杂度的时间。那么我们的努力也就徒劳。怎么找到更好的转换算法就成了解决大整数乘法的关键所在。

五：用快速傅立叶变换和傅氏逆变换对多项式的两种形式进行 的转换

我们利用单位虚根之间的内部联系通过在这些特殊点上进行抽样(取值)以降低采样的时间复杂度。为了后面使用分治策略解决这个问题我们不失一般性的假设多项式的次数为n－1且n为2的幂(低次多项式可以看作高次系数为0的高次多项式)

六：利用FFT的大整数乘法

只要将多项式中的变量看成一个固定的常数N。那么各项系数就对应着N进制各位的数字。所以我们进行大整数相乘的方法如下：

1：把一个乘数(整数)的个,十,百,千。。。位分别放入a数组中的第0,1,2,3…项

2：把另一个乘数按同样的方法对应到b数组中。

3：对于a，b数组分别求FFT得到虚数数组A，B。（按n位处理，n为2的幂，不够的话高位补0）

4：虚数数组C中C[i]=A[i]*B[i]。（该乘法为虚数相乘需要自己实现运算符重载）

5：对C做傅立叶逆变换得到c。

（本来是虚数组，但由于整数乘法这个实际问题的解一定是整数故最终c是整数数组，虚部舍掉，实部四舍五入）

6：这个c即对应着大整数成绩的结果

七：复杂度问题

每次运算时a被分为两个和原来相同大小的数组，经FFT后在以O(n)的速度进行重组得到A，整个傅立叶变换（FFT）和傅立叶逆变换的时间复杂度如为:O(n*lg(n)).

八：实现过程中应该考虑的问题

1：整数溢出问题

由于傅立叶变换法类似于先做不进位乘法，然后处理进位问题。故当两个数位数最大(出于现有计算机计算速度等方面的原因，这里考虑1,000,000以内的数 )且都为N－1时，中间结果最有可能溢出。在C＋＋语言中，N<=43时中间结果不超出int型，N<=61时中间结果不会超过unsigned型。

2：精度问题

由于只进行log(n)次迭代，而double型有18位有效数字，因而误差远远小于0.5。因而不会对运算结果产生影响。

3：空间问题

由于处理位数比较大，并且存储效率比较低（每一位占4个字节）占用内存间比较大。因而当位数接近处理上限（这里考虑1,000,000）时，局部变量已经不能使用了。只能用new去申请空间。否则会在运行时出错。

注：源码可在http://victorcheng.wy8.net下载

【参考文献】
[1]:Thomas H.Cormen, Charles E.Leiserson, Ronald L.Rivest and Clifford Stein. Introduction to Algorithm. The MIT Press, second edition.

[2]:李庆阳，王能超，易大义，《数值分析》