基于kary消减的快速最大公约数算法

文章时发现，所有的文章都是Jebelean作为主

要作者，而且国外文章有共同一作的说法，所以，这里按照固有习惯如此安排。改为文献1

中指出，在典型的代数计算中有一多半的时间都花在大数的GCD计算上。在密码学中，GCD算法主要出现在大数分解和素性判别中，GCD算法的扩展算法普遍被用于求逆。

GCD算法从提出至今已经被广泛研究。在GCD算法中比较著名的算法是辗转相除法，也叫作Euclidean算法，该经典算法被D.Knuth称为所有算法的祖先[2]。二进制GCD（Binary GCD）算法主要适合于二进制算数，由Stein[3]于1961年提出，对于小整数的GCD计算有很好的实现效率。PM GCD（PlusMinus GCD）补充中英全称没有对应中文，英文全称为PlusMinus GCD。

算法[4]是Binary GCD算法的一个推广，它在硬件实现上有很大的优势。对于输入规模为n比特的两个整数，这些算法平均复杂度都是O（n2）。

GCD计算的第一个准线性时间算法由D. Knuth于1970年提出，这个GCD算法基于快速傅里叶变换（Fast Fourier Transformation，FFT）乘法，能在O（n log5 n log log n）的时间内计算两个n位整数的GCD，该算法的复杂度被Schnhage改进为O（n log2 n log log n）[5]。虽然KnuthSchnhage算法被认为是最快的GCD算法，但是在软件实现时过于复杂不实用且由于基于FFT乘法导致只能对大规模输入加速。2004年，Stehlé等[6]提出一种迭代GCD算法，该算法在复杂度方面上没有提高，但实现了对KnuthSchnhage算法的简化。

近二十年来，许多针对GCD算法的改进被提出。1994年，Sorenson[7]提出了对Binary GCD算法思想的推广的右移kary 消减方法，它有许多不同的版本，包括Sorenson[7]的左移、右移，以及Jebelean[8]和Weber[9]的右移kary GCD算法。

基于右移kary消减思想提出的算法中，Weber[9]的算法加速效果最好，可以对Biay GCD算法提速4.5倍，然而这种提速只是针对4096比特规模的输入。在考虑混合Binary GCD算法和Euclidean算法的算法中，文献[10]提出了结合两种算法优势的算法；但是该算法并不能对算法的复杂度O（n2）进行描述是否有误，是否是降低复杂度优化，而且在解决较大规模的输入上并不能做到实现上的优化，更适合于并行实现。Mohamed[11]提出一种混合算法，该算法的优点是便于实现且当减法操作与除法操作时间花费相比较小时该算法表现最优；但是该算法在速度方面上并不比已有算法有显著提高。前面已有这句话，是否重复，应删除这里前面是介绍算法性能等，这里是说明该算法存在的不足之处，以此来说明本文算法的优势。

对于更大规模的输入，KnuthSchnhage算法复杂度较高；但是在软件实现时不实用且由于基于FFT乘法导致只能对大规模输入加速。同样地，虽然文献[6]中算法实现了对KnuthSchnhage算法的简化，但是在软件实现上依然比较复杂。本文基于右移kary消减思想提出一种实现起来相对简单的快速GCD算法，而且该算法能对20万比特以上的输入实现对Binary GCD算法的加速。

右移kary消减想法是寻找整数x、y使得对于给定的输入B和C满足Bx-Cy≡0 mod k，要求B、C和k都互素，k是正整数。对于这样的x、y，由于Bx-Cy被k整除，如果x、y小于k，那么（Bx-Cy）/k的比特长度比B减少约0.5 lb k比特。本文的思路是在Binary GCD的设计基础上进一步利用右移kary消减，可以假定输入是奇数，如果输入不是奇数则可以通过简单的预处理转换为求两个奇数的最大公约数，所以本文考虑的输入都是奇数，在计算过程中产生的偶数都可以右移转换为奇数继续计算下去，这一点可以由Binary GCD的理论分析保证。那么当k为2的方幂时就满足了B、C和k互素的条件。

Euclid算法是从最高比特位优先的角度减小输入整数的规模，Binary GCD则是从最低比特位考虑。本文将结合Binary GCD和Euclid算法的基本思想，同时考虑输入整数的最低位和最高位，给出一个快速GCD算法，为方便说明，记作fast_gcd算法。fast_gcd算法的主要难点在于寻找辅助因子x和y的算法，记作findxy_gcd算法。findxy_gcd算法结合模算法提出fast_gcd算法。本文的算法不同于Binary GCD算法和PM GCD算法的一次循环内有一次相减或相加操作和一次右移操作，在一次循环内计算Bx-Cy并将其右移以减少更多比特。对于相同规模的输入，fast_gcd算法能大规模减少循环量，并且利用模算法来提高算法速度，这是因为对于两个不同长度的输入，模算法能快速降低输入的规模。

本文所有操作都在二进制上进行。对任意正整数A，记Val（A）为满足A=（（A>>x）<

1 findxy_gcd算法

