查看“︁頌哈吉-施特拉森演算法”︁的源代码

[[File:Integer multiplication by FFT.svg|thumb|350px|頌哈吉-施特拉森演算法演算法是以整數乘法的FFT方法為基礎。圖中所列的是用單純FFT法計算1234乘5678 = 7006652的過程。使用了模數337的[[數論轉換]]，選擇85為1的8次方根。為了說明方便，其中數字使用十進制表示，不用二進制表示。頌哈吉-施特拉森以此方式為基礎，再加上negacyclic convolutions]]
[[File:StrassenSchonhage1979 MFO8655.jpg|thumb|頌哈吉（右）和施特拉森（左）在[[上沃尔法赫数学研究所]]下西洋棋，1979年]]
'''頌哈吉-施特拉森演算法'''（{{lang-en|'''Schönhage–Strassen algorithm'''}}）是漸近快速的大[[整数]][[乘法算法]]。是由{{le|阿諾德·頌哈吉|Arnold Schönhage}}和[[沃爾克·施特拉森]]在1971年發明<ref name="schönhage">A. Schönhage and V. Strassen, "[https://link.springer.com/article/10.1007/BF02242355 Schnelle Multiplikation großer Zahlen] {{Wayback|url=https://link.springer.com/article/10.1007/BF02242355 |date=20210602215804 }}", ''Computing'' '''7''' (1971), pp. 281–292.</ref>。若針對二個n位元的整數，其運行的{{le|位元複雜度|bit complexity}}，若以[[大O符号]]表示，是<math>O(n \cdot \log n \cdot \log \log n)</math>。演算法使用在有2<sup>''n''</sup>+1個元素的[[环 (代数)|环]]上的迭代[[快速傅里叶变换]]，這是一種特別的[[數論轉換]]。

