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

ショーンハーゲ・ストラッセン法は高速フーリエ変換に基づく乗算アルゴリズムである。この図は1234 × 5678 = 7006652を高速フーリエ変換を用いて計算する過程を表している。337を法とする整数環高速フーリエ変換（数論変換）を用いており、85はmod 337 において1の8乗根である。人間が理解しやすいように、2^wを基数とする代わりに、10進法で表す。ショーンハーゲ・ストラッセン法はこのように、畳込みを用いて掛け算を効率的に行う。

ショーンハーゲ・ストラッセン法は、大きな整数の掛け算を高速に計算できる乗算アルゴリズムである。1971年にアーノルド・ショーンハーゲとフォルカー・シュトラッセンによって発見された。2進法で n 桁の整数同士の掛け算の時間計算量はランダウの記号を用いて $O(n\cdot \log (n)\cdot \log (\log (n)\left)\right)$ である。このアルゴリズムでは数論変換の1つである 2ⁿ+1個の要素を持つ整数環での高速フーリエ変換を用いる。

ショーンハーゲ・ストラッセン法は、発見された1971年から、フューラーのアルゴリズムが発見される2007年まで十分大きな整数の掛け算の最速アルゴリズムであった。しかし、ショーンハーゲ・ストラッセン法よりフューラーのアルゴリズムの性能が上回るのは天文学的に大きな整数同士の掛け算に限るため、フューラーのアルゴリズムはほとんど実用されていない。

実際には、ショーンハーゲ・ストラッセン法がカラツバ法やToom-3のような従来の手法より性能が上回り始めるのは、2²¹⁵ - 2²¹⁷（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を筆算する場合、以下のようになる。

		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) の線形畳み込みもしくは直線畳み込みと呼ぶ。2つの数列の線形畳み込みが計算できれば、元の数の積の計算は単に繰り上がりを実行するだけで容易であり、例えば上の例では右端の8を保持したまま1を27に加算するなどの作業を繰り返すだけである。この例では、繰り上がりを実行することで123 × 456 の積 56088 が得られる。

他にも、便利な畳込みが2種類存在する。入力する数列がn要素である場合の（ここでは n=3）線形畳み込みを考える。上の例の（繰り上がり処理をする前の５つの要素について）右端からn要素を左端からn-1要素に足す。これが巡回畳み込み（cyclic convolution）である。

	28	27	18
+		4	13
	28	31	31

巡回畳み込みを実行すると、入力の積と Bⁿ − 1 を法として合同な結果が得られる.この例では、10³ − 1 = 999であり、(28, 31, 31) に繰り上がりを適用した 3141に対して 3141 ≡ 56088 (mod 999) が成り立つ.

同様に（上の例の繰り上がり処理をする前の５つの要素について）逆に、右端から n 要素から左端の n−1 要素を引く逆向きの巡回畳み込み(negacyclic convolution)を行えば、

	28	27	18
−		4	13
	28	23	5

となる。この結果は Bⁿ + 1 つまり 10³ + 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 の原始根であるとする(つまり、a²ⁿ = 1 )。このとき、Aを重みベクトルとして以下のように定義する

A = (a^j), 0 ≤ j < n

A⁻¹ = (a^−j), 0 ≤ j< n

よって、

NegacyclicConvolution(X, Y) = A⁻¹ · IDFT(DFT(A · X) · DFT(A · Y))

といえる。離散フーリエ変換前にAが掛けられ、逆離散フーリエ変換後にA⁻¹が掛けられることを除けばほぼ同じ形である。

環の選択

離散フーリエ変換は、任意の代数環で実行できる抽象的な概念である。一般的には複素数において計算されるが、乗算の正確な結果を保証するため十分に高い精度で計算することは複雑な処理が必要になり、実行速度は遅くなる。ここではその代わりに、数論変換を行い、整数 N を用いて mod N において変換を実行する。

複素平面内にすべての1の冪根とその冪乗が存在するのと同じように、与えられた任意の位数 n に対して、うまく整数 N を選んで　b が法 N における１の原始根となるようにできる （つまり、bⁿ ≡ 1 (mod N) であり、かつ b の冪乗で指数が正でかつ n 未満のものはどれも mod N で 1 と合同にならない。）

