施特拉森算法

**頌哈吉 - 施特拉森演算法**（英語：**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_ 夠大，因此計算中不會有進位的情形。結果會是：

數列 (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）：

若在循環卷積的結果再考慮其進位，所得結果等效於原來兩數字的乘積 mod Bn− 1。在此例中，103− 1 = 999，將 (28, 31, 31) 計算其進位後得到 3141，而 3141 ≡ 56088 (mod 999)。

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

若在負循環卷積的結果再考慮其進位，所得結果等效於原來兩數字的乘積 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_ 是 2_n_ 階的單位根（意思是 _a_2_n_ = 1，其他 _a_ 的較小冪次都不等於 1），可以定義第三個向量 _A_，稱為加權向量：

　
: _A_ = (_a_j), 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_ 階的原根（_b_n ≡ 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_ 位數，若用位值計數法表示，可以寫成 (_d_k−1,...,_d_1,_d_0)，代表數值 $\scriptstyle \sum _{i=0}^{k-1}d_{i}\cdot (2^{n})^{i}$ ，針對每一個 _d_i 可得 0 ≤ _d_i < 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_ 位元後最多有 2_n_−1 位元，因此只需要一次位移和減法（之後可能要正規化）。

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

　
: 22_n_ ≡ 1 (mod 2n + 1)

因此

　
: A/2k = A(2−_k_) ≡ A(22_n_ − _k_) = Shift\_Left(A, 2_n_ − _k_) (mod 2n + 1).

### 演算法

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

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

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

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

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

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

## 參考資料