頌哈吉-施特拉森演算法

頌哈吉-施特拉森演算法(英語: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.