Karatsuba Multiplication (GNU MP 6.1.2)

Karatsuba乗算アルゴリズムは，Knuth本の4.3.3節のパートAや他のテキストで解説されていますので，ここでは簡単な紹介に留めます。

入力値がx と yである時，それぞれ，二つに均等分割します（Nが奇数の時は大きい桁の方を1リムだけ短くします）。

 high              low
+----------+----------+
|    x1    |    x0    |
+----------+----------+

+----------+----------+
|    y1    |    y0    |
+----------+----------+

今，bを分割ポイントである2のべき乗とします。即ち，もしx0がk個のリムであるとすれば，(y0も同様)，b=2^(k*mp_bits_per_limb)となります。 x=x1*b+x0，y=y1*b+y0であれば，下記の等式が成立します。

x*y = (b^2+b)*x1*y1 - b*(x1-x0)*(y1-y0) + (b+1)*x0*y0

この等式は，NxNリムの乗算を計算する際，筆算アルゴリズムで実行すると4回の(N/2)x(N/2)乗算が必要になるところ，3回の(N/2)x(N/2)で済んでしまうということを示しています。等式の係数である(b^2+b)，-b，(b+1)は，3つの項がどの位置で加算されるかを表現しているだけで，実際に乗算を行う必要はありません。

 high                              low
+--------+--------+ +--------+--------+
|      x1*y1      | |      x0*y0      |
+--------+--------+ +--------+--------+
          +--------+--------+
      add |      x1*y1      |
          +--------+--------+
          +--------+--------+
      add |      x0*y0      |
          +--------+--------+
          +--------+--------+
      sub | (x1-x0)*(y1-y0) |
          +--------+--------+

(x1-x0)*(y1-y0)の項は絶対値として計算し，符号に応じて加算するか減算するかを選ぶようにします。high(x0*y0)+low(x1*y1)の加算は2回必要になりますが，これは6*kリム長ではなく，5*kリム長の加算を実行することで済み，余計な関数呼び出しによるオーバーヘッド減らすことができます。

2乗の計算も普通の乗法と同様のアルゴリズムで実行しますが，x=yであることを利用すると，2乗計算を3回実行するだけで済むようになります。

x^2 = (b^2+b)*x1^2 - b*(x1-x0)^2 + (b+1)*x0^2

最終演算結果は上記の乗法の場合と同様に，3つの2乗計算を足し込んで得られます。2乗の場合は中間の項(x1-x0)^2は常に正になります。

乗法も2乗計算も中間項を(x1+x0)*(y1+y0)という和の形で構成することも可能ですが，kリムを超える加算が必要になりますし，差の中間項を使った場合に比べて桁上がりとその和の処理が必要になる可能性が出てきてしまいます。

Karatsuba乗算アルゴリズムは漸近的にO(N^1.585)の計算量を必要とします。ここで指数1.585...はlog(3)/log(2)であり，長さ1/2の入力値の乗算を3回実行する，ということを意味しています。筆算アルゴリズムがO(N^2)ですから，かなりの改善がなされたことになり，Karatsuba乗算を実行することで必要となる加算のコストを上回る効果が期待できます。MUL_TOOM22_THRESHOLDは10リム程度に設定しており，SQR（2乗計算)の閾値は必ずこの約2倍になるように設定されています。

筆算アルゴリズムの計算時間はM(N) = a*N^2 + b*N + cという形で与えられ，Karatsuba乗算アルゴリズムではK(N) = 3*M(N/2) + d*N + eという形になり，これを展開するとK(N) = 3/4*a*N^2 + 3/2*b*N + 3*c + d*N + eが得られます。aについている3/4という係数は，筆算アルゴリズムの中の積を減少させた分量を意味しており， M(N)がK(N)よりも小さくなる閾値を押し上げます。逆に，bの係数3/2は比例的に計算量を増加させることを意味しており，K(N)の方がM(N)より小さくなる閾値を押し上げます。後者のケースは，最適化されたmpn_sqr_diagonal関数をmpn_sqr_basecase関数に持ち込んだ時にも当てはまります。どんなスピードアップ法を適用しても計算時間は減少しますので，アルゴリズム切り替えの閾値をどのように決めるかは学術的な課題に過ぎません。

15.1.2 Karatsuba乗算