颂哈吉-施特拉森算法

颂哈吉-施特拉森算法(英语:Schönhage–Strassen algorithm)是渐近快速的大整数乘法算法。是由阿诺德·颂哈吉英语Arnold Schönhage沃尔克·施特拉森在1971年发明[1]。若针对二个n位元的整数,其运行的位元复杂度英语bit complexity,若以大O符号表示,是。算法使用在有2n+1个元素的上的迭代快速傅里叶变换,这是一种特别的数论转换

颂哈吉-施特拉森算法算法是以整数乘法的FFT方法为基础。图中所列的是用单纯FFT法计算1234乘5678 = 7006652的过程。使用了模数337的数论转换,选择85为1的8次方根。为了说明方便,其中数字使用十进制表示,不用二进制表示。颂哈吉-施特拉森以此方式为基础,再加上negacyclic convolutions
颂哈吉(右)和施特拉森(左)在上沃尔法赫数学研究所下国际象棋,1979年

颂哈吉-施特拉森算法是1971年至2007年之间,渐近最快的乘法算法,2007年时有一个新的乘法算法Fürer算法英语Fürer's algorithm,其渐近复杂度较低[2],不过Fürer算法只有在大到天文数字的程度时,其速度才会比颂哈吉-施特拉森算法快,实务上不会使用Fürer算法(参考银河式算法)。

在实务上,当数字大于2215至2217(十进制下的10,000到40,000位数)时,颂哈吉-施特拉森算法的速度会比较早期的卡拉楚巴算法图姆-库克算法要快[3][4][5]GNU多重精度运算库用这个算法计算至少1728至7808个64位元字节的数值(十进制下的33,000至150,000位数),依硬件架构而定[6],有一个用Java实现的颂哈吉-施特拉森算法,可以计算十进制超过74,000位数[7]

颂哈吉-施特拉森算法的应用包括数学哲学(例如互联网梅森素数大搜索以及计算圆周率近似值),实务的应用包括克罗内克代入英语Kronecker substitution,其中将整系数多项式的乘法有效的简化为大数的乘法,GMP-ECM中用此算法来计算Lenstra椭圆曲线分解英语Lenstra elliptic curve factorization[8]

细节

编辑

此段会说明颂哈吉-施特拉森算法实现的细节,主要是依Crandall和Pomerance在《Prime Numbers: A Computational Perspective》书中的简介[9]。其中说明的和颂哈吉原始的方法有些差异:其中用离散加权转换(discrete weighted transform)来快速计算负循环卷积英语negacyclic convolution。另一个细节资讯的来源是高德纳的《计算机程序设计艺术[10]。此章节从建立blocks开始,最后会一步一步的说明算法。

卷积

编辑

假设要将B基底的123和456用直式乘法相乘,B够大,因此计算中不会有进位的情形。结果会是:

1 2 3
× 4 5 6

6 12 18
5 10 15
4 8 12

4 13 28 27 18

数列(4, 13, 28, 27, 18)称为二个原始数列(1,2,3)和(4,5,6)的线性卷积(linear convolution)或非循环卷积(acyclic convolution)。只要有二个数列的线性卷积,要算出原数字在十进制以下的乘积就很简单了:只需要处理进位(例如,最右栏的数字,只保留8,将1进到上一栏的27)。此例中可以得到正确的乘积56088。

有另外几种也很有用的卷积。假设输入数列有n个元素(假设3个),线性卷积会有n+n−1个元素,若将最右边的n个元素加上最左边的 n−1个元素,可以产生循环卷积(cyclic convolution):

28 27 18
+ 4 13

28 31 31

若在循环卷积的结果再考虑其进位,所得结果等效于原来两数字的乘积mod Bn − 1。在此例中,103 − 1 = 999,将(28, 31, 31)计算其进位后得到 3141,而3141 ≡ 56088 (mod 999)。

相对的,若将最右边的n个元素减去最左边的n−1个元素,会产生负循环卷积英语negacyclic convolution(negacyclic convolution):

28 27 18
4 13

28 23 5

若在负循环卷积的结果再考虑其进位,所得结果等效于原来两数字的乘积mod Bn + 1。此例中,103 + 1 = 1001。将(28, 23, 5)计算进位后得到3035,而3035 ≡ 56088 (mod 1001)。负循环卷积可能会有负数,可以用借位的方式消除其负数。

卷积定理

编辑

