高精度运算精讲-优快云博客

高精度运算(转)

整数

整数表示进制转换四则运算

一般原理 Karatsuba乘法 Toom-Cook乘法 FFT乘法

我们所熟知的科学计算一般就是指数值计算,数值计算是计算数学的一个主要部分,它研究用计算机求解各种数学问题的数值计算方法及其理论与软件实现[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$ 进制表示(向量)来表示整数,为了更好的理解这一点,首先对整数的进制表示作一个严格的定义.

定义1 设 $B/in/mathbb{N}_+$ , $B/ge 2$ ,称记号 $(a_n/ldots a_1a_0)_B$ 为 $B$ 进制整数,如果 $a_k/in/{0,1,/ldots,B-1/},a_n/ne 0,0/le k/le n.$ $(a_n/ldots a_1a_0)_B$ 具有的整数值由如下规则确定:

$(a_0)_B = a_0$ ;
$(a_n/ldots a_1a_0)_B = B/cdot(a_n/ldots a_1)_B+a_0$ .

即 $(a_n/ldots a_1a_0)_B=/sum_{k=0}^na_k/cdot B^k.$

事实上,给定进制 $B$ 后,整数的进制表示可以看成由有限集 $/{0,1,/ldots,B-1/}$ 的直和之并形成的无穷维向量空间(类似于一维多项式空间)到整数集的同构.也就是说,每个整数都可以唯一地表示成一个 $B$ 进制整数,这就是进制转换的基础,具体的计算公式见下面的定理.

定理1 给定整数 $n$ ,存在非负整数 $m$ 使得 $B^m /le n < B^{m+1}$ ,并且存在唯一的 $m$ 位 $B$ 进制整数满足 $(a_m /ldots a_1a_0)_B=n$ ,其中 $a_k=[{n/over B^k}]/bmod B,/quad k=0,1,/ldots,m.$

从计算机科学的角度来看,整数的进制表示包含着数据结构的思想,它将高精度整数所包含的信息用一列较小的单元来共同存储和表示.在现有的符号计算软件中,高精度整数典型的存储格式[3] 如下: $/begin{center} /begin{tabular}{|c|p{8em}|} /hline $/pm e$ & /hfil $d$ /hfil// /hline /end{tabular} /end{center}$ 其中 $d$ 就是整数的 $B$ 进制表示(向量), $e$ 是 $B$ 进制表示的位数,即这个整数在 $B$ 进制下是多少位数.进制 $B$ 在软件设计之初就要确定下来,高精度整数四则运算的时间随着整数的位数增加而增加,所以在不超过机器精度整数溢出界的范围内,应该选取尽量大的基数 $B$ ,以便充分利用基于硬件的整数指令.在现有的软件中,常常取 $B$ 为10或2的方幂,前者不必在输入输出时进行进制转换,适用于对输入输出速度要求特别高的场合;后者用于计算时有着天然的优势,在后面将会看到,当进制 $B$ 为2的方幂时,很容易通过移位完成一部分乘除法运算.

进制转换

高精度整数输入输出时,常常需要做进制转换,即给定正整数 $n$ 的 $B$ 进制表示 $n=(a_s /ldots a_1a_0)_B$ ,需要得到 $n$ 在 $B'$ 进制下的表示.计算方法有如下两种:

在 $B$ 进制下除以 $B'$ .设 $n$ 的 $B'$ 进制表示为 $(b_t /ldots b_1b_0)_{B'}$ ,根据进制表示的定义,在 $B$ 进制下依次计算出 $b_0=n/bmod B',b_1=[{n/over B'}]/mod B',/ldots.$
在 $B'$ 进制下乘以 $B$ .设 $B'$ 进制整数 $x=0$ ,利用Horner法则,在 $B'$ 进制下依次计算 $x=Bx+a_k$ , $k$ 从 $s$ 到0.最后得到 $x$ 即为 $n$ 在 $B'$ 进制下的表示.

当转换极长的数时,最方便的是由转换数字块开始,这些数字块可以用单精度技术来处理,而后用简单的多精度技术把这些块组合在一起[3] .以第一种方法为例,反复地用 $B'^m$ 除它,于是得到 $n$ 的 $B'^m$ 进制表示,然后对于 $B'^m$ 进制表示的每一位通过单精度操作给出 $m$ 位 $B'$ 进制数字.

四则运算

本章的开头时曾经强调过,数是作为一种符号出现在符号计算中的,这就要求关于数的运算应该是绝对精确的.在机器精度的范围内,现有的计算机可以轻松的完成整数的四则运算,这主要得益于计算机设计师们的工作.而接下来要讨论的高精度整数四则运算,其基本原理其实和基于硬件的整数指令是一致的,事实上,这就是阿拉伯计数法的发明者们很早就发展出的一整套借助纸笔进行的四则运算理论,我们从小就开始学习熟练地利用它们来操作整数.需要特别指出的是,这里将要提到的加减乘除与作为群运算的加减乘除是不同的,因为加减乘除作为群运算是一个抽象的二元映射,而我们现在要做的其实是基于这些群运算和上一节中进制表示所做的"表面"工作,说得具体一些,就是输入两个数的 $B$ 进制表示,输出群运算结果的 $B$ 进制表示,关键点集中在求出 $B$ 进制表示上.

整数加法

首先来看高精度整数的加法,为了简单起见,假定相加的两个数的位数相同,如果位数不同首先通过添零补齐.注意,这样的假定不会影响算法的效率.

算法1 (整数加法)

输入:整数 $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)$

$c_k=(a_k+b_k+d_k)/bmod B$
$d_{k+1}=[{a_k+b_k+d_k/over B}]$

容易看出,这里定义的辅助序列 $d_k$ 其实就是加到第 $k$ 位的进位,由于不用输出各位上的进位,算法实际实现时可以只用一个变量来存储 $d_k$ .尽管我们对这种笔式加法早已烂熟于胸,但是算法的有效性还是严格地证明一下比较好.

证明根据进制表示和 $d_k$ 的定义,我们有 $/begin{align*} a_k+b_k+d_k&=([{a/over B^k}]-B[{a/over B^{k+1}}])+([{b/over B^k}]-B[{b/over B^{k+1}}])// &+([{a+b/over B^k}]-[{a/over B^k}]-[{b/over B^k}])// &=[{a+b/over B^k}]-B[{a/over B^{k+1}}]-B[{b/over B^{k+1}}], /end{align*}$ 由此推出 $/begin{align*} (a_k+b_k+d_k)/bmod B&=[{a+b/over B^k}]/mod B=c_k,// [{a_k+b_k+d_k/over B}]&=[{a+b/over B^{k+1}}]-[{a/over B^{k+1}}]-[{b/over B^{k+1}}]=d_{k+1}. /end{align*}$ 证毕. □

现在我们将这个算法改写成适合于编程的更加紧凑的形式.

算法2 (整数加法)

输入:两个 $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$ .

令 $j/leftarrow 0,k/leftarrow 0$ .
令 $w_j/leftarrow (u_j+v_j+k)/bmod B$ , $k/leftarrow (u_j+v_j+k/ge B)$ .
$j/leftarrow j+1$ .若 $j<n$ ,转到 $2$ ;否则令 $w_n/leftarrow k$ ,返回 $w$ .

整数减法

高精度整数减法的方法与加法完全类似,我们直接给出算法.

算法3 (整数减法)

输入:整数 $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)$

$c_k=(a_k-b_k+d_k)/bmod B$
$d_{k+1}=[{a_k-b_k+d_k/over B}]$

辅助序列 $d_k$ 对应于到第 $k$ 位的借位.

同样地,也可以将这个算法改写成适合于编程的更加紧凑的形式.

算法4 (整数减法)

输入:两个 $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$ .

令 $j/leftarrow 0,k/leftarrow 0$ .
令 $w_j/leftarrow (u_j-v_j+k)/bmod B$ , $k/leftarrow -(u_j-v_j+k<0)$ .
$j/leftarrow j+1$ .若 $j<n$ ,转到2;否则返回 $w$ .

如果高精度整数有求补操作,那么仿照机器精度整数补码表示的思想,还可以通过补码加法来做减法:用 $n^*$ 来表示 $m$ 位 $B$ 进制数 $n$ 的补码,则 $n^*=(B^m-1)-n,/quad a-b=(a^*+b)^*.$

高精度整数的加减法算法实际上就是对我们熟悉的十进制笔算加减法的精确化描述并推广到 $B$ 进制记数系统后得到的,用 $T(n)$ 表示用它们计算 $n$ 位 $B$ 进制整数和或差时所需要的 $B$ 进制个位数加减法的次数,可以证明 $T(n)=O(n)$ , 这和必要的输入输出操作相比是都是线性级别的,因此有理由认为这就是理想的算法.尽管如此,在后面有关多项式模算术的讨论中我们将看到,作为多项式模算术的一个直接应用,利用有限域上的整数模算术可以把高精度整数加减法分解成在多个处理器上并行的独立任务,这样一来它就会比经典算法要快.

整数乘法

接下来讨论一位数乘多位数的乘法,这里用到的方法和加减法仍然是相似的,同样可以直接给出算法.

算法5 (一位数乘多位数)

输入:整数 $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)$

$c_k=(a_kb_0+d_k)/bmod B$
$d_{k+1}=[{a_kb_0+d_k/over B}]$

对于多位数乘多位数的情形,如 $b=(b_t/ldots b_1b_0)_B$ ,则 $ab=/sum_{k=0}^tab_kB^k,$ 每一个 $ab_k$ 都是一位数乘多位数,而乘上 $B^k$ 只需在右边添上 $k$ 个0,或者说左移 $k$ 位,于是可以通过 $t$ 次一位数乘多位数, $t$ 次加法和添零(移位)求出 $ab$ 来.除此之外,在后面的"快速乘法"一节中将会看到,对于多位数乘多位数,我们还有更快的算法.

同样地,也可以将这个算法改写成适合于编程的更加紧凑的形式.

算法6 (整数乘法)

输入: $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$ .

令 $w_{m-1},/ldots,w_0/leftarrow 0$ , $j/leftarrow 0$ .
若 $v_j=0$ ,则 $w_{j+m}/leftarrow 0$ ,转到6.
令 $i/leftarrow 0$ , $k/leftarrow 0$ .
令 $t/leftarrow u_i/times v_j+w_{i+j}+k$ , $w_{i+j}/leftarrow t/bmod B$ , $k/leftarrow [/frac{t}{B}]$ .
$i/leftarrow i+1$ .若 $i<m$ ,转到4;否则令 $w_{j+m}/leftarrow k$ .
$j/leftarrow j+1$ .若 $j<n$ ,转到2;否则返回 $w$ .

注1 如果 $w_{m-1},/ldots,w_0$ 不清零,将有 $w=u/times v+(w_{m-1}/cdots w_0)_B$ .

商为一位数的除法

在四则运算的笔算方法中,除法可能是最复杂的,因为除法需要试商,试商包含着猜测的成分,而不能形成有效的算法.机器精度整数除法也要试商,但得益于二进制表示,对它们而言商只有0,1两种情况,只需要比较大小就可以求出.为了有效地计算高精度整数的除法,我们需要将试商的过程算法化[3] ,首先来看商为一位数的除法,即假定商 $/left[/frac{a}{b}/right]$ 满足 $0/le/left[/frac{a}{b}/right]/le B-1$ ,笔算除法的经验告诉我们,被除数和除数的最高几位数对试商是很重要的,我们常常凭借目测(口算)最高几位数就能基本上把商给确定下来.

定理2 设整数 $a=(a_{s+1}a_s/ldots a_1a_0)_B,b=(b_s/ldots b_1b_0)_B$ ,则商 $[{a/over b}]$ 满足不等式 $/left[/frac{a_{s+1}B+a_s}{b_s+1}/right]/le/left[/frac{a}{b}/right]/le/min/left/{/left[/frac{a_{s+1}B+a_s}{b_s}/right],B-1/right/}.$

证明首先 $/frac{a}{b}>/frac{a_{s+1}B^{s+1}+a_sB^s}{(b_s+1)B^s}=/frac{a_{s+1}B+a_s}{b_s+1},$ 故 $/left[/frac{a_{s+1}B+a_s}{b_s+1}/right] /le /left[/frac{a}{b}/right]$ ,不等式的前半部分得证.

而由 $/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.

定理3 若 $b$ 的最高位不小于 $/left[/frac{B}{2}/right]$ ,即 $b_s/ge/left[/frac{B}{2}/right]$ ,则 $q-2/le/left[/frac{a}{b}/right]/le q.$

证明用反证法.假设 $/left[/frac{a}{b}/right]<q-2$ ,则 $/left[/frac{a}{b}/right]/le/left[/frac{a_{s+1}B+a_s}{b_s}/right]-3,/quad/left[/frac{a}{b}/right] /le B-4.$ 因为 $/left[/frac{a}{b}/right]>/frac{a}{b}-1,/quad/frac{a_{s+1}B+a_s}{b_s}</frac{a}{b-B^s},$ 所以 $/frac{a}{b}+2</frac{a}{b-B^s}$ ,推出 $/frac{a}{b}>2(/frac{b}{B^s}-1)/ge 2(b_s-1).$ 因此 $b_s /le /frac{B}{2}-1</left[/frac{B}{2}/right]$ ,与题设条件矛盾,证毕. □

上面只考虑了被除数的头两位与除数的首位,如果允许做更精细的考察,譬如说考虑被除数的头三位与除数的头两位,我们还可以进一步接近商的真值.

定理4 设 $t=(a_{s+1}B+a_s)-qb_s$ ,若 $B/cdot t/ge qb_{s-1}-a_{s-1}$ ,则 $q-1/le/left[/frac{a}{b}/right]/le q.$

证明将 $t$ 代入不等式中得 $B/cdot((a_{s+1}B+a_s)-qb_s)/ge qb_{s-1}-a_{s-1},$ 化简得到 $/begin{align*} a&/ge a_{s+1}B^2+a_sB+a_{s-1}// &/ge q(b_sB+b_{s-1})// &>q(b-B^{s-1}) /end{align*}$ 因为 $q(b-B^{s-1})>qb-B^s/ge(q-1)b$ ,所以 ${a/over b}>q-1$ ,证毕. □

如果定理中的条件不满足,将 $q$ 减一,这时我们总可以说 $q$ 比商的真值至多大1.综合以上这些想法,现在可以写出商为一位数的除法的详细过程了.

算法7 (商为一位数的除法)

输入:整数 $a=(a_{s+1}/ldots a_1a_0)_B$ , $b=(b_s/ldots b_1b_0)_B$ .

输出: $a,b$ 的商 $[{a/over b}]$ .

若 $b_s<[{B/over 2}]$ , $a,b$ 同乘以 $/left[/frac{B}{b_s+1}/right]$ .
计算 $q = /min/left/{/left[/frac{a_{s+1}B+a_s}{b_s}/right],B-1/right/}.$
计算 $t=(a_{s+1}B+a_s)-qb_s$ ,如果 $B/cdot t<qb_{s-1}-a_{s-1}$ , $q$ 减一.
计算余数 $r=a-qb$ ,如果 $r/ge 0$ ,返回 $q$ ,否则返回 $q-1$ .

整数除法

注意到 $r=a-qb=a/bmod b$ 在算法的最后一步也同时被算出,以此为基础可以写出一般的商为多位数的除法算法.

算法8 (整数除法)

输入:整数 $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)$