本章首先提出递归算法（recursive_gcd算法），该算法能寻找至少一组合适的x、y，使得Bx-Cy的低比特部分为0。接着该算法以寻找一组合适的x、y，使得Bx-Cy尽量多的低比特部分为0，这是为了能让Bx-Cy尽量右移，大规模减少B的长度，虽然Bx-Cy越多的低比特部分为0，x、y的长度也随之增加，但是相比x、y长度的增加，Bx-Cy能右移规模增加更多，这样就达到了更快地降低输入规模的目的。

虽然recursive_gcd算法能返回至少一组辅助因子x和y，但是由于其失败情形的出现，导致不能最快返回辅助因子，所以才需要findxy_gcd算法寻找辅助因子x和y算法以最快速度返回最合适的x和y。

1.1 递归算法：recursive_gcd算法

1.2 寻找辅助因子x和y算法：findxy_gcd算法

考虑recursive_gcd算法返回一组辅助因子x和y，另一组均为0，也就是失败的情形出现时，如果该结果再被recursive_gcd算法调用则不能返回更优的辅助因子，那么此时由findxy_gcd算法来判断何时停止递归调用以提高算法效率。那么findxy_gcd算法不仅返回了最优的辅助因子，而且省掉了只有recursive_gcd算法时花费在失败递归上面的时间。

在recursive_gcd算法基础上，考虑寻找一组合适的x、y使得Bx-Cy尽量多的低比特部分为0。记（x，y）= findxy_gcd（B， C， L（B）），输入为B、C和 L（B），B > C，输出为（x，y），且满足Val（Bx-Cy）≥2·max{L（x），L（y）}。类似于recursive_gcd算法可以容易得到findxy_gcd算法。