このアルゴリズムは、小さな数値の再帰的な乗算にほとんどの時間を費やす。それは素朴な実装では、以下のように多くの場所で発生する。

1. 高速フーリエ変換において、原始根 b は繰り返して、累乗・自乗・（他の値との）乗算がされる。
2. 重みベクトル A か逆の重みベクトル A⁻¹ を他のベクトルに掛けるときには、重みベクトル A を形成するために原始根 a の累乗が計算される。
3. 変換された２つのベクトルの間の成分ごとの積の計算。

ショーンハーゲ・ストラッセン法の鍵となる洞察は法 N の選択であり、ある整数 n を用いてN=テンプレート:Mathであるものを選ぶ。 ここで n は変換を受けるブロックの数 P=2^p の倍数である。このことから大きな整数を二進表現する標準的なシステムにおいては以下のような多くの利点が生じる：

• 任意の整数を法 テンプレート:Math で合同な数に還元することは、シフト演算と加算演算だけでできる。
• 環における全ての1の冪根は、 2^k のかたちになる。そのことから、任意の数の1の冪根による乗算や除算は二進数のシフト演算になり、１の冪根の累乗や平方の計算は冪指数についての計算だけでできる。
• 変換された２つのベクトルの間の要素ごとの積の再帰的な計算は、逆向きの巡回畳み込みによってできる。この計算は線形畳み込みよりも高速であり、 mod (テンプレート:Math )での結果の還元は既にできているので手間が省ける。

再帰的な乗算の実行を便利に行うためには、ショーンハーゲ・ストラッセン法を、2つの整数の積を計算するだけではなくて、与えられた n に対するmod ( テンプレート:Math) での2つの整数の積の計算を行うものとして組み立てる。これによって一般性を失うことはない、なぜならば n を十分大きく選べば、mod(テンプレート:Math )での積が単に2つの整数の積となるようにできるからである。

シフト最適化

アルゴリズムの途中で、2の累乗（1の冪根を含む）との乗算・除算は多くの場合少数のシフト演算・加算によって置き換えられる。これは次の性質を利用している

(2ⁿ)^k ≡ (−1)^k mod (2ⁿ + 1)

ここで 2ⁿ を基数とする位取り記数法において k 桁の数は $(d_{k-1},\dots ,d_1,d_0)$ と表され、${\sum }_{i=0}^{k-1}{d}_i\cdot (2^n{)}^i$ の数を示し、それぞれは $d_i$ $0\le d_i<2^n$ である。

この表現により、数を mod テンプレート:Math において簡単に削減可能である。最下位である右端から n ビットを取り上げ、次の n ビットを引き、さらに次の n ビットを足す…と全てのビットに対して行う。その結果が 0 から 2ⁿ の範囲にない場合は、テンプレート:Math の倍数を足すことで正規化する。例えば、 n = 3 であり法が2³+1 = 9 である場合に 656 は以下のように減らされる。

656 = 1 010 010 000₂ ≡ 000₂ − 010₂ + 010₂ − 1₂ = 0 − 2 + 2 − 1 = −1 ≡ 8 (mod 2³ + 1).

さらに、シフト演算の結果を構築することなく、非常に大きなシフトを行える。0 から 2ⁿ の範囲のAを考え、それに 2^kを乗じたい場合には、k を n で割って テンプレート:Math ( r < n )を得る。このとき

A(2^k) = A(2^{qn + r}) = A[(2ⁿ)^q(2^r)] ≡ (−1)^q · テンプレート:Abbr(A, r) (mod 2ⁿ + 1).

が成り立つ。ここで、shl(A, r)は A を r ビット左シフトしたものである。A は 2ⁿ 以下で、 r < n であるため、r ビット左シフトされた A は高々 2n−1 ビットであり、1回のシフト演算と正規化のための減算のみが必要である。

最後に、 2^k で割る場合には、2nが原始根であることを用いれば

2²ⁿ ≡ 1 (mod 2ⁿ + 1)

である。したがって

A/2^k = A(2^−k) ≡ A(2^{2n − k}) = テンプレート:Abbr(A, (2n − k) (mod 2ⁿ + 1))である。

アルゴリズム

このアルゴリズムはカラツバ法やToom-3と同様に、分割・評価（高速フーリエ変換）・アダマール積・補間（逆高速フーリエ変換）・結合の順に行われる。