$q_k=[{B/cdot r_{k+1}+a_k/over b}]$
$r_k=(B/cdot r_{k+1}+a_k)/bmod b$

为了利用快速乘法,还可以利用牛顿迭代来计算整数除法,具体来说就是先利用数值计算中浮点数的牛顿迭代求出除数倒数具有一定精度的近似值,然后利用整数乘法将被除数乘上去,当除数的位数较多时,这种方法是很有效的.

快速乘法

一般原理

高精度整数是计算机代数系统的内置基本类型,四则运算的快慢对系统的性能好坏有着决定性的影响,目前最著名的高精度整数运算库是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乘法

一次多项式的乘法

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算法

综合以上这些想法,现在可以写出Karatsuba乘法的详细过程了[2] .

算法9 (Karatsuba乘法)

输入: $n_a$ 位整数 $a$ 和 $n_b$ 位整数 $b$ .

输出: $a,b$ 的积 $c=a/cdot b$ .

令 $n=/max/{n_a,n_b/}$ ,若 $n=1$ ,利用普通乘法计算并返回 $a/cdot b$ ;若 $n$ 为奇数,令 $n=/frac{n+1}{2}$ ,否则令 $n=/frac{n}{2}$ .
若 $n_a/le n$ ,令 $a_1=0,a_0=a$ ;否则令 $a_1=a$ 的高 $n_a-n$ 位, $a_0=a$ 的低 $n$ 位.
若 $n_b/le n$ ,令 $b_1=0,b_0=b$ ;否则令 $b_1=b$ 的高 $n_b-n$ 位, $b_0=b$ 的低 $n$ 位.
使用Karatsuba乘法计算 $c_0=a_0b_0,c_1=a_1b_1,c_2=(a_0-a_1)(b_0-b_1)$ ,返回 $(c_1/cdot B^{/frac{n}{2}}+c_0+c_1-c_2)/cdot B^{/frac{n}{2}}+c_0.$

