Next: , Previous: , Up: Multiplication Algorithms   [Index]

15.1.6 FFT乗算

巨大な数の乗算を行うには，FermatスタイルのFFT乗算を使用します。これは Schönhage と Strassen (see References)のアルゴリズムに基づいています。 FFTは色々な形でたくさんのテキストで解説されており，例えばKnuth section 4.3.3 part C や Lipson chapter IX 等があります。ここではGMPで使用されている 形式のものを簡単に解説します。

実行する乗算はNが与えられた時のx*y mod 2^N+1です。 積の完成形は x*yですが，これはN>=bits(x)+bits(y) となるようにNを選び，x と yの大きい方の有効リムにゼロを詰め込むと得ることができます。 この剰余積は，このアルゴリズムの基本なので，このような操作が不可欠になります。

このアルゴリズムは，分割(split), 値評価(evaluate), 点ごとの乗算(pointwise multiply), 補間(interpolate)，KaratsubaアルゴリズムやToom-Cook 3-wayアルゴリズム同様の結合処理(combine) をこの順で実行します。 パラメータkは分割処理を制御し，FFT-k 分割を行って，M=N/2^kビットごとに2^k個に分割します。Nは (2^k)*mp_bits_per_limbなので， この分割でリムの境界内に収めることができ，分割や結合処理におけるビット操作が不要になります。

値評価，点ごとの乗算，補間はすべてモジュロ2^N'+1上で実行されます。 ここでN' は 2M+k+3を2^k及び mp_bits_per_limbの倍数に切り上げられます。 補間した結果， 下記のように入力された分割パーツに対して準巡回畳み込みになります。 N'はこの和が打ち切られない程度に長く取ります。 

           ---
           \         b
w[n] =     /     (-1) * x[i] * y[j]
           ---
       i+j==b*2^k+n
          b=0,1

値評価を行う点はg^i ( i=0 から2^k-1)となります。ここでg=2^(2N'/2^k)です。g は，2^N'+1を法とする1の2^k'乗根で，補間の段階で必要となる 計算の削減ができるようになります。また，これは2のべき乗なので，値の評価と補間に対して実行される Fourier変換はビットシフト，2進加算と符号反転だけで行えます。

この点ごとに行われる乗算は2^N'+1を法として実行され，再帰的にFFTを繰り返し， 最終的には通常の乗算を行って完了します（Toom-3, Karatsuba，筆算アルゴリズムを使用)。 どれを行うかは，サイズN'に対して最適なアルゴリズムを選択します。補間処理は逆FFT になります。結果として出てくるx[i]*y[j]の和の集合は，最終結果にふさわしいオフセット で加算されます。

2乗の計算も同様に実行できますが，入力値xは一つだけになるので，値評価は一つで済み，点毎の乗算は2乗になります。補間についても同様です。

2^N+1を法とする積に対しては，FFT-kアルゴリズムは O(N^(k/(k-1)))となり，この指数は2^k回，各1/2^(k-1)のサイズの入力値に対して乗じた ものになります。 各 kに対しては漸近的に改善が行われますが，オーバーヘッドも，サイズに応じて大きく なっていきます。このコードでは，各kが使われるところで，MUL_FFT_TABLEとSQR_FFT_TABLEを閾値として使用しています。新しいkに対しては効率的に乗算を，シフト，加算，オーバーヘッドに置き換えていきます。

2^N+1を法とする積は通常の，NxN->2N bit 積と減算を行って 求められますので，FFTと Toom-3 等のアルゴリズムは直接比較ができます。k=4 FFT は O(N^1.333)なので，O(N^1.465)となるToom-3よりは高速になることが期待できます。 但し，実際には，MUL_FFT_MODF_THRESHOLD と SQR_FFT_MODF_THRESHOLDは 300 ～ 1000 リムの間で設定されており，CPU ごとに異なります。このサイズを越えた時のみ， 大規模なFFTの再帰が各点における乗算ごとに実行されるようになっています。

FFTを用いて完全積を求めるときにはNを2Nに変更するだけで済み， 理論的な計算量は与えられたkが変わらない限りは同じになります。しかし，FFTが真っ先に使われるような 状況を考えると，FFTは通常の積を求めるために再帰呼び出しされ，その基盤上で，入力サイズ1/2^(k-2)毎に2^k回乗算を実行することになり，結果としてO(N^(k/(k-2)))となります。 これはk=7の場合はO(N^1.4)となり，Toom-3より高速になることが期待できます。 実際には，MUL_FFT_THRESHOLD と SQR_FFT_THRESHOLDはk=8の範囲で，大体3000 ～ 10000リムに設定されています。

Nを 2^k個に分割すると，2M+k+3は2^k と mp_bits_per_limbの倍数に切り上げられますが，これは2^k>=mp\_bits\_per\_limbである時に，Nが効果的になるのは2^(2k-1) bitの倍数の時だからです。この+k+3の意味は，Nが倍数よりも小さい時に，一番近い倍数に切り上げられるということを意味しています。このアルゴリズムの計算量は，望ましいサイスであり，丸めに際してもパディングが発生せず，付加的な+k+3 bitは通常のFFTサイズでは無視できる程度の大きさということを仮定しています。

2^(2k-1)に制限した結果，階段効果(step-effect)が実際の計算速度に影響してきます。 例えばk=8はNを32768 bitの倍数に切り上げられ，32-bitリムの場合は512 リムの グループに分割され，mpn_mul_nは同程度の計算速度になります。k=9の場合は 2048リムのグループに，k=10の場合は8192リムのグループ，といった具合になります。 実際には，各kは小さいサイズの倍数で使用されるので，階段効果の影響が計算時間 vs. サイズの グラフに如実に現れてきます。

閾値の値は現在の実装では各サイズステップの中央に設定されていますが，これはそこそこの チューニングにしかなっていません。というのは，新しいステップの出発地においては， 前のkに戻った方がいいという場合があるからです。もう少し洗練された MUL_FFT_TABLEと SQR_FFT_TABLEの設定方法が知りたいところです。