颂哈吉-施特拉森算法和其他以快速傅立叶变换为基础的乘法算法一样,在本质上都是以卷积定理为基础,可以快速的计算二个数列的循环卷积,卷积定理如下:

两个向量的循环卷积可以将这两个向量求其离散傅里叶变换(DFT),将所得向量依元素相乘,将所得结果再进行反离散傅里叶变换(IDFT)。

若依符号表示,可表示如下:

CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y))

若用快速傅里叶变换(FFT)来计算DFT及IDFT,再用递回的方式处理乘法算法,来将DFT(X)和DFT(Y)相乘,就可以得到计算循环卷积的快速算法。

在此算法中,若计算负循环卷积会更有效:使用修改后的卷积定理版本(参考数论转换)也可以有相同效果。假设向量X和Y的长度是n,而a是2n单位根(意思是a2n = 1,其他a的较小幂次都不等于1),可以定义第三个向量A,称为加权向量:

A = (aj), 0 ≤ j < n
A−1 = (aj), 0 ≤ j < n

负循环卷积可以表示为下式:

NegacyclicConvolution(X, Y) = A−1 · IDFT(DFT(A · X) · DFT(A · Y))

这和循环卷积的公式类似,只是输入要先乘以A,输出要乘以A−1

代数环的选用

编辑

离散傅里叶变换是抽象的运算,因此可以在任意的代数结构下运作。一般而言是在复数下进行,但实际上为了要有准确的结果,要处理复数运算到足够的精度,这种乘法不但慢,而且容易有错误。因此会用数论转换的方式进行,是在整数模NN为特定整数)的底下进行。

就像在复数平面上,每一阶都可以找到原根一样,给定任何的阶数n,可以选择适当的N,使得b是在整数模N的系统下,n阶的原根(bn ≡ 1 (mod N),b其他较小的次幂,都不会是1 mod N)。

此算法会花很多时间演算小数字的次幂,若用一个天真的方式来处理算法,这会出现许多次。

  1. 在FFT算法中,原根b会一直的计算幂次,平方,并且和其他的数相乘。
  2. 在计算原根a的次方,以计算加权向量A时,以及将A或A−1和其他向量相乘的时候。
  3. 在进行向量项对项的相乘时。

颂哈吉-施特拉森算法的关键是选择模数N等于2n + 1,其中的n是要转换件数P=2p的倍数。若标准系统,用二进制的形式表示大数值,有许多的好处:

  • 所有的数字计算modulo 2n + 1,都可以只用移位和加法快速求得(这会在下一节解释。)
  • 此转换会用到的所有原根都可以写成2k,因此原根和其他数字的相乘只需要用到移位,原根的平方和幂次都只是在其指数上的变化。
  • 转换向量的项对项递回乘法可以用逆循环卷积来计算,这比卷积要快,而且其结果会自动计算mod 2n + 1。

为了让递回乘法比较方便,会将颂哈吉-施特拉森算法定位为二个数字的乘积mod 2n + 1(n为某整数),而不是一般的乘法。这不会失去算法的一般性,因为可以选择够大的n,使得乘积mod 2n + 1就是乘积本身。

移位的优化

编辑

在此算法中,有许多和二个次幂相乘或相除(包括所有的原根)的情形是可以用较小数字的移位和相加达成。原因是因为观察到以下的算式:

(2n)k ≡ (−1)k mod (2n + 1)

其中2n进制下的k位数,若用位值计数法表示,可以写成(dk−1,...,d1,d0),代表数值 ,针对每一个di可得0 ≤ di < 2n

因此可以简化表示为二进制数字进行mod 2n + 1的计数:取最右边(最小位)的n位元,减去较高位的n位元,再加上再高位的n位元,一直计算到所有位元都已处理为止。若所得的数超过0到2n的范围,可以加或减模数2n + 1的倍数以进行正规化。例如,若n=3(因此模数是23+1 = 9),要化简的数字是656,可以得到:

656 = 10100100002 ≡ 0002 − 0102 + 0102 − 12 = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 23 + 1).

甚至,针对很多位元的移位,还有简化的方式。假设数字A介于0到2n之间,希望乘以2k。将k除以n可得到k = qn + r,其中r < n。因此可得:

A(2k) = A(2qn + r) = A[(2n)q(2r)] ≡ (−1)q 左位移英语Shift_Left(A, r) (mod 2n + 1).

因为A≤ 2n,且r < n。左移r位元后最多有2n−1位元,因此只需要一次位移和减法(之后可能要正规化)。