与普通乘法类似,在第4步中乘上 $B^k$ 时只需要在右边添上 $k$ 个0,或者说左移 $k$ 位.Karatsuba乘法实现起来并不困难,因此在算法理论中也常常被拿来用作分治算法或者二分法的例子.

Toom-Cook乘法

Toom-3乘法

沿着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})$ .

Lagrange插值公式

更一般地,考虑 $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})$ .

定理5 (Lagrange插值公式) 设 $C(x)=/sum/limits_{k=0}^{n-1}c_kx^k$ , $x_i/neq x_j,0/le i/le j<n$ ,则 $/begin{equation*} /begin{bmatrix} c_0// c_1// /vdots// c_{n-1} /end{bmatrix} = /begin{bmatrix} 1 & x_0 & /cdots & x_0^{n-1}// 1 & x_1 & /cdots & x_1^{n-1}// /vdots & /vdots & /ldots & /vdots// 1 & x_{n-1} & /cdots & x_{n-1}^{n-1} /end{bmatrix} ^{-1} /begin{bmatrix} C(x_0)// C(x_1)// /vdots// C(x_{n-1}) /end{bmatrix}. /end{equation*}$

Lagrange插值公式右端的矩阵是Vandermonde阵,其行列式为 $/prod_{0/le i/le j<n}(x_i-x_j)/neq 0.$ 它告诉我们多项式的系数表示与点值表示之间存在着简单的线性关系.类似地,接下来推广Karatsuba乘法就可以得到Toom-r乘法,简单地说来是这样的,如果 $a,b$ 都是 $rn$ 位 $B$ 进制整数,那么在 $B^n$ 进制下 $a,b$ 都是 $r$ 位数 $a=a_rB^{rn}+/cdots+a_0,/quad b=b_rB^{rn}+/cdots+b_0,$ 其中 $a_i,b_i(i=0,1,/ldots,r)$ 都是 $B^n$ 进制下的个位数,也就是 $B$ 进制下不超过 $n$ 位的数,记 $T(n)$ 是递归地利用Toom-r乘法将两个n位 $B$ 进制数相乘所需要的 $B$ 进制个位数乘法次数,那么有 $T(rn)/le (2r+1)T(n)$ ,由此推出 $T(n)=O(n^{/log_{r+1}{2r+1}})$ .