入力である x と y 、そして整数 N が与えられれば、次のアルゴリズムは xy mod  テンプレート:Math を計算する。N が充分大きい場合、単に xy である。

1. それぞれの入力を 2^k 個の部分に分割し、X と Y とする。(e.g. 12345678 → (12, 34, 56, 78))
2. 再帰的な乗算のために、小さな N を用意する。そのために、2N/2^k + k 以上で 2^k で割り切れる最小の整数を n とする。
3. 逆向きの巡回畳み込みによって、mod テンプレート:Math における X と Y の積を計算する
   1. シフト演算を用いて、X と Y に重みベクトル A を乗ずる。
   2. 数論変換高速フーリエ変換を用いてXとYの離散フーリエ変換を計算する。ここで、全ての乗算はシフト演算で行われる。
   3. 再帰的にこのアルゴリズムを適用して、変換後の X と Y の要素を掛ける（ドット積）。
   4. 3.の結果の逆離散フーリエ変換を計算し、ベクトルCを得る。ここでも全ての積はシフト演算で行われる。これは補間に対応する。
   5. シフト演算を用いて、ベクトルCに重みベクトルの逆行列 A⁻¹ を乗ずる。
   6. 符号を調整する：幾つかの要素は負になる。Cのj番目の最大のありえる要素 テンプレート:Math を計算し テンプレート:Math を超えていればそれで引く。
4. 最後に、mod テンプレート:Math において繰り上がりを実行する。

分割数の最適な数は入力する数の桁数 N に対して$\sqrt{N}$であり、O(N log N log log N) である。k は入力サイズとアーキテクチャに基づいて経験的に設定され、典型的には 4 から 16の値が用いられる。

ステップ2において、次の条件を満たしているかを確認する。

• 入力ベクトルの各要素は n/2^k ビット以下か
• 入力ベクトルの要素の積は 2n/2^k 以下か
• 畳み込みの各要素は 2^k 個以下の総和であるため 2n/2^k + k ビット未満か。
• n は 2^k を割り切れるか

最適化

実際のシステムでショーンハーゲ・ストラッセン法を実装する際に用いられた実用的な最適化について述べる。主に2007年のゴードリー、Kruppa、ジマーマンのGNU Multi-Precision Libraryの機能拡張によるものである。

特定のカットオフポイント以下ではToom-3などの他のアルゴリズムを使用して再帰的に乗算を実行するほうが効率的である。結果を  mod 2ⁿ + 1 を用いて削減する必要があるが、シフト最適化で説明した通り効率的に計算可能である。

逆離散フーリエ変換には各要素を 2^2n/2k  で割る計算が含まれる。この演算は同様に2の累乗で割る操作を含む A⁻¹ との乗算と組み合わせることが多い。

大きな数を 2^w ビットのワードで表すシステムでは、ベクトルのサイズ 2^k が1ワードあたりのビット数の倍数とすると便利である。(例えば32ビットのコンピュータでは k ≥ 5 、64ビットのコンピュータでは k ≥ 6 )。これにより、ビットシフト無しで入力を分割可能であり、高位のワードが0か1の場合に mod テンプレート:Math で一様な表現ができる。

正規化には テンプレート:Math の加算か減算が含まれる。この値は2つのビット以外0であり、平均計算時間は定数時間で行える。

クーリー–テューキー型FFTアルゴリズムのような反復高速フーリエ変換アルゴリズムは、複素数のベクトルを用いて高速フーリエ変換を行うが、ショーンハーゲ・ストラッセン法では参照の局所性はあまり生じない。高速フーリエ変換の直接的で再帰的な非現実的な実装は、全ての演算が再帰呼び出しの深さがある一定以上より深くなると参照の局所性が生じる。ゴードリー・Kruppa、ジマーマンはベイリーの4ステップアルゴリズムと、複数の再帰ステップを組み合わせたより高次の基数変換を組み合わせた手法を用いた。

ショーンハーゲが説明した「ルート2のトリック」では、 2^3n/4−2^n/4 が mod 2ⁿ+1 でのルート2になることを用いており、 2²ⁿ ≡ 1より、1の4n乗根でもある。これによって変換の長さを 2^k から 2^{k + 1}まで長くすることが可能となった。