頌哈吉-施特拉森演演算法

<image>

<image>

頌哈吉-施特拉森演演算法（Schönhage–Strassen algorithm）是漸近快速的大整數乘法算法。是由阿諾德·頌哈吉和沃爾克·施特拉森在1971年發明。若針對二個n位元的整數，其運行的位元複雜度，若以大O符號表示，是$O(n \cdot \log n \cdot \log \log n)$。演演算法使用在有2n+1個元素的環上的迭代快速傅立葉變換，這是一種特別的數論轉換。

頌哈吉-施特拉森演演算法是1971年至2007年之間，漸近最快的乘法演演算法，2007年時有一個新的乘法演演算法Fürer演演算法，其漸近複雜度較低，不過Fürer演演算法只有在大到天文數字的程度時，其速度才會比頌哈吉-施特拉森演演算法快，實務上不會使用Fürer演演算法（參考銀河式算法）。

在實務上，當數字大於2215至2217（十進制下的10,000到40,000位數）時，頌哈吉-施特拉森演演算法的速度會比較早期的卡拉楚巴算法和圖姆-庫克演算法要快。GNU多重精度運算庫用這個演演算法計算至少1728至7808個64位元字節的數值（十進制下的33,000至150,000位數），依硬體架構而定，有一個用Java實現的頌哈吉-施特拉森演演算法，可以計算十進位超過74,000位數。

頌哈吉-施特拉森演演算法的應用包括數學哲學（例如網際網路梅森素數大搜索以及計算圓周率近似值），實務的應用包括克羅內克代入，其中將整係數多項式的乘法有效的簡化為大數的乘法，GMP-ECM中用此演算法來計算Lenstra橢圓曲線分解。

## 細節

此段會說明頌哈吉-施特拉森演演算法實現的細節，主要是依Crandall和Pomerance在《Prime Numbers: A Computational Perspective》書中的簡介。其中說明的和頌哈吉原始的方法有些差異：其中用離散加權轉換（discrete weighted transform）來快速計算負循環卷積。另一個細節資訊的來源是高德納的《電腦程式設計藝術》。此章節從建立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）：

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 = (a−j), 0 ≤ j < n

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

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

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

### 代數環的選用

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

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

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

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

頌哈吉-施特拉森演演算法的關鍵是選擇模數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)，代表數值$\scriptstyle \sum_{i=0}^{k-1}d_{i}\cdot(2^{n})^{i}$，針對每一個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 左位移(A, r) (mod 2n + 1).

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

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

22n ≡ 1 (mod 2n + 1)

因此

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

### 演演算法

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

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

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

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

在步驟2，已觀察到以下的現象：

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