这样得到的一系列快速乘法算法统称为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乘法

快速傅立叶变换(FFT)[6] 被认为是20世纪数值计算和算法领域最重要的成果之一,它常常被用于快速计算卷积,这对于信号处理等工程领域尤其重要.对于符号计算来讲,FFT实际上是多项式快速多点求值和快速插值算法一个高度优化的特例.不过,在这里我们只需要利用多项式乘法和向量卷积之间的相似性,就可以很好地利用FFT来加速多项式乘法和整数乘法.在上一节中,Toom-Cook乘法通过充分地利用分治与递归有效地减少了多项式多点求值和多项式插值操作的消耗;下面将要讨论的FFT 乘法则是通过精心地挑选求值点,把系数表示与点值表示之间的转换所需的操作压缩为次线性级别,即 $O(n/log{n})$ ,其中 $n$ 代表相乘的两个整数的位数.

原根

顾名思义,FFT是用来快速计算离散傅立叶变换(DFT)的方法,我们可以将DFT简单的理解成以单位复根作为求值点时多项式系数表示到点值表示之间的转换,关于它的原始定义及详细介绍可以在任何一本有关数值计算的教科书中找到.

为了更好地理解DFT,首先来看一下什么是原根,这是一个代数数论的概念.

定义2 (原根) 设 $n/in/mathbb{N}_+$ , $R$ 是一个环, $/omega/in R$ .