頌哈吉-施特拉森演算法是1971年至2007年之間，漸近最快的乘法演算法，2007年時有一個新的乘法演算法{{le|Fürer演算法|Fürer's algorithm}}，其漸近複雜度較低<ref>Martin Fürer, "[http://www.cse.psu.edu/~furer/Papers/mult.pdf Faster integer multiplication] {{webarchive|url=https://web.archive.org/web/20130425232048/http://www.cse.psu.edu/~furer/Papers/mult.pdf |date=2013-04-25 }}", STOC 2007 Proceedings, pp. 57–66.</ref>，不過Fürer演算法只有在大到天文數字的程度時，其速度才會比頌哈吉-施特拉森演算法快，實務上不會使用Fürer演算法（參考[[银河式算法]]）。

在實務上，當數字大於2<sup>2<sup>15</sup></sup>至2<sup>2<sup>17</sup></sup>（十進制下的10,000到40,000位數）時，頌哈吉-施特拉森演算法的速度會比較早期的[[卡拉楚巴算法]]和[[图姆-库克算法]]要快<ref>{{cite journal |first1=Rodney |last1=Van Meter |first2=Kohei M. |last2=Itoh |title=Fast Quantum Modular Exponentiation |journal=Physical Review |series=A |volume=71 |year=2005 |issue=5 |pages= 052320|doi=10.1103/PhysRevA.71.052320 |arxiv=quant-ph/0408006 |s2cid=14983569 }}</ref><ref>[http://magma.maths.usyd.edu.au/magma/Features/node86.html Overview of Magma V2.9 Features, arithmetic section] {{webarchive|url=https://web.archive.org/web/20060820053803/http://magma.maths.usyd.edu.au/magma/Features/node86.html |date=2006-08-20 }}: Discusses practical crossover points between various algorithms.</ref><ref>Luis Carlos Coronado García, "[http://www.cdc.informatik.tu-darmstadt.de/~coronado/Vortrag/MoraviaCrypt-talk-s.pdf  Can Schönhage multiplication speed up the RSA encryption or decryption?] [https://web.archive.org/web/20110719095830/http://www.cdc.informatik.tu-darmstadt.de/~coronado/Vortrag/MoraviaCrypt-talk-s.pdf Archived]", ''University of Technology, Darmstadt'' (2005)</ref>。[[GNU多重精度运算库]]用這個演算法計算至少1728至7808個64位元字節的數值（十進制下的33,000至150,000位數），依硬體架構而定<ref>{{cite web|title=MUL_FFT_THRESHOLD|url=http://gmplib.org/devel/MUL_FFT_THRESHOLD.html|work=GMP developers' corner|access-date=3 November 2011|url-status=dead|archive-url=https://web.archive.org/web/20101124062232/http://gmplib.org/devel/MUL_FFT_THRESHOLD.html|archive-date=24 November 2010}}</ref>，有一個用Java實現的頌哈吉-施特拉森演算法，可以計算十進位超過74,000位數<ref>{{cite web |title=An improved BigInteger class which uses efficient algorithms, including Schönhage–Strassen |url=https://github.com/tbuktu/bigint/raw/master/src/main/java/java/math/BigInteger.java |publisher=Oracle |access-date=2014-01-10}}</ref>。

頌哈吉-施特拉森演算法的應用包括[[数学哲学]]（例如[[互联网梅森素数大搜索]]以及計算[[圆周率近似值]]），實務的應用包括{{le|克羅內克代入|Kronecker substitution}}，其中將整係數多項式的乘法有效的簡化為大數的乘法，GMP-ECM中用此算法來計算{{le|Lenstra橢圓曲線分解|Lenstra elliptic curve factorization}}<ref name="Gaudry">Pierrick Gaudry, Alexander Kruppa, and Paul Zimmermann. [http://www.loria.fr/~gaudry/publis/issac07.pdf A GMP-based Implementation of Schönhage–Strassen’s Large Integer Multiplication Algorithm][https://web.archive.org/web/20110720104333/http://www.loria.fr/~gaudry/publis/issac07.pdf Archived]. Proceedings of the 2007 International Symposium on Symbolic and Algebraic Computation, pp.167–174.</ref>。

==細節==
此段會說明頌哈吉-施特拉森演算法實現的細節，主要是依Crandall和Pomerance在《Prime Numbers: A Computational Perspective》書中的簡介<ref name="crandall">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}}</ref>。其中說明的和頌哈吉原始的方法有些差異：其中用離散加權轉換（discrete weighted transform）來快速計算{{le|負循環卷積|negacyclic convolution}}。另一個細節資訊的來源是[[高德纳]]的《[[计算机程序设计艺术]]》<ref>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.</ref>。此章節從建立blocks開始，最後會一步一步的說明演算法。

===卷積===

假設要將''B''基底的123和456用直式乘法相乘，''B''夠大，因此計算中不會有進位的情形。結果會是：

{|width=300 style="text-align:right"
|   ||        || 1 || 2 || 3
|-
|   || &times;|| 4 || 5 || 6
|-
|colspan=5|<hr/>
|-
|   ||        || 6  || 12 || 18
|-
|   ||      5 || 10 || 15 ||
|-
| 4 ||      8 || 12 ||    ||
|-
|colspan=5|<hr/>
|-
| 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）：

{|width=300 style="text-align:right"
|   || 28 || 27 || 18
|-
| + ||    ||  4 || 13
|-
|colspan=5|<hr/>
|-
|   || 28 || 31 || 31
|}

若在循環卷積的結果再考慮其進位，所得結果等效於原來兩數字的乘積mod B<sup>''n''</sup>&nbsp;−&nbsp;1。在此例中，10<sup>3</sup>&nbsp;−&nbsp;1 = 999，將(28, 31, 31)計算其進位後得到 3141，而3141 ≡ 56088 (mod 999)。

相對的，若將最右邊的''n''個元素減去最左邊的''n''−1個元素，會產生{{le|負循環卷積|negacyclic convolution}}（negacyclic convolution）：

{|width=300 style="text-align:right"
|         || 28 || 27 || 18
|-
| − ||    ||  4 || 13
|-
|colspan=5|<hr/>
|-
|         || 28 || 23 ||  5
|}

若在負循環卷積的結果再考慮其進位，所得結果等效於原來兩數字的乘積mod B<sup>''n''</sup>&nbsp;+&nbsp;1。此例中，10<sup>3</sup>&nbsp;+&nbsp;1 = 1001。將(28, 23, 5)計算進位後得到3035，而3035 ≡ 56088 (mod 1001)。負循環卷積可能會有負數，可以用借位的方式消除其負數。
===卷積定理===
頌哈吉-施特拉森演算法和其他以快速傅立葉變換為基礎的[[乘法算法]]一樣，在本質上都是以[[卷积定理]]為基礎，可以快速的計算二個數列的循環卷積，卷积定理如下：

:兩個向量的循環卷積可以將這兩個向量求其[[离散傅里叶变换]]（DFT），將所得向量依元素相乘，將所得結果再進行反离散傅里叶变换（IDFT）。

若依符號表示，可表示如下：

:CyclicConvolution(''X, Y'') = IDFT(DFT(''X'') &middot; DFT(''Y''))

若用[[快速傅里叶变换]]（FFT）來計算DFT及IDFT，再用遞迴的方式處理乘法演算法，來將DFT(''X'')和DFT(''Y'')相乘，就可以得到計算循環卷積的快速演算法。

在此演算法中，若計算負循環卷積會更有效：使用修改後的卷積定理版本（參考[[數論轉換]]）也可以有相同效果。假設向量X和Y的長度是''n''，而''a''是2''n''[[階 (群論)|階]]的[[单位根]]（意思是''a''<sup>2''n''</sup> = 1，其他''a''的較小幂次都不等於1），可以定義第三個向量''A''，稱為加權向量：

:''A'' = (''a''<sup>''j''</sup>), 0 &le; ''j'' < ''n''
:''A''<sup>−1</sup> = (''a''<sup>−''j''</sup>), 0 &le; ''j'' < ''n''

負循環卷積可以表示為下式：

:NegacyclicConvolution(''X'', ''Y'') = ''A''<sup>−1</sup> &middot; IDFT(DFT(''A'' &middot; ''X'') &middot; DFT(''A'' &middot; ''Y''))

這和循環卷積的公式類似，只是輸入要先乘以''A'',輸出要乘以''A''<sup>−1</sup>。

===代數環的選用===

离散傅里叶变换是抽象的運算，因此可以在任意[[环 (代数)|环]]的代數結構下運作。一般而言是在複數下進行，但實際上為了要有準確的結果，要處理複數運算到足夠的精度，這種乘法不但慢，而且容易有錯誤。因此會用[[數論轉換]]的方式進行，是在整數模''N''（''N''為特定整數）的底下進行。

就像在複數平面上，每一階都可以找到原根一樣，給定任何的階數''n''，可以選擇適當的N，使得''b''是在整數模''N''的系統下，''n''階的原根（''b''<sup>''n''</sup> ≡ 1 (mod ''N'')，''b''其他較小的次幂，都不會是1 mod ''N''）。

此演算法會花很多時間演算小數字的次幂，若用一個天真的方式來處理演算法，這會出現許多次。

# 在FFT演算法中，原根''b''會一直的計算幂次，平方，並且和其他的數相乘。
# 在計算原根''a''的次方，以計算加權向量A時，以及將A或A<sup>−1</sup>和其他向量相乘的時候。
# 在進行向量項對項的相乘時。

頌哈吉-施特拉森演算法的關鍵是選擇模數N等於2<sup>n</sup>&nbsp;+&nbsp;1，其中的''n''是要轉換件數''P''=2<sup>''p''</sup>的倍數。若標準系統，用二進位的形式表示大數值，有許多的好處：

* 所有的數字計算modulo 2<sup>''n''</sup>&nbsp;+&nbsp;1，都可以只用移位和加法快速求得（這會在[[#移位的優化|下一節]]解釋。）
* 此轉換會用到的所有原根都可以寫成2<sup>''k''</sup>，因此原根和其他數字的相乘只需要用到移位，原根的平方和幂次都只是在其指數上的變化。
* 轉換向量的項對項遞迴乘法可以用逆循環卷積來計算，這比卷積要快，而且其結果會自動計算mod 2<sup>''n''</sup>&nbsp;+&nbsp;1。

為了讓遞迴乘法比較方便，會將頌哈吉-施特拉森演算法定位為二個數字的乘積mod 2<sup>''n''</sup>&nbsp;+&nbsp;1（''n''為某整數），而不是一般的乘法。這不會失去演算法的一般性，因為可以選擇夠大的''n''，使得乘積mod 2<sup>''n''</sup>&nbsp;+&nbsp;1就是乘積本身。
===移位的優化===
在此演算法中，有許多和二個次幂相乘或相除（包括所有的原根）的情形是可以用較小數字的移位和相加達成。原因是因為觀察到以下的算式：

:(2<sup>''n''</sup>)<sup>k</sup> &equiv; (−1)<sup>k</sup> mod (2<sup>''n''</sup> + 1)

其中2<sup>''n''</sup>進制下的''k''位數，若用[[进位制|位值計數法]]表示，可以寫成(''d''<sub>''k''&minus;1</sub>,...,''d''<sub>1</sub>,''d''<sub>0</sub>)，代表數值<math>\scriptstyle \sum_{i=0}^{k-1}d_{i}\cdot(2^{n})^{i}</math>，針對每一個''d''<sub>''i''</sub>可得0 ≤ ''d''<sub>''i''</sub> < 2<sup>''n''</sup>。

因此可以簡化表示為二進制數字進行mod {{math|2<sup>''n''</sup>&nbsp;+&nbsp;1}}的計數：取最右邊（最小位）的''n''位元，減去較高位的''n''位元，再加上再高位的''n''位元，一直計算到所有位元都已處理為止。若所得的數超過0到2<sup>''n''</sup>的範圍，可以加或減模數{{math|2<sup>''n''</sup>&nbsp;+&nbsp;1}}的倍數以進行正規化。例如，若''n''=3（因此模數是2<sup>3</sup>+1 = 9），要化簡的數字是656，可以得到：

:656 = 1010010000<sub>2</sub> &equiv; 000<sub>2</sub> − 010<sub>2</sub> + 010<sub>2</sub> − 1<sub>2</sub> = 0 − 2 + 2 − 1 = −1 &equiv; 8 (mod 2<sup>3</sup> + 1).

甚至，針對很多位元的移位，還有簡化的方式。假設數字A介於0到2<sup>''n''</sup>之間，希望乘以2<sup>''k''</sup>。將''k''除以''n''可得到''k'' = ''qn'' + ''r''，其中''r'' < ''n''。因此可得：

:A(2<sup>''k''</sup>) = A(2<sup>''qn'' + ''r''</sup>) = A[(2<sup>''n''</sup>)<sup>''q''</sup>(2<sup>''r''</sup>)] &equiv; (−1)<sup>''q''</sup> {{le|左位移|Shift_Left}}(A, ''r'') (mod 2<sup>''n''</sup> + 1).

因為A≤ 2<sup>''n''</sup>，且''r'' < ''n''。左移''r''位元後最多有2''n''−1位元，因此只需要一次位移和減法（之後可能要正規化）。

最終，若要除以2<sup>''k''</sup>，可將此段的第一個式子平方，可得：

:2<sup>2''n''</sup> &equiv; 1 (mod 2<sup>''n''</sup> + 1)

因此

:A/2<sup>''k''</sup> = A(2<sup>−''k''</sup>) &equiv; A(2<sup>2''n'' − ''k''</sup>) = Shift_Left(A, 2''n'' − ''k'') (mod 2<sup>''n''</sup> + 1).

===演算法===
此演算法有分割、評估（前向FFT）、各項相乘、內插（反FFT）以及合併的階段，類似[[Karatsuba算法]]和[[Toom-Cook算法]]。

假設輸入數字是''x''和''y''，以及整數''N''，以下演算法可以計算''xy'' mod 2<sup>''N''</sup>&nbsp;+&nbsp;1的結果。假設N夠大，使得乘積取模數後的結果仍是乘積本身。

# 將輸入的數字分割為向量X和Y，各有2<sup>''k''</sup>個元素，其中2<sup>''k''</sup>是''N''的因素（例如將12345678 &rarr; (12, 34, 56, 78)）。
# 為了運算方便，需要用較小的''N''以進行遞迴乘法。因此選擇''n''是大於2''N''/2<sup>''k''</sup> + ''k''，而且可被2<sup>''k''</sup>整除的最小數字。
# 利用負循環卷積計算XY mod 2<sup>''n''</sup>&nbsp;+&nbsp;1：
## 將X和Y和加權向量A相乘，作法是使用移位（將第''j''項左移''jn''/2<sup>''k''</sup>）
## 用數論FFT計算X和Y的DFT（將所有乘法用移位表示：針對第2<sup>''k''</sup>次原根，使用2<sup>2''n''/2<sup>''k''</sup></sup>）。
## 遞迴的應用此演算法，將轉換後的X和Y對應項相乘。
## 計算所得向量的IDFT，得到結果向量C（用移位代替乘法），這對應內插階段。
## 將結果向量C和A<sup>−1</sup>相乘，用移位來進行。.
## 調整符號：結果的一些項可能是負的。可以計算C的第''j''項的最大可能正值{{math|(''j'' + 1)2<sup>2''N''/2<sup>''k''</sup></sup>}}，若有超過，再減去模數2<sup>''n''</sup>&nbsp;+&nbsp;1。
# 最後，處理進位mod 2<sup>''N''</sup>&nbsp;+&nbsp;1，得到最終的結果。

輸入數字最佳的分割組數和<math>\sqrt{N}</math>成比例，其中''N''是輸入數字的位元數，此設法會讓執行時間的複雜度為O(''N'' log ''N'' log log ''N''),<ref name="schönhage"/><ref name="crandall"/>，因此需依規則設定參數''k''。實務上，會依輸入數字大小以及演算法架構而定，一般會介於4到16之間<ref name="Gaudry"/>。

在步驟2，已觀察到以下的現象：
* 輸入向量的每一項最多只有''n''/2<sup>''k''</sup>個位元。
* 二個輸入向量項的乘積最多有2''n''/2<sup>''k''</sup>個位元。
* 卷積的每一個元素是最多2<sup>''k''</sup>個乘積的和，因此不會超過 2''n''/2<sup>''k''</sup> + ''k''個位元。
* ''n''一定要可被2<sup>''k''</sup>整除，以確保在步驟1的遞迴呼叫中2<sup>''k''</sup> 整除''N''"的條件會成立。

==參考資料==
{{reflist}}

{{数论算法}}

{{DEFAULTSORT:Schonhage-Strassen Algorithm}}
[[Category:计算机算术算法]]
[[Category:乘法]]