我们所熟知的科学计算一般就是指数值计算,数值计算是计算数学的一个主要部分,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现[1],关于数值计算的研究在计算机被发明之前就已经有了相当的基础,它涉及到的内容包括函数的数值逼近,数值微分与数值积分,非线性方程数值解,数值线性代数,常微分方程与偏微分方程数值解等.数值计算中处理的对象并不仅仅是数值,还包括由数值构成的简单数据结构,例如一般的多项式,无穷级数,矩阵等,数值计算处理问题的一般方法是通过数学推导将问题化归到这些数学对象的运算上.
作为应用数学,数值计算的主要目标是解决来自于生产实践的工程学问题.与此同时,数学工作者做数学研究本身也是一种生产实践,数学研究过程中同样会产生许多问题,与工程学问题不同,这些问题往往只能用抽象的符号来表达,仅用数值计算的方法是不易解决的,对于这类问题解决方案的研究逐渐形成了应用数学的一个新的分支,为了与数值计算相区别,常常称之为符号计算.类似的,我们可以给符号计算下一个简单的定义:符号计算是一门研究用计算机求解各种数学问题的符号计算方法及其理论与软件实现的科学,它是数学(家)的计算数学.符号计算中处理的数据和结果都是符号,这种符号可以是字母,公式,也可以是数.与数值计算不同的是,数是作为一种符号出现在符号计算中的,这就要求关于数的运算应该是绝对精确的,我们接下来就要讨论关于数的高精度运算.
在基于硬件的整数指令中,计算机能够处理的整数是有界的,在目前典型的计算机中整数的溢出界都不超过$2^{64}$,而符号计算中常常需要处理更大的整数,例如阶乘,斐波拉契数列这样简单的数论函数计算,另一个不平凡的例子是所谓的中间表示膨胀[2],例如采用Euclid算法计算两个整系数多项式的最大公约数时,即使输入的两个多项式和输出的最大公约数多项式都具有较小的系数,计算过程中的中间结果仍然很可能出现非常大的系数.设 \begin{align*} F&=7x^7+2x^6-3x^5-3x^3+x+5,\\ G&=9x^5-3x^4-4x^2+7x+7, \end{align*} 在计算过程中将有理数化为整数,我们将得到如下的多项式序列 \begin{align*} &1890x^4-4572x^3-6930x^2-846x+4527,\\ &294168996x^3+257191200x^2-20614662x-142937946,\\ &-103685278369841305200x^2-32576054233115610000x\\ &+122463167842311670000,\\ &2956790833503849546789342057565207098291763520000x\\ &+555325261806247996966034784074025291687620160000,\\ &1092074685733031219201041602791259862659169966184593803518\\ &602418777140682884334769647063543607737698426880000000000, \end{align*} 最后的那个整数达到了118位.除此之外,高精度浮点数的表示和运算也是直接依赖于高精度整数的.
为了简单起见,在下面的讨论中,如果没有明确指出,我们所说的整数都是正整数.关于高精度整数的表示方法,一个简单的想法是选定进制$B$,用整数的$B$进制表示(向量)来表示整数,为了更好的理解这一点,首先对整数的进制表示作一个严格的定义.
即$$(a_n\ldots a_1a_0)_B=\sum_{k=0}^na_k\cdot B^k.$$
高精度整数输入输出时,常常需要做进制转换,即给定正整数$n$的$B$进制表示$n=(a_s \ldots a_1a_0)_B$,需要得到$n$在$B'$进制下的表示.计算方法有如下两种:
当转换极长的数时,最方便的是由转换数字块开始,这些数字块可以用单精度技术来处理,而后用简单的多精度技术把这些块组合在一起[3] .以第一种方法为例,反复地用$B'^m$除它,于是得到$n$的$B'^m$进制表示,然后对于$B'^m$进制表示的每一位通过单精度操作给出$m$位$B'$进制数字.
本章的开头时曾经强调过,数是作为一种符号出现在符号计算中的,这就要求关于数的运算应该是绝对精确的.在机器精度的范围内,现有的计算机可以轻松的完成整数的四则运算,这主要得益于计算机设计师们的工作.而接下来要讨论的高精度整数四则运算,其基本原理其实和基于硬件的整数指令是一致的,事实上,这就是阿拉伯计数法的发明者们很早就发展出的一整套借助纸笔进行的四则运算理论,我们从小就开始学习熟练地利用它们来操作整数.需要特别指出的是,这里将要提到的加减乘除与作为群运算的加减乘除是不同的,因为加减乘除作为群运算是一个抽象的二元映射,而我们现在要做的其实是基于这些群运算和上一节中进制表示所做的"表面"工作,说得具体一些,就是输入两个数的$B$进制表示,输出群运算结果的$B$进制表示,关键点集中在求出$B$进制表示上.
首先来看高精度整数的加法,为了简单起见,假定相加的两个数的位数相同,如果位数不同首先通过添零补齐.注意,这样的假定不会影响算法的效率.
输入:整数$a=(a_s\ldots a_1a_0)_B$,$b=(b_s\ldots b_1b_0)_B$.
输出:$a,b$的和$a+b=(c_{s+1}c_s\ldots c_1c_0)_B$.
定义辅助序列$$d_k=[{a+b\over B^k}]-[{a\over B^k}]-[{b\over B^k}],\quad 0\le k\le s.$$ 按递推公式依次求出$c_k,d_k(0\le k\le s)$
现在我们将这个算法改写成适合于编程的更加紧凑的形式.
输入:两个$n$位整数$u=(u_{n-1}\cdots u_1u_0)_B$,$v=(v_{n-1}\cdots v_1v_0)_B$.
输出:$u,v$的和$w=(w_nw_{n-1}\cdots w_1w_0)_B$,$w_n=0$或$1$.
高精度整数减法的方法与加法完全类似,我们直接给出算法.
输入:整数$a=(a_s\ldots a_1a_0)_B$,$b=(b_s\ldots b_1b_0)_B$.
输出:$a,b$的差$a-b=(c_s\ldots c_1c_0)_B$.
定义辅助序列$$d_k=[{a-b\over B^k}]-[{a\over B^k}]+[{b\over B^k}],\quad 0\le k\le s.$$ 按递推公式依次求出$c_k,d_k(0\le k\le s)$
同样地,也可以将这个算法改写成适合于编程的更加紧凑的形式.
输入:两个$n$位整数$u=(u_{n-1}\cdots u_1u_0)_B$,$v=(v_{n-1}\cdots v_1v_0)_B$,$u\ge v$.
输出:$u,v$的差$w=(w_{n-1}\cdots w_1w_0)_B$.
高精度整数的加减法算法实际上就是对我们熟悉的十进制笔算加减法的精确化描述并推广到$B$进制记数系统后得到的,用$T(n)$表示用它们计算$n$位$B$进制整数和或差时所需要的$B$进制个位数加减法的次数,可以证明$T(n)=O(n)$,这和必要的输入输出操作相比是都是线性级别的,因此有理由认为这就是理想的算法.尽管如此,在后面有关多项式模算术的讨论中我们将看到,作为多项式模算术的一个直接应用,利用有限域上的整数模算术可以把高精度整数加减法分解成在多个处理器上并行的独立任务,这样一来它就会比经典算法要快.
接下来讨论一位数乘多位数的乘法,这里用到的方法和加减法仍然是相似的,同样可以直接给出算法.
输入:整数$a=(a_s\ldots a_1a_0)_B$,$b_0=(b_0)_B$.
输出:$a,b_0$的积$ab_0=(c_{s+1}\ldots c_1c_0)_B$.
定义辅助序列$$d_k=[{ab_0\over B^k}]-b_0[{a\over B^k}],\quad 0\le k\le s.$$ 按递推公式依次求出$c_k,d_k(0\le k\le s)$
同样地,也可以将这个算法改写成适合于编程的更加紧凑的形式.
输入:$m$位整数$u=(u_{m-1}\cdots u_0)_B$,$n$位整数$v=(v_{n-1}\cdots v_0)_B$.
输出:$u,v$的积$w=(w_{m+n-1}\cdots w_1w_0)_B$.
在四则运算的笔算方法中,除法可能是最复杂的,因为除法需要试商,试商包含着猜测的成分,而不能形成有效的算法.机器精度整数除法也要试商,但得益于二进制表示,对它们而言商只有0,1两种情况,只需要比较大小就可以求出.为了有效地计算高精度整数的除法,我们需要将试商的过程算法化[3],首先来看商为一位数的除法,即假定商$\left[\frac{a}{b}\right]$满足$0\le\left[\frac{a}{b}\right]\le B-1$,笔算除法的经验告诉我们,被除数和除数的最高几位数对试商是很重要的,我们常常凭借目测(口算)最高几位数就能基本上把商给确定下来.
而由$$\left[\frac{a_{s+1}B+a_s}{b_s}\right] > \frac{a_{s+1}B+a_s}{b_s}-1$$可以得到 \begin{align*} \left[\frac{a_{s+1}B+a_s}{b_s}\right]+1&\ge\frac{a_{s+1}B+a_s+1}{b_s}\\ &=\frac{a_{s+1}B^{s+1}+(a_s+1)B^s}{b_sB^s}\\ &>\frac{a}{b}, \end{align*} 故$\left[\frac{a}{b}\right]\le\left[\frac{a_{s+1}B+a_s}{b_s}\right]$,证毕. □
记$$q = \min\left\{\left[\frac{a_{s+1}B+a_s}{b_s}\right],B-1\right\},$$则$q$是商$\left[{a\over b}\right]$一个很好的上界,下面将看到它已经很接近商的真值:$q$比商的真值至多大2.
上面只考虑了被除数的头两位与除数的首位,如果允许做更精细的考察,譬如说考虑被除数的头三位与除数的头两位,我们还可以进一步接近商的真值.
输入:整数$a=(a_{s+1}\ldots a_1a_0)_B$,$b=(b_s\ldots b_1b_0)_B$.
输出:$a,b$的商$[{a\over b}]$.
注意到$$r=a-qb=a\bmod b$$在算法的最后一步也同时被算出,以此为基础可以写出一般的商为多位数的除法算法.
输入:整数$a=(a_m\ldots a_1a_0)_B$,$b=(b_n\ldots b_1b_0)_B$.
输出:$a,b$的商$[{a\over b}]=(q_{m-n}\ldots q_1q_0)_B$.
定义辅助序列$$r_k=[{a\over B^k}]\bmod b,\quad m-n+1\ge k\ge 0.$$ 按递推公式依次求出$q_k,r_k(m-n+1\ge k\ge 0)$
高精度整数是计算机代数系统的内置基本类型,四则运算的快慢对系统的性能好坏有着决定性的影响,目前最著名的高精度整数运算库是GNU的GMP[4],许多著名的计算机代数系统都是基于GMP构建的.加法和减法的复杂度关于整数位数是线性的,考虑到输入输出的复杂度关于整数位数也是线性的,因此从算法上来看加减法已经达到了复杂度的下界.高精度除法总可以通过Newton迭代归结为高精度乘法,所以四则运算的主要问题集中在乘法上,设$n$为乘数的位数,就目前已知的情况而言,不同乘法算法的复杂度可以从平凡的$O(n^2)$(普通乘法),$O(n^{\log_2{3}})$(Karatsuba乘法),$O(n^{\log_3{5}})$(Toom-3乘法),$O(n\log^*{n})$(复数域上的FFT)直到$O(n(\log{n})(\log{\log{n}}))$(有限域上的FFT),最后一个已经相当接近线性复杂度了.复杂度较低的算法往往有较大的常数因子,例如当两个相乘数都很小时,普通乘法反而是最快的,实用中常常将这些不同的乘法算法结合起来使用,每次做乘法时都根据相乘两数的大小动态地选择具体采用哪一种算法,选择时用到的参数一般通过实验来确定.
在前面的利用一位数乘多位数来做整数乘法时,已经提到过整数乘法还有更快的算法.在介绍这些高级算法之前,先来分析一下,笔算乘法的算法是如此地直截了当,以至于如果单从整数的角度考虑,应该不可能有更快的算法了.我们现在换一个角度,把目光投向多项式.
对于一个次数界为$n$的多项式$$A(x)=\sum_{k=0}^{n-1}a_kx^k,$$其系数表示法就是一个由系数组成的向量$\mathbf{a}=(a_0,a_1,\ldots,a_{n-1})$,这是我们习以为常的一种表示形式,并且是唯一的;而点值表示法[5]是多项式在$n$个不同点处的点值对构成的集合$$\{(x_0,y_0),(x_1,y_1),\ldots,(x_{n-1},y_{n-1})\},$$其中$y_k=A(x_k),k=0,1,\ldots,n-1$.一个多项式可以有很多种不同的点值表示,这是因为每一组$x_k$都决定着一种点值表示.系数表示与点值表示之间是可以相互转换的,系数表示到点值表示称为多项式多点求值,点值表示到系数表示称为多项式插值,关于这两个主题的高级讨论将留到与多项式相关的章节中,下面关于整数快速乘法的讨论只需要这几个概念就足够了.
已知两个次数界为$n$的多项式的系数表示分别为$$A(x)=\sum_{k=0}^{n-1}a_kx^k,\quad B(x)=\sum_{k=0}^{n-1}b_kx^k,$$设乘积的系数表示为$C(x)=\sum\limits_{k=0}^{2n-2}c_kx^k$,那么$$c_k=\sum_{i+j=k}a_ib_j,$$这样计算乘积时总共需要$n^2$次系数乘法.现在我们首先取定$2n-1$个不同点$x_k$,然后利用多项式多点求值计算出$A(x)$与$B(x)$在这组点上的点值表示,如果$C(x)=A(x)\cdot B(x)$,那么显然有$C(x_k)=A(x_k)\cdot B(x_k)$,只需要$2n-1$次系数乘法就可以计算出乘积$C(x)$的点值表示,再通过多项式插值就可以计算出其系数表示.从以上过程不难看出,利用多项式的点值表示做多项式乘法来得简洁明了,如果多项式的默认表示形式就是点值表示,那么连多项式多点求值和多项式插值这两个步骤也可以省去了,多项式乘法就和向量逐点相乘一样简单,从模方法这个更高的角度来看,多项式的点值表示对应于模取为一次多项式组$x-x_k$的模表示,多项式多点求值和多项式插值则分别对应于中国剩余定理的逆定理与原定理.不幸的是,我们现在讨论的是整数快速乘法,而整数与多项式之间的联系依赖于整数的进制表示$$(a_n\ldots a_1a_0)_B=\sum_{k=0}^na_k\cdot B^k,$$这对应于多项式的系数表示(如果采用整数的模表示,我们也可以让整数对应于多项式的点值表示,不过现在我们不讨论它).这样一来,整数快速乘法就必须解决多项式系数表示与点值表示之间的转换问题,可以把它看成三个步骤:多项式多点求值,向量逐点相乘,多项式插值.
Karatsuba算法[3]是1962年发现的,它使用了一个惊人简单但十分有效的策略.首先研究一次多项式的乘法,设$$A(x)=a_1x+a_0,\quad B(x)=b_1x+b_0,$$选取插值点组为$x_0=-1,x_1=0,x_2=\infty$,则$A(x),B(x)$的点值表示分别为$\{(-1,a_0-a_1),(0,a_0),(\infty,a_1)\}$,$\{(-1,b_0-b_1),(0,b_0),(\infty,b_1\}$,如果$C(x)=A(x)\cdot B(x)$,那么$C(-1)=(a_0-a_1)\cdot(b_0-b_1)$,$C(0)=a_0\cdot b_0$,$C(\infty)=a_1\cdot b_1$,利用简单的多项式插值算法,譬如Lagrange插值公式,可以计算出$C(x)$的系数表示$$(c_0,c_1,c_2)=(C(0),C(0)+C(\infty)-C(-1),C(\infty)),$$即$$C(x)=a_1b_1x^2+(a_0b_0+a_1b_1-(a_0-a_1)(b_0-b_1))x+a_0b_0.$$在这个例子中,多项式多点求值与多项式插值都只用到系数加减法,这表明可以只用三次系数乘法完成一次多项式的乘法.前面这些推理都是多项式的点值表示理论的自然应用,俄罗斯数学家Karatsuba的贡献在于将它应用到了一般的多项式乘法和整数乘法,简单地说来是这样的,如果$a,b$都是$2n$位$B$进制整数,那么在$B^n$进制下$a,b$都是两位数$$a=a_1B^n+a_0,\quad b=b_1B^n+b_0,$$其中$a_i,b_i(i=0,1)$都是$B^n$进制下的个位数,也就是$B$进制下不超过$n$位的数,记$T(n)$是递归地利用Karatsuba乘法将两个n位$B$进制数相乘所需要的$B$进制个位数乘法次数,那么有$T(2n)\le 3T(n)$,由此推出$T(n)=O(n^{\log_2{3}})\approx O(n^{1.59})$.
综合以上这些想法,现在可以写出Karatsuba乘法的详细过程了[2].
输入:$n_a$位整数$a$和$n_b$位整数$b$.
输出:$a,b$的积$c=a\cdot b$.
沿着Karatsuba的思路,考虑二次多项式的乘法,设$$A(x)=a_2x^2+a_1x+a_0,\quad B(x)=b_2x^2+b_1x+b_0,$$选取插值点组为$$x_0=-1,x_1=0,x_2=1,x_3=2,x_4=\infty,$$如果$C(x)=A(x)\cdot B(x)$,类似地,利用简单的多项式插值算法,譬如Lagrange插值公式,可以计算出$C(x) $的系数表示$(c_0,c_1,c_2,c_3,c_4)$,其中 \begin{align*} c_0&=C(0),\\ c_1&=C(-1)+C(0)+C(1)+C(2)+C(\infty),\\ c_2&=C(-1)+C(0)-C(1)-C(2)+C(\infty),\\ c_3&=4C(-1)+C(0)+2C(1)+8C(2)+16C(\infty),\\ c_4&=C(\infty). \end{align*} 如果不计乘以较小常数的乘法和加减法,我们可以只用5次系数乘法完成二次多项式的乘法,以此为基础仿照Karatsuba乘法而得到的算法称为Toom-3乘法.记$T(n)$是递归地利用Toom-3乘法将两个n位$B$进制数相乘所需要的$B$进制个位数乘法次数,那么有$T(3n)\le 5T(n)$,由此推出$T(n)=O(n^{\log_3{5}})\approx O(n^{1.46})$.
更一般地,考虑$r$次多项式的乘法,设$$A(x)=a_rx^r+\cdots+a_0,\quad B(x)=b_rx^r+\cdots+b_0,$$选取插值点组为$$x_0=-r,x_1=-(r-1),\ldots,x_r=0,x_{r+1}=1,\ldots,x_{2r}=r,$$如果$C(x)=A(x)\cdot B(x)$,这一次我们实际地用一下Lagrange插值公式来计算出$C(x) $的系数表示$(c_0,c_1,\ldots,c_{2r})$.
这样得到的一系列快速乘法算法统称为Toom-Cook乘法[3],其中Toom-2乘法与Karatsuba乘法大致相同,只是选取的插值点不一样而已.由于$$\lim_{r\rightarrow\infty}\log_{r+1}{2r+1}=1,$$所以使用Toom-Cook乘法进行$n$位整数乘法所需要的$B$进制个位数乘法次数的下界趋近于$O(n)$,这是一个很好的结果.但是由于多项式多点求值和多项式插值所需要的线性级操作会随着$r$的增加而迅速增加,以至于一般说来仅仅当$r=2,3$时,Toom-Cook乘法才是实用的.
快速傅立叶变换(FFT)[6]被认为是20世纪数值计算和算法领域最重要的成果之一,它常常被用于快速计算卷积,这对于信号处理等工程领域尤其重要.对于符号计算来讲,FFT实际上是多项式快速多点求值和快速插值算法一个高度优化的特例.不过,在这里我们只需要利用多项式乘法和向量卷积之间的相似性,就可以很好地利用FFT来加速多项式乘法和整数乘法.在上一节中,Toom-Cook乘法通过充分地利用分治与递归有效地减少了多项式多点求值和多项式插值操作的消耗;下面将要讨论的FFT乘法则是通过精心地挑选求值点,把系数表示与点值表示之间的转换所需的操作压缩为次线性级别,即$O(n\log{n})$,其中$n$代表相乘的两个整数的位数.
顾名思义,FFT是用来快速计算离散傅立叶变换(DFT)的方法,我们可以将DFT简单的理解成以单位复根作为求值点时多项式系数表示到点值表示之间的转换,关于它的原始定义及详细介绍可以在任何一本有关数值计算的教科书中找到.
为了更好地理解DFT,首先来看一下什么是原根,这是一个代数数论的概念.
从定义容易导出原根的一些特殊性质.
另外一个需要用到的概念是循环卷积.
现在可以正式地给出DFT的定义了.
在下面的定理我们将会发现,原根的良好性质保证了DFT的逆变换和DFT具有相同的结构.
如果采用普通的多项式多点求值方法,譬如说Horner法则,计算$n$阶DFT大约需要$n^2$次系数乘法,为了寻找更快的方法,需要对DFT的结构进行仔细分析.先来看向量$DFT_\omega(f)$的第$j$个分量$f(\omega^j)=\sum\limits_{k=0}^{n-1}f_k\omega^{jk}$,因为$\omega$是$n$次单位原根,所以$\omega^{jk}$至多只有$n$个不同的值$\omega^0,\omega^1,\ldots,\omega^{n-1}$,特别地,如果$j$是$n$的因子,那么$\omega^j$是${n\over j}$次单位根,因此$(\omega^j)^k$至多只有${n\over j}$个不同的值,如果我们已经知道了哪些$\omega^{jk}$具有相同的值,就可以像合并同类项一样将对应于相同的$\omega^{jk}$的系数$f_k$先加在一起,然后再和$\omega^{jk}$相乘,这样就可以有效地减少乘法的次数,跟计算$ab+ac$时转而计算$a(b+c)$是一个道理.设$0\le k_1<k_2<n$,那么$$\omega^{jk_1}=\omega^{jk_2}\Leftrightarrow jk_1\equiv jk_2\mod n,$$即$n|j(k_1-k_2)$,一般地,$n$是素数的整数次幂,不妨设$n=p^m$,其中$p$为素数,设$j$的素因子分解中$p$的幂次为$m'<m$,那么应该有$p^{m-m'}|k_1-k_2$,这可以当作合并$\omega^{jk}$的系数时的依据.广义的快速傅立叶变换可以处理$p$为任意素数的情形,我们这里接下来要讨论的的FFT则特指基$p=2$的快速傅立叶变换.
为了高效地多点求值,FFT方法运用了分治策略,设$n=2^m$次多项式$f$的系数表示向量为$f=(f_0,f_1,\ldots,f_{n-1})$,将它依下标的奇偶性分成两组得到 \begin{align*} f^{[0]}&=(f_0,f_2,\ldots,f_{n-2}),\\ f^{[1]}&=(f_1,f_3,\ldots,f_{n-1}), \end{align*} 那么向量$DFT_\omega(f)$的第$j$个分量$$f(\omega^j)=\sum_{k=0}^{n-1}f_k\omega^{jk}=\sum_{k=0}^{n/2-1}f_{2k}(\omega^{2j})^k+\omega^j\sum_{k=0}^{n/2-1}f_{2k+1}(\omega^{2j})^k,$$注意到$\omega^{j+n/2}=-\omega^j,\omega^{j+n}=\omega^j$,因此可以改写成 \begin{align*} DFT_\omega(f)[j]=DFT_{\omega^2}(f^{[0]})[j]+\omega^jDFT_{\omega^2}(f^{[1]})[j],\\ DFT_\omega(f)[j+n/2]=DFT_{\omega^2}(f^{[0]})[j]-\omega^jDFT_{\omega^2}(f^{[1]})[j], \end{align*} 其中$0\le j<n/2$.如果从多项式的角度来看,这其实就是$f(x)=f^{[0]}(x^2)+xf^{[ 1]}(x^2)$.
经过以上的分析,现在可以写出FFT的详细过程了.
输入:$n-1$次多项式的系数表示向量$f=(f_0,f_1,\ldots,f_{n-1})$,$n$次单位原根$\omega_n$,其中$n=2^m$.
输出:$f$的离散傅立叶变换$DFT_{\omega_n}(f)$.
如果$R=\mathbb{C}$,即使$m$很大,也很容易求出$n=2^m$次单位原根.但是对于有限域,前面已经提到过,如果$q$是素数的整数次幂,那么当且仅当$n|q-1$时有限域$\mathbb{F}_q$才存在$n$次单位原根,当$m$很大时,并非总能轻易地找到这样的素数$q$.尽管如此,相比于复数域而言,使用有限域上的原根做FFT不会涉及到浮点操作,因此也不用考虑中间精度的问题,正因为这样,在利用FFT计算整系数多项式的乘积或者计算高精度整数乘法时,我们还是希望使用有限域.注意到有限域$F_q=\mathbb{Z}/q\mathbb{Z}$上的算术都是以$q$为模的,因此利用$F_q$上的FFT求出的点值向量的各项系数都不应该超过$q$,否则最后无法将结果还原到$\mathbb{Z}$中[6].
设$a=(a_s\ldots a_0)_B$,$b=(b_s\ldots b_0)_B$,那么$a,b$分别对应于多项式$$a(B)=\sum_{k=0}^sa_kB^k,\quad b(B)=\sum_{k=0}^sb_kB^k,$$ $a,b$的乘积$c$对应于多项式$C(B)=\sum\limits_{k=0}^{2s}c_kB^k$,其中$$c_k=\sum_{i+j=k}a_ib_j<\sum_{i+j=k}B^2\le(s+1)\cdot B^2,\quad 0\le k\le 2s.$$在典型的计算机代数系统中,机器精度整数的溢出界$B=2^{64}$,整数位数$s$的界也取为$2^{64}$,那么根据$F_q$的限制应该有$q>(s+1)B^2\approx 2^{192}$.这样大的素数$q$以及有限域$F_q$是不容易使用的,为了解决这个问题,可以借用模方法中的技术,选取三个机器精度的素数$2^{63}\le q_0,q_1,q_2<2^{64}$,分别在$F_{q_i}$中计算出模乘积$a\cdot b\bmod q_i$,然后利用中国剩余定理计算出$\mathbb{Z}$中的乘积$a\cdot b$.
考虑到中国剩余定理的计算,利用有限域上的FFT进行整数乘法的复杂度为$O(n(\log{n})(\log{\log{n}}))$.
[1]数值分析, 清华大学出版社, 2001.
[2]计算机代数, 清华大学出版社, 北京, 2004.
[3]The art of computer programming, volume 2 (3rd ed.): seminumerical algorithms, Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1997.
[4]The GNU MP Bignum Library, http://gmplib.org/.
[5]算法导论, 机械工业出版社, 2006.
[6]Modern Computer Algebra, Cambridge University Press, 2002.