Schönhage Strassen Algorithm

SSA こおしいず

テーマ • Schönhage–Strassen algorithm • 巨大整数同士の乗算アルゴリズム • 計算量𝑂 𝑛 log 𝑛 log log 𝑛 • 以下、SSA

アルゴリズムの流れ

0.入力 • N bitの非負整数2つを入力として受け取るものとする。 • 当然、出力は2N bitとなる。 • Sampleとして、 • 0x1252 × 0x2223 = 0x02716536を考える。

1.分割 • 入力をそれぞれ適当な要素数𝑘の数列に分割する • このとき、上位半分は0にする。 今回はk=8とした。 8bitの入力を4要素に振り分ける。 2 5 2 1 0 0 0 0 𝑓 𝑥 = 2𝑥 0 + 5𝑥 1 + 2𝑥 2 + 1𝑥 3 3 2 2 2 0 0 0 0 𝑔 𝑥 = 3𝑥 0 + 2𝑥 1 + 2𝑥 2 + 2𝑥 3

2.フーリエ変換 • 数列全体に対し、フーリエ変換を行う。 • 多項式の値を求めることに相当する。 • 法は適当に大きくしなければならない。 𝑓 𝑥 = 2𝑥 0 + 5𝑥 1 + 2𝑥 2 + 1𝑥 3 𝑔 𝑥 = 3𝑥 0 + 2𝑥 1 + 2𝑥 2 + 2𝑥 3 𝑓(20 ) 𝑓(22) 𝑓(24) 𝑓(26 ) 𝑓(28 ) 𝑓(210 ) 𝑓(212 ) 𝑓(214 ) mod 28 + 1 𝑔(20) 𝑔(22 ) 𝑔(24 ) 𝑔(26) 𝑔(28) 𝑔(210 ) 𝑔(212 ) 𝑔(214 ) mod 28 + 1

2.フーリエ変換 • 数列全体に対し、フーリエ変換を行う。 • 多項式の値を求めることに相当する。 • 法は適当に大きくしなければならない。 𝑓 𝑥 = 2𝑥 0 + 5𝑥 1 + 2𝑥 2 + 1𝑥 3 𝑔 𝑥 = 3𝑥 0 + 2𝑥 1 + 2𝑥 2 + 2𝑥 3 0a 76 40 25 ff cf c1 a0 mod 28 + 1 09 ab 01 7b 01 9c 01 5c mod 28 + 1

3.畳み込み • 数列の各要素の積を取る。 • この例では8bit同士の積だが、一般にはここも大きな整数になるた 0 + 5𝑥 1 + 2𝑥 2 + 1𝑥 3 𝑓 𝑥 = 2𝑥 め、再帰的にSSAをすることになる。 𝑔 𝑥 = 3𝑥 0 + 2𝑥 1 + 2𝑥 2 + 2𝑥 3 𝑓(20 ) 𝑓(22) 𝑓(24) 𝑓(26 ) 𝑓(28 ) 𝑓(210 ) 𝑓(212 ) 𝑓(214 ) mod 28 + 1 𝑔(20) 𝑔(22 ) 𝑔(24 ) 𝑔(26) 𝑔(28) 𝑔(210 ) 𝑔(212 ) 𝑔(214 ) mod 28 + 1

3.畳み込み • 数列の各要素の積を取る。 • この例では8bit同士の積だが、一般にはここも大きな整数になるた 0 + 5𝑥 1 + 2𝑥 2 + 1𝑥 3 𝑓 𝑥 = 2𝑥 め、再帰的にSSAをすることになる。 𝑔 𝑥 = 3𝑥 0 + 2𝑥 1 + 2𝑥 2 + 2𝑥 3 𝑓𝑔(20 ) 𝑓𝑔(22 ) 𝑓𝑔(24 ) 𝑓𝑔(26 ) 𝑓𝑔(28 ) 𝑓𝑔(210 ) 𝑓𝑔(212 ) 𝑓𝑔(214 ) mod 28 + 1 𝑔(20) 𝑔(22 ) 𝑔(24 ) 𝑔(26) 𝑔(28) 𝑔(210 ) 𝑔(212 ) 𝑔(214 ) mod 28 + 1

3.畳み込み • 数列の各要素の積を取る。 • この例では8bit同士の積だが、一般にはここも大きな整数になるた め、再帰的にSSAをすることになる。 5a 84 40 68 ff a7 c1 47 mod 28 + 1 09 ab 01 7b 01 9c 01 5c mod 28 + 1

4.逆フーリエ変換 • 逆フーリエ変換を行う • 多項式補間に相当する 5a 84 40 68 ff a7 c1 47 mod 28 + 1 06 13 14 15 10 06 02 00 mod 28 + 1 ℎ 𝑥 = 06𝑥 0 + 13𝑥 1 + 14𝑥 2 + 15𝑥 3 + 10𝑥 4 + 06𝑥 5 + 02𝑥 7