最终,若要除以2k,可将此段的第一个式子平方,可得:

22n ≡ 1 (mod 2n + 1)

因此

A/2k = A(2k) ≡ A(22nk) = Shift_Left(A, 2nk) (mod 2n + 1).

算法

编辑

此算法有分割、评估(前向FFT)、各项相乘、内插(反FFT)以及合并的阶段,类似Karatsuba算法Toom-Cook算法

假设输入数字是xy,以及整数N,以下算法可以计算xy mod 2N + 1的结果。假设N够大,使得乘积取模数后的结果仍是乘积本身。

  1. 将输入的数字分割为向量X和Y,各有2k个元素,其中2kN的因素(例如将12345678 → (12, 34, 56, 78))。
  2. 为了运算方便,需要用较小的N以进行递回乘法。因此选择n是大于2N/2k + k,而且可被2k整除的最小数字。
  3. 利用负循环卷积计算XY mod 2n + 1:
    1. 将X和Y和加权向量A相乘,作法是使用移位(将第j项左移jn/2k
    2. 用数论FFT计算X和Y的DFT(将所有乘法用移位表示:针对第2k次原根,使用22n/2k)。
    3. 递回的应用此算法,将转换后的X和Y对应项相乘。
    4. 计算所得向量的IDFT,得到结果向量C(用移位代替乘法),这对应内插阶段。
    5. 将结果向量C和A−1相乘,用移位来进行。.
    6. 调整符号:结果的一些项可能是负的。可以计算C的第j项的最大可能正值(j + 1)22N/2k,若有超过,再减去模数2n + 1。
  4. 最后,处理进位mod 2N + 1,得到最终的结果。

输入数字最佳的分割组数和 成比例,其中N是输入数字的位元数,此设法会让执行时间的复杂度为O(N log N log log N),[1][9],因此需依规则设定参数k。实务上,会依输入数字大小以及算法架构而定,一般会介于4到16之间[8]

在步骤2,已观察到以下的现象:

  • 输入向量的每一项最多只有n/2k个位元。
  • 二个输入向量项的乘积最多有2n/2k个位元。
  • 卷积的每一个元素是最多2k个乘积的和,因此不会超过 2n/2k + k个位元。
  • n一定要可被2k整除,以确保在步骤1的递回呼叫中2k 整除N"的条件会成立。

参考资料

编辑
  1. ^ 1.0 1.1 A. Schönhage and V. Strassen, "Schnelle Multiplikation großer Zahlen页面存档备份,存于互联网档案馆)", Computing 7 (1971), pp. 281–292.
  2. ^ Martin Fürer, "Faster integer multiplication 互联网档案馆存档,存档日期2013-04-25.", STOC 2007 Proceedings, pp. 57–66.
  3. ^ Van Meter, Rodney; Itoh, Kohei M. Fast Quantum Modular Exponentiation. Physical Review. A. 2005, 71 (5): 052320. S2CID 14983569. arXiv:quant-ph/0408006 . doi:10.1103/PhysRevA.71.052320. 
  4. ^ Overview of Magma V2.9 Features, arithmetic section 互联网档案馆存档,存档日期2006-08-20.: Discusses practical crossover points between various algorithms.
  5. ^ Luis Carlos Coronado García, "Can Schönhage multiplication speed up the RSA encryption or decryption? Archived", University of Technology, Darmstadt (2005)
  6. ^ MUL_FFT_THRESHOLD. GMP developers' corner. [3 November 2011]. (原始内容存档于24 November 2010). 
  7. ^ An improved BigInteger class which uses efficient algorithms, including Schönhage–Strassen. Oracle. [2014-01-10]. 
  8. ^ 8.0 8.1 Pierrick Gaudry, Alexander Kruppa, and Paul Zimmermann. A GMP-based Implementation of Schönhage–Strassen’s Large Integer Multiplication AlgorithmArchived. Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation, pp.167–174.
  9. ^ 9.0 9.1 R. Crandall & C. Pomerance. Prime Numbers – A Computational Perspective. Second Edition, Springer, 2005. Section 9.5.6: Schönhage method, p. 502. ISBN 0-387-94777-9
  10. ^ Donald E. Knuth, The Art of Computer Programming(计算机程序设计艺术), Volume 2: Seminumerical Algorithms (3rd Edition), 1997. Addison-Wesley Professional, ISBN 0-201-89684-2. Section 4.3.3.C: Discrete Fourier transforms, pg.305.