若 $/omega^n=1$ ,则称 $/omega$ 为 $R$ 的 $n$ 次单位根.
若 $/omega$ 是 $n$ 次单位根,并且对于 $n$ 的每个素因子 $p$ , $/omega^{n/p}-1$ 不是 $R$ 的零因子,则称 $/omega$ 为 $R$ 的 $n$ 次单位原根.

按照这个定义,任意素数的整数次幂单位复根都是复数域 $/mathbb{C}$ 的单位原根,例如 $/omega=e^{{2/pi i/over 8}}/in R=/mathbb{C}$ 就是一个8次单位原根.除此之外,对于模整数环 $/mathbb{Z}_{17}$ , $/overline{3}/in R=/mathbb{Z}_{17}$ 是16次单位原根,与此同时,虽然 $/overline{2}^{16}=/overline{1}/in R=/mathbb{Z}_{17}$ ,但 $/overline{2}$ 并不是 $/mathbb{Z}_{17}$ 的16次单位原根,因为 $2^8-1/equiv 0/mod 17$ .如果 $q$ 是素数的整数次幂,那么当且仅当 $n|q-1$ 时有限域 $/mathbb{F}_q$ 存在 $n$ 次单位原根.

从定义容易导出原根的一些特殊性质.

定理6 如果 $/omega$ 是环 $R$ 的 $n$ 次单位原根,那么 $/forall 1<k<n$ , $/omega^k-1$ 不是 $R$ 的零因子,并且 $/sum_{j=0}^{n-1}/omega^{jk}=0.$