由表2可知Qmb（B，C，x，y）∈[10，20）的出现次数为1658，也就是在findxy_gcd算法的哪个是第4步？为何区间a与第4步相关。为何不用第4）行第4步是：“if x2b==0，

返回（x1b，y1b）；”。之所以区间a和第4步相关，是因为第4步直接给出返回值时

就是区间a的[10，20）情形。第4步表示findxy_gcd算法的4）。

第4）步返回值的概率约为0.1658。另由表2得到比特长度消减量的平均值约为40怎么得出的由比特长度消减量的定义可以直观地看到后面的描述就是比特长度消减量

的定义。

，也就是计算一次Bx-Cy和一次右移可以消减输入的40个比特。虽然这里一次消减40比特相比于输入规模并不大，但是，由下文分析可以看到Binary GCD算法一个循环内只能消减不到3个比特，由此可以看出该算法相比Binary GCD算法在输入规模较大时降低循环次数还是非常明显的，并以此提高Binary GCD的速度。怎么看？文中该句后面紧接着就对这句话给出了例子及解释。

由表2还可以看出来findxy_gcd算法能返回的一组非零值大小与其出现的频率。比特长度消减量以2.1%的概率属于[100，206]，那就说明该算法返回的x、y比特长度大于100的概率约为2.1%，这是因为一般只有当x、y比特长度大于100比特时，才能使得Val（Bx-Cy）≈2L（x）≥200，从而使得Qmb（B，C，x，y）≥100。

2 fast_gcd算法

下面就可以利用findxy_gcd算法得到的x、y来进一步结合Euclid算法思想来消减输入的比特长度。对于求偶数的最大公约数，可以知道容易转化为求奇数的最大公约数，这里就不叙述。fast_gcd算法，即fast_gcd算法描述如下：

对于输入为op1、op2，fast_gcd算法通过求gcd（op1，op2，rop）来得到正确的最大公约数。由于rop长度相对较小，那么就可以两次利用Euclid算法通过求gcd（op1，rop）和gcd（op2， gcd（op1，rop））来得到gcd（op1，op2，rop） = gcd（op1，op2），这是由于rop能被gcd（op1，op2）整除，所以gcd（op1，op2，rop） = gcd（op1，op2）。那么fast_gcd算法，就能够快速且准确地得到最大公约数。

一般情况下findxy_gcd算法返回值的大小在100比特内，则相对于n比特规模输入而言，findxy_gcd算法花费时间可以看作常数cf。 fast_gcd算法花费的时间主要在该算法中的第6）行、第7）行和第11）行，其中第6）行有4个操作，包括2个乘法操作，一般情况为一个100比特的数乘以n比特的数，其复杂度为O（n），1个减法操作和1个右移操作，复杂度为O（n），可见第6）步复杂度为O（n）；第7）步采用的模算法虽然复杂度为O（n2），但是这里采用模算法是为了加速对于Binary GCD算法中1次相减、1次右移的操作，这是因为对于一定规模的输入，一次模算法约消减输入的40比特，而Binary GCD算法要消减40比特，则需要约13次相减和13次右移操作，实际实现时不如模算法花费时间少，而当输入规模达到一定程度导致模算法时间代价更高时则换用一次相减一次右移的方法。则第7）行的复杂度也可以看作O（n）。那么由循环次数约为n/40，再结合第11）行中算法的复杂度为O（n2）得到fast_gcd算法的复杂度为O（n2）。

假如每次调用findxy_gcd算法，输入输出表达为（x，y）= findxy_gcd（B，C，L（B）），都满足L（B）≈2L（x）。那么对于n比特规模的两个输入，findxy_gcd算法所需时间，记为f（n），对于n比特规模的输入，可以得到SchnhageStrassen乘法复杂度[13]记为M（n） = O（n log n log log n）。

4 fast_gcd算法实现结果

在Intel Core i53470，主频为3.2GHz的CPU的PC上，利用GNU MP 4.1.2版本大数包和VC 6.0分别实现Binary GCD和fast_gcd算法得出实验结果见图1和图2。

图1里的Binary GCD算法的一次循环是指进行一次减法操作和一次右移操作，而fast_gcd算法的1次循环是指进行1次findxy_gcd算法的调用、2次乘法操作、2次减法操作和1次右移操作。

由图1看出fast_gcd算法进入循环的次数比Binary GCD算法要少，这也是fast_gcd算法比Binary GCD算法快的最主要原因。由表2知道fast_gcd算法一次循环调用一次findxy_gcd算法至少消减输入10比特，平均消减的比特数约为40，再加上模算法消减的约40比特，那么一次循环平均能消减约80比特。而Binary GCD算法一次循环消减输入的比特数约为3比特，这一点可以从图1看出来，而且文献[14]指出对lb N比特的随机输入的循环次数的期望是k lb N（N→∞， k ≈ 0.706）。那么可以看出fast_gcd算法明显少于Binary GCD算法进入循环的次数。但是由于一个循环内fast_gcd算法花费时间更多，所以对于小规模的输入fast_gcd算法并没有很大程度的提高。

由图2可以看出，当输入规模为20万比特时，两个算法所需时间一样，随着比特规模的增加，fast_gcd算法的速度比Binary GCD算法越来越快。当输入规模达到100万比特规模时，fast_gcd算法的速度是Binary GCD算法的两倍。

5 结语

本文主要提出了寻找辅助因子x和y算法findxy_gcd，并利用findxy_gcd算法对Binary GCD算法进行加速。由于fast_gcd算法的循环次数相比Binary GCD算法循环次数有大幅度的降低，所以使得fast_gcd算法速度更快。在以后的工作中考虑提高recursive_gcd算法返回两组非零值的概率以提高Qmb（B， C， x， y）得到更大数值的概率，或者考虑解决对于一组返回值如何继续进行递归返回两组非零值的问题，这样在面对更大规模的输入时就可以更大幅度地提高fast_gcd算法速度。

参考文献：

文献已查

[1]BUCHBERGER B， JEBELEAN T. Parallel rational arithmetic for computer algebra systems： motivating experiments [M]. Vienna： ACPC-Austrian Center for Parallel Computation， 1993： 1-3.查不到该文献，作者确认是否有误。书名[M].出版地：出版社，出版年：引用页码。

[2]LIOYD E K. The art of computer programming， vol. 2， seminumerical algorithms[J]. Software： Practice and Experience， 1982， 12（9）： 883-884.

[3]STEIN J. Computational problems associated with Racah algebra [J]. Journal of Computational Physics， 1967， 1（3）： 397-405.

[4]BRENT R P， KUNG H T. Systolic VLSI arrays for polynomial GCD computation [J]. IEEE Transactions on Computers， 1984， C33 （8）： 731-736.

[5]YAP C K. Fundamental problems in algorithmic algebra [M]. Oxford： Oxford University Press， 核实是2000版还是1999？2000： 43-76.

[6]STEHL D， ZIMMERMANN P. A binary recursive gcd algorithm [C]// Proceedings of the 2004 6th International Symposium on Algorithmic Number Theory

， LNCS 3076. Berlin： Springer， 2004： 411-425.

[7]SORENSON J. Two fast GCD algorithms [J]. Journal of Algorithms， 1994， 16（1）： 110-144.

[8]JEBELEAN T. A generalization of the binary GCD algorithm [C]// Proceedings of the 1993 International Symposium on Symbolic and Algebraic Computation. New York： ACM， 1993： 111-116.

[9]WEBER K.The accelerated integer GCD algorithm [J]. ACM Transactions on Mathematical Software， 1995， 21（1）： 111-122.

[10]SEDJELMACI S M. The mixed binary Euclid algorithm [J]. Electronic Notes in Discrete Mathematics， 2009， 35查过，无期： 169-176.

[11]MOHAMED F K. A novel fast hybrid GCD computation algorithm [J]. International Journal of Computing Science and Mathematics， 2014， 5（1）： 37-47.

[12]COHN H. Advanced number theory [M]. New York： Courier Dover Publications， 网上是19802012： 55-56.

[13]SCHNHAGE A， STRASSEN V. Schnelle multiplikation grosser zahlen [J]. Computing， 1971， 查不到卷期，核实卷期7（3/4）： 281-292.

[14]BRENT R P. Analysis of the binary Euclidean algorithm [EB/OL]. [2014-12-02]. http：//mathspeople.anu.edu.au/～brent/pd/rpb037a.pdf. Pittsburgh： Department of Computer Science， CarnegieMellon University， 1976： 321-355.

推荐访问： 最大公约数 消减 算法 快速 kary