ショーンハーゲ・ストラッセン法

ショーンハーゲ・ストラッセン法は、大きな整数の掛け算を高速に計算できる乗算アルゴリズムである。1971年にアーノルド・ショーンハーゲとフォルカー・シュトラッセンによって発見された。2進法で n 桁の整数同士の掛け算の時間計算量はランダウの記号を用いて O ( n ⋅ log ⁡ ( n ) ⋅ log ⁡ ( log ⁡ ( n ) ) ) {\displaystyle O{\Big (}n\cdot \log(n)\cdot \log {\big (}\log(n){\big )}{\Big )}} である。このアルゴリズムでは数論変換の1つである 2n+1個の要素を持つ整数環での高速フーリエ変換を用いる。 ショーンハーゲ・ストラッセン法は、発見された1971年から、フューラーのアルゴリズムが発見される2007年まで十分大きな整数の掛け算の最速アルゴリズムであった。しかし、ショーンハーゲ・ストラッセン法よりフューラーのアルゴリズムの性能が上回るのは天文学的に大きな整数同士の掛け算に限るため、フューラーのアルゴリズムはほとんど実用されていない。 実際には、ショーンハーゲ・ストラッセン法がカラツバ法やToom-3をのような従来の手法より性能が上回り始めるのは、2215 - 2217（10進法で10 000 - 40 000桁）同士の掛け算である。GNU Multi-Precision Libraryはアーキテクチャに依存するが 、64ビット演算において少なくとも1728 - 7808ワード（10進法で33 000 - 150 000桁）の大きさの計算において初めてショーンハーゲ・ストラッセン法を用いる。10進法で74 000桁以上においてショーンハーゲ・ストラッセン法を用いるJavaの実装も存在する。 ショーンハーゲ・ストラッセン法の応用には、GIMPSや円周率の近似計算などの数学の哲学や、整数係数多項式の乗算を整数の掛け算に帰着して効率的に行うクロネッカー置換などの実用的な応用がある。これは楕円曲線法のためにGMP-ECMで用いられている。 ここではショーンハーゲ・ストラッセン法のクランドールとポメランスの Prime Numbers: A Computational Perspective に基づいた実装について主に説明する。ショーンハーゲのオリジナルの実装では Discrete Weighted Transform (DWT) を用いてより効率的に畳み込みを行っている。 クヌースの The Art of Computer Programming でも本アルゴリズムが紹介されている。そこでは、ブロックの構成法について説明し、アルゴリズム全体の手順を順序立てて説明している。 B 進法において、繰り上がりを無視し123 × 456を筆算する場合、以下のようになる。 こうして計算された最下部の数列 (4, 13, 28, 27, 18) を (1,2,3)と(4,5,6) の線形畳み込みもしくは直線畳み込みと呼ぶ。2つの数列の線形畳み込みが計算できれば、元の数の積の計算は単に繰り上がりを実行するだけで容易であり、例えば上の例では右端の8を保持したまま1を27に加算するなどの作業を繰り返すだけである。この例では、繰り上がりを実行することで123 × 456 の積 56088 が得られる。 他にも、便利な畳込みが2種類存在する。入力する数列がn要素である場合の（ここでは n=3）線形畳み込みを考える。右端からn要素を左端からn-1要素に足す。これが巡回畳み込みである。 巡回畳み込みを実行すると、入力の積と Bn − 1 を法として合同な結果が得られる.この例では、103 − 1 = 999であり、(28, 31, 31) に繰り上がりを適用した 3141に対して 3141 ≡ 56088 (mod 999) が成り立つ. 逆に、右端から n 要素から左端の n−1 要素を引く逆向きの巡回畳み込みを行えば、 となる。この結果は Bn + 1 つまり 103 + 1 = 1001 を法として入力の積と合同である。 (28, 23, 5) に繰り上がりを適用した 3035 に対して、 3035 ≡ 56088 (mod 1001) が成り立つ。逆向きの巡回畳み込みは負の数を出力しうるが、繰り上がり・繰り下がりを適切に行うことで除去される。 ショーンハーゲ・ストラッセン法も、他の高速フーリエ変換を用いる乗算と同じように、畳み込み定理の巡回畳み込みを効率的に計算できる性質を用いている。具体的には、 2つのベクトルの巡回畳み込みは、それぞれを一度離散フーリエ変換し、その結果の積を逆離散フーリエ変換することで得られる。 数式で表現すると CyclicConvolution(X, Y) = IDFT(DFT(X) · DFT(Y)) 入力を変換した DFT(X) と DFT(Y) の積を計算するためにも高速フーリエ変換を用いて離散フーリエ変換と逆離散フーリエ変換を行い、乗算アルゴリズムを再帰的に呼び出すことで、巡回畳み込みを効率的に計算できる。 このアルゴリズムは、逆向きの巡回畳み込みを用いれば重みの付いた変換である DWT に対応する畳み込み定理も適用でき、より有用なアルゴリズムとなる。ベクトル X と Y の長さが n であり、 aが 位数 2n の原始根であるとする(つまり、a2n = 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 を用いて mod N において変換を実行する。 複素平面内に1の冪根の累乗の全てが存在するように、任意の位数 n が与えられれば b が mod N において原始根となるような整数 N を選ぶことができる。（つまり、bn ≡ 1 (mod N) でありより次数の少ない b の累乗は mod N において 1と合同でない。） このアルゴリズムは、ほとんどの時間を小さな数値の再帰的な乗算に費やす。ナイーブな実装では、以下のように多くの場所で発生する。 高速フーリエ変換において、原始根 b は繰り返し、累乗・自乗・（他の値との）乗算がされる。 重みベクトル A かその逆行列 A−1 を他のベクトルに掛けるとき、重みベクトル A を形成するために原始根 a の累乗が計算される。 変換されたベクトルのドット積の計算 ショーンハーゲ・ストラッセン法の鍵は、 n と 2N + 1が合同となる法 N の選択である。ここで n は変換されるブロックの数 P=2p の倍数である。これにより、バイナリ形式において大きな整数を表す標準的なシステムに対して多くの利点が生じる。 任意の数は法 2N + 1 における合同な数に削減するために、シフト演算と可算演算しか必要としない。 環における全ての1の冪根は、 2k のかたちで表される。その結果、シフト演算によって任意の数を1の冪根との乗算を計算可能であり、それ自身の指数の演算のみによって累乗と平方根が計算可能である。 変換されたベクトルの、要素ごとの再帰的な乗算は、逆向きの巡回畳み込みによって計算できる。この計算は線形畳み込みより高速であり、「自由に」 mod 2N + 1 において結果を減らせる効果を有する。 再帰的な乗算を便利にするために、ショーンハーゲ・ストラッセン法は単に2つの数の積を計算するのではなく、与えられた n に対して2つの数の mod 2N + 1 での答えを計算する枠組みを作る。これは mod 2N + 1 での積が単に2つの数の積となるほど大きな n を選ぶことが可能であるため一般性を失わない。 アルゴリズムの途中で、2の累乗（1の冪根を含む）との乗算・除算は多くの場合少数のシフト演算・加算によって置き換えられる。これは次の性質を利用している (2n)k ≡ (−1)k mod (2n + 1) ここで 2n を基数とする位取り記数法において k 桁の数は  ( d k − 1 , … , d 1 , d 0 ) {\displaystyle \scriptstyle (d_{k-1},\dots ,d_{1},d_{0})} と表され、 ∑ i = 0 k − 1 d i ⋅ ( 2 n ) i {\displaystyle \scriptstyle \sum _{i=0}^{k-1}d_{i}\cdot (2^{n})^{i}} の数を示し、それぞれは  d i {\displaystyle \scriptstyle d_{i}}   0 ≤ d i < 2 n {\displaystyle \scriptstyle 0\leq d_{i}<2^{n}} である。 この表現により、数を mod 2N + 1 において簡単に削減可能である。最下位である右端から n ビットを取り上げ、次の n ビットを引き、さらに次の n ビットを足す…と全てのビットに対して行う。その結果が 0 から 2n の範囲にない場合は、2N + 1 の倍数を足すことで正規化する。例えば、 n = 3 であり法が23+1 = 9 である場合に 656 は以下のように減らされる。 656 = 1 010 010 0002 ≡ 0002 − 0102 + 0102 − 12 = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 23 + 1). さらに、シフト演算の結果を構築することなく、非常に大きなシフトを行える。0 から 2n の範囲のAを考え、それに 2kを乗じたい場合には、k を n で割って k = qn + r ( r < n )を得る。このとき A(2k) = A(2qn + r) = A[(2n)q(2r)] ≡ (−1)q · shl(A, r) (mod 2n + 1). が成り立つ。ここで、shl(A, r)は A を r ビット左シフトしたものである。A は 2n 以下で、 r < n であるため、r ビット左シフトされた A は高々 2n−1 ビットであり、1回のシフト演算と正規化のための減算のみが必要である。 最後に、 2k で割る場合には、2nが原始根であることを用いれば 22n ≡ 1 (mod 2n + 1) である。したがって A/2k = A(2−k) ≡ A(22n − k) = shl(A, (2n − k) (mod 2n + 1))である。 このアルゴリズムはカラツバ法やToom-3と同様に、分割・評価（高速フーリエ変換）・アダマール積・補間（逆高速フーリエ変換）・結合の順に行われる。 入力である x と y 、そして整数 N が与えられれば、次のアルゴリズムは xy mod  2N + 1 を計算する。N が充分大きい場合、単に xy である。 それぞれの入力を 2k 個の部分に分割し、X と Y とする。(e.g. 12345678 → (12, 34, 56, 78)) 再帰的な乗算のために、小さな N を用意する。そのために、2N/2k + k 以上で 2k で割り切れる最小の整数を n とする。 逆向きの巡回畳み込みによって、mod 2N + 1 における X と Y の積を計算する シフト演算を用いて、X と Y に重みベクトル A を乗ずる。 数論変換高速フーリエ変換を用いてXとYの離散フーリエ変換を計算する。ここで、全ての乗算はシフト演算で行われる。 再帰的にこのアルゴリズムを適用して、変換後の X と Y の要素を掛ける（ドット積）。 3.の結果の逆離散フーリエ変換を計算し、ベクトルCを得る。ここでも全ての積はシフト演算で行われる。これは補間に対応する。 シフト演算を用いて、ベクトルCに重みベクトルの逆行列 A−1 を乗ずる。 符号を調整する：幾つかの要素は負になる。Cのj番目の最大のありえる要素 (j + 1)22N/2k を計算し 2N + 1 を超えていればそれで引く。 最後に、mod 2N + 1 において繰り上がりを実行する。 分割数の最適な数は入力する数の桁数 N に対して N {\displaystyle {\sqrt {N}}} であり、O(N log N log log N) である。k は入力サイズとアーキテクチャに基づいて経験的に設定され、典型的には 4 から 16の値が用いられる。 ステップ2において、次の条件を満たしているかを確認する。 入力ベクトルの各要素は n/2k ビット以下か 入力ベクトルの要素の積は 2n/2k 以下か 畳み込みの各要素は 2k 個以下の総和であるため 2n/2k + k ビット未満か。 n は 2k を割り切れるか 実際のシステムでショーンハーゲ・ストラッセン法を実装する際に用いられた実用的な最適化について述べる。主に2007年のゴードリー、Kruppa、ジマーマンのGNU Multi-Precision Libraryの機能拡張によるものである。 特定のカットオフポイント以下ではToom-3などの他のアルゴリズムを使用して再帰的に乗算を実行するほうが効率的である。結果を  mod 2n + 1 を用いて削減する必要があるが、シフト最適化で説明した通り効率的に計算可能である。 逆離散フーリエ変換には各要素を 22n/2k  で割る計算が含まれる。この演算は同様に2の累乗で割る操作を含む A−1 との乗算と組み合わせることが多い。 大きな数を 2w ビットのワードで表すシステムでは、ベクトルのサイズ 2k が1ワードあたりのビット数の倍数とすると便利である。(例えば32ビットのコンピュータでは k ≥ 5 、64ビットのコンピュータでは k ≥ 6 )。これにより、ビットシフト無しで入力を分割可能であり、高位のワードが0か1の場合に mod 2n + 1 で一様な表現ができる。 正規化には 2n + 1 の加算か減算が含まれる。この値は2つのビット以外0であり、平均計算時間は定数時間で行える。 クーリー–テューキー型FFTアルゴリズムのような反復高速フーリエ変換アルゴリズムは、複素数のベクトルを用いて高速フーリエ変換を行うが、ショーンハーゲ・ストラッセン法では参照の局所性はあまり生じない。高速フーリエ変換の直接的で再帰的な非現実的な実装は、全ての演算が再帰呼び出しの深さがある一定以上より深くなると参照の局所性が生じる。ゴードリー・Kruppa、ジマーマンはベイリーの4ステップアルゴリズムと、複数の再帰ステップを組み合わせたより高次の基数変換を組み合わせた手法を用いた。 ショーンハーゲが説明した「ルート2のトリック」では、 23n/4−2n/4 が mod 2n+1 でのルート2になることを用いており、 22n ≡ 1より、1の4n乗根でもある。これによって変換の長さを 2k から 2k + 1まで長くすることが可能となった。 ^ A. Schönhage and V. Strassen, "Schnelle Multiplikation großer Zahlen", Computing 7 (1971), pp. 281–292. ^ Martin Fürer, "Faster integer multiplication", STOC 2007 Proceedings, pp. 57–66. ^ Rodney Van Meter and Kohei M. Itoh, "Fast quantum modular exponentiation", Physical Review A, Vol. 71 (2005). ^ Overview of Magma V2.9 Features, arithmetic section: Discusses practical crossover points between various algorithms. ^ Luis Carlos Coronado García, "Can Schönhage multiplication speed up the RSA encryption or decryption?" ^ “MUL_FFT_THRESHOLD”. GMP developers' corner. 2011年11月3日閲覧。 ^ “An improved BigInteger class which uses efficient algorithms, including Schönhage–Strassen”. Oracle. 2014年1月10日閲覧。 ^ R. Crandall & C. Pomerance. ^ Donald E. Knuth, The Art of Computer Programming, Volume 2: Seminumerical Algorithms (3rd Edition), 1997. ^ Pierrick Gaudry, Alexander Kruppa, and Paul Zimmermann.