5.繰り上げ • 分割のときの基数を使って整数に戻す。 ℎ 24 = 02716536 6 3 5 6 1 7 2 0 06 13 14 15 10 06 02 00 mod 28 + 1 ℎ 𝑥 = 06𝑥 0 + 13𝑥 1 + 14𝑥 2 + 15𝑥 3 + 10𝑥 4 + 06𝑥 5 + 02𝑥 7

実装

値の表現 • 最小の単位はもちろん1byte=8bit • mod 28𝑘 + 1をしたいときには、kを覚えておき、k byte単位でみなす ことで表現 • 28𝑘 と0の区別がつかなくなる問題については、これを検知して続行不可能 なエラーを返すことで対処 • 多項式は、前節の例でみたように、添字iに𝑥 𝑖 の係数が入っているよ うな表現

プリミティブ • 剰余つき加算 • 繰り上げを考えつつひたすら加算を行い、全体で繰り上げ（高々1）が発生し た際には、1の2の補数(0xFF…F)を加算、その結果繰り上げが発生することを 確認 • 0xFF…Fを加算したのに繰り上げが発生しないということは、ちょうど28𝑘 が出 現したということ→続行不可能 • 剰余付き減算 • 28𝑘 + 1 − 𝑥 > 0を足すことと同値なので、2の補数的に考えて、2+~xを加 算するようなコードを書く。

プリミティブ • 剰余つきシフト演算⇔ 2𝑚 の乗算 • シフト幅をmとする。 • 22∗8𝑘 mod 28𝑘 = 1であることに注意して、mを2*8kで割った余りにしてよい。 • 剰余を考えないシフト結果が [0] … [k-1] [k] … [2k-1] … [2k] [3k-1] であるとき、剰余を考えた結果は、少し考えると右下の図のようになる。 もちろん、シフト演算の時は、図の縦に揃う位置で [0] … 同時に値がでてくることはないので、 [k] … 実質数回のmemmoveとNOT、繰り上げだけで済む。 ー ＋ [2k] … [k-1] [2k-1] [3k-1]

高速フーリエ変換 • 前に書いたプリミティブを揃えれば、複素浮動小数点数FFTを少し修 正するだけ。 • 法28𝑘 + 1においては、 2が位数2𝑘をもつ。 • (2𝑙 )−1 = 22∗8𝑘−𝑙 • 原始𝑙乗根は22∗8𝑘/𝑙 • これらがでてくるのも掛け算ばかりなので、指数部を覚えてシフト演 算をするだけでよい。

パラメタ • 今回の実装では、決定しなければならない重要なパラメタが2つある。 1. 分割数…多倍長整数があるとき、何要素からなる数列に分割するか(例：8) • ほとんど制約はなく、自由に決めるパラメタ。 • 高速フーリエ変換をする要素数に一致し、その法にも関係する • アルゴリズム全体の速度をかなり左右する • “入力サイズとアーキテクチャに基づいて経験的に設定され、典型的には 4 から 16の 値が用いられる。”（Wikipedia） • 今回は全体の実装が終了後、各入力幅で様々な値を試して一番良いものを表に保存。 2. 基数…分割した各要素には、何bit∝何byte幅をもたせるか(例:8bit) • とうぜん少ないほうがいいので、分割数を定めるとだいたい機械的に決まる • FFT、畳込み、IFFTをした時点で、この幅をはみ出すとアルゴリズムが破綻するので、そ うならないように定める。

SSA • 入力幅が小さければKaratsuba法 • そうでなければ、パラメタを定め、それに基づいて「アルゴリズムの 流れ」のとおりに処理を進める 1. 2. 3. 4. 5. 分割 FFT 畳み込み（再帰） IFFT 繰り上げ

結果 • 比較対象およびSSAの再帰のベースとしてKaratsuba法を別途実装し た • Karatsuba法の再帰のベース・verifierとして筆算法も実装した • テストは、乱数を入力とし、小さい入力で筆算法とKaratsuba法の一 致を確認、それ以上はKaratsuba法とSSAでの一致を確認した。

結果 速度比較 SSA vs Karatsuba vs Classical vs GMP

2526sec Classical 2650sec Karatsuba 182sec SSA 268435456 bit x 268435456 bit = 161614249桁くらい

My SSA 182sec 8.73sec GMP 6.1.2-1 mpz_mul 268435456 bit x 268435456 bit = 161614249桁くらい

結果 optimal k 参考文書 今回の実験結果 Optimal factor k 16 14 12 10 8 6 4 2 0 11 Theo Kortekaas: Multiplying large numbers and the Schönhage-Strassen Algorithm 12 13 14 15 16 17 18 19 factor k 20 21 22 23 24 25

参考文書 • Theo Kortekaas: Multiplying large numbers and the SchönhageStrassen Algorithm • Wikipedia英語版 Schönhage–Strassen_algorithm • Wikipedia日本語版 ショーンハーゲ・ストラッセン法