DFT

另外一个需要用到的概念是循环卷积.

定义3 (循环卷积) 设 $f(x)=/sum/limits_{k=0}^{n-1}f_kx_k,g(x)=/sum/limits_{k=0}^{n-1}g_kx_k$ ,则称 $f(x)*g(x)=/sum/limits_{k=0}^{n-1}h_kx_k$ 为 $f(x),g(x)$ 的循环卷积,其中 $h_k=/sum/limits_{i+j/equiv k/bmod n}f_ig_j$ .

不难发现,循环卷积的定义跟多项式乘积的定义很相似,我们可以简单地把它理解成先扩充多项式次数再做多项式乘法的过程.设 $f,g$ 的多项式乘积 $fg$ 是 $n-1$ 次多项式,则 $f$ 和 $g$ 的次数可能小于 $n-1$ ,如果在它们的最高次项之前添上系数为0的高次项并把它们看成 $n-1$ 次多项式,这时候会发现 $f,g$ 的循环卷积 $f*g$ 恰好等于 $fg$ ,这就是循环卷积和普通多项式乘积之间的关系.

现在可以正式地给出DFT的定义了.

定义4 (离散傅立叶变换) 设 $f(x)=/sum/limits_{k=0}^{n-1}f_kx_k$ , $/omega$ 是环 $R$ 的 $n$ 次单位原根,则称 $DFT_{/omega}(f)=(f(1),f(/omega),/ldots,f(/omega^{n-1})$ 为 $f$ 的离散傅立叶变换.

在下面的定理我们将会发现,原根的良好性质保证了DFT的逆变换和DFT具有相同的结构.

定理7 设 $/omega$ 是环 $R$ 的 $n$ 次单位原根,用 $(DFT_/omega)^{-1}$ 表示 $DFT_/omega$ 的逆变换,那么 $(DFT_/omega)^{-1}=/frac{1}{n}DFT_{/omega^{-1}}.$

证明 $DFT_/omega$ 是 $R^n$ 上的线性变换，它在自然基下的矩阵表示为 $/vec{y}=V_n(/omega)/vec{a}$ ，即 $/begin{align*} /begin{bmatrix} y_0// y_1// y_2// /vdots// y_{n-1} /end{bmatrix} = /begin{bmatrix} 1 & 1 & /cdots & 1// 1 & /omega & /cdots & /omega^{n-1}// 1 & /omega^2 & /cdots & /omega^{2n-2}// /vdots & /vdots & & /vdots// 1 & /omega^{n-1} & /cdots & /omega^{(n-1)^2} /end{bmatrix} /begin{bmatrix} a_0// a_1// a_2// /vdots// a_{n-1} /end{bmatrix}, /end{align*}$ 其中 $V_n(/omega)$ 为 $n$ 阶范德蒙(Vandermonde)方阵. $V_n(/omega)$ 的 $(i, j)$ 位元素是 $/omega^{ij}$ ， $V_n(/omega)/cdot V_n(/omega^{-1})$ 的 $(i, j)$ 位元素为 $/begin{align*} u_{ij}&=/sum_{0/le k<n}/omega^{ik}(/omega^{-1})^{kj}// &=/sum_{0/le k<n}(/omega^{i-j})^k// &=n/delta_{ij}. /end{align*}$ 因此 $V_n(/omega)/cdot V_n(/omega^{-1})=n(/delta_{ij})=nI_n$ ,即 $(V_n(/omega))^{-1}=/frac{1}{n}V_n(/omega^{-1})$ ,亦即 $(DFT_/omega)^{-1}=/frac{1}{n}DFT_{/omega^{-1}}$ . □

这样一来,系数表示和点值表示之间的正逆转换可以用相同的办法去处理,也就是说可以采用统一的方法来计算单位原根处的多点求值和插值,这是适当地选取插值点带来的好处之一,现在 $f$ 和 $g$ 的循环卷积可以表示成 $f*g={1/over n}DFT_{/omega^{-1}}(DFT_/omega(f)/cdot DFT_/omega(g)),$ 其中 $/cdot$ 表示一维向量之间逐点相乘.

FFT算法

如果采用普通的多项式多点求值方法,譬如说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的详细过程了.

算法10 (FFT)

输入: $n-1$ 次多项式的系数表示向量 $f=(f_0,f_1,/ldots,f_{n-1})$ , $n$ 次单位原根 $/omega_n$ ,其中 $n=2^m$ .

输出: $f$ 的离散傅立叶变换 $DFT_{/omega_n}(f)$ .

若 $n$ 等于1,返回 $f$ .
令 $f^{[0]}=(f_0,f_2,/ldots,f_{n-2})$ , $f^{[ 1]}=(f_1,f_3,/ldots,f_{n-1})$ , $/omega=1$ .
计算 $y^{[0]}=DFT_{/omega_n^2}(f^{[0]})$ , $y^{[ 1]}=DFT_{/omega_n^2}(f^{[ 1]})$ .
$k$ 从0到 $n/2-1$ ,顺次计算 $t=/omega y_k^{[ 1]}$ , $y_k=y_k^{[0]}+t$ , $y_{k+n/2}=y_k^{[0]}-t$ , $/omega=/omega/omega_n$ .
返回 $y$ .

设 $T(n)$ 是利用FFT进行多点求值时所需要的系数乘法的次数,那么 $T(n)=O(n/log{n})$ ,如果环R为复数域,那么利用FFT进行整数乘法时需要考虑浮点运算的截断误差和截断位数,这种乘法算法的复杂度为 $O(n/log^*{n})$ ,其中 $/log^*{n}=(/log{n})(/log{/log{n}})/cdots$ .采用FFT计算DFT比直接按照定义计算要快得多,它成功的主要原因是利用了单位原根 $/omega$ 的特殊性质和分治策略.快速傅立叶变换在信号处理,多媒体数据压缩等领域被广泛使用,人们对基本算法进行了许多改进以获得更高的速度和适应于特定的硬件[6] .

有限域上的FFT

如果 $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] Donald E. Knuth, 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] 潘金贵等译，H. Corman, 算法导论, 机械工业出版社, 2006.

[6] Joachim von zur Gathen and Jürgen Gerhard, Modern Computer Algebra, Cambridge University Press, 2002.

转自:http://mathmu.cn/doc/Arithmetic.html#ref-3