3 Karatsuba乗算の一般化

3.1 数と多項式

$b$進$n$桁の数 $u=(u_{n-1}u_{n-2} \ldots u_0)_b$の値は

$$u_{n-1}b^{n-1}+u_{n-2}b^{n-2}+ \cdots +u_0$$

であるが、これは多項式

$$u(x)=u_{n-1}x^{n-1}+u_{n-2}x^{u-2}+ \cdots +u_0$$ (3.28)

において、$x=b$を代入したときの値に他ならない。このようにして数を多項式に対応付けると、数の加減乗除は多項式の加減乗除と考えることができる。例えば加算$w=u+v$には $w(x)=u(x)+v(x)$が対応する。これも$x=b$を代入すれば

$$w(b)=u(b)+v(b)=u+v=w$$

となることからすぐに分かる。例えば$u=128,v=186$の加算$w=u+v=314$を多項式で行うと、

$$\begin{eqnarray*} u(x)&=&x^2+2x+8 \\ v(x)&=&x^2+8x+6 \\ w(x)&=&u(x)+v(x)=2x^2+10x+14 \end{eqnarray*}$$

となり、$w(10)=314$である。しかしここで得られた$w(x)$は$314$をそのまま多項式に置き換えた$3x^2+x+4$とは異なっている。これは多項式の演算には繰り上げが存在しないためである。数の演算の道具として多項式を使う限り、必要な値は$w(b)$だけであるから、上の例($b=10$)の場合は

$$w(x)=3x^2+x+4$$

のように $w(x)$の係数を「繰り上げ」れば良い。

次に乗算を考えてみる。例えば

$$u(x)=u_2x^2+u_1x+u_0,v(x)=v_2x^2+v_1x+v_0$$

のような簡単な例で多項式の乗算$w(x)=u(x)v(x)$を行ってみると、

$$w(x)=u_2v_2x^4+(u_1v_2+u_2v_1)x^3+(u_0v_2+u_1v_1+u_2v_0)x^2+(u_0v_1+u_1v_0)x+u_0v_0$$

のようになる。このことから考えて、2つの多項式

$$u(x)=\sum_{n=0}^{\infty}u_nx^n,v(x)=\sum_{n=0}^{\infty}v_nx^n$$ (3.29)

の積$w(x)=u(x)v(x)$が

$$w(x)=\sum_{n=0}^{\infty}w_nx^n$$ (3.30)

のように表されているとき、$w_n$は

$$w_n=\sum_{k=0}^{n}u_kv_{n-k}$$ (3.31)

で与えられることが分かる。このような形の和を畳み込み(convolution)と呼ぶ。多倍長数の乗算とは要するにこの畳み込みの計算である。筆算ではこの畳み込みをそのまま計算し、Karatsuba乗算では

$$u_0v_1+u_1v_0=u_1v_1+u_0v_0-(u_1-u_0)(v_1-v_0)$$

という関係を使って乗算回数を減らしているのである。

3.2 連立方程式の解としての畳み込み

畳み込みは次のようにしても得られる。多項式

$$\begin{eqnarray*} u(x)&=&u_1x+u_0 \\ v(x)&=&v_1x+v_0 \\ w(x)&=&w_2x^2+w_1x+w_0 \end{eqnarray*}$$

の間に$w(x)=u(x)v(x)$という関係があり、$u(x),v(x)$の係数から$w(x)$の係数を決定したいとすると、未知数が3個であるから、乗算3回を含む連立方程式

$$w(x_i)=u(x_i)v(x_i) \qquad (i=0,1,2)$$ (3.32)

を解けばよいことになる。具体的に $x_0=0,x_1=1,x_2=-1$とすれば

$$\begin{eqnarray*} w_0&=&u_0v_0 \\ w_2+w_1+w_0&=&(u_1+u_0)(v_1+v_0) \\ w_2-w_1+w_0&=&(u_1-u_0)(v_1-v_0) \end{eqnarray*}$$

となって

$$w_0 = u_0v_0$$

$$w_1 = u_1v_1+u_0v_0-(u_1-u_0)(v_1-v_0)$$ (3.33)

$$w_2 = u_1v_1$$

が求まり、Karatsuba乗算のアルゴリズムが導かれる。

この方法は容易に拡張できる。例えば $u(x),v(x),w(x)$として

$$u(x) = u_2x^2+u_1x+u_0$$

$$v(x) = v_2x^2+v_1x+v_0$$ (3.34)

$$w(x) = w_4x^4+w_3x^3+w_2x^2+w_1x+w_0$$

を考え、 $w(x_i)=u(x_i)v(x_i)$ $(i=0,1,\ldots,4)$として乗算5回を含む連立方程式を立てた場合、必要な時間$T(n)$は

$$T(n) = 5T(n/3)+cn/3$$

$$= 5^{\log_{3}n}T(1)+cn((5/3)^{\log_{3}n}-1)$$

$$= n^{\log_{3}5}(T(1)+c)-cn$$

$$= \Theta(n^{\log_{3}5}) \qquad(\log_{3}5 \approx 1.465)$$ (3.35)

となり、Karatsuba乗算よりも速いことが分かる。

この方法を一般に拡張して多倍長数を$p$分割した場合、上と同様の議論から必要な時間は $T(n)=\Theta(n^{\log_{p}(2p-1)})$となる。$p$を大きくすれば次数をかなり下げることができるが、その場合には解を求める式が煩雑になるので、機械的に解を求めるアルゴリズムが必要となる。

3.3 多項式補間

問題は$(x_i,y_i)$ $(i=0,\ldots,n-1)$という点を通る多項式 $w(x)=\sum_{k=0}^{n-1}w_kx^k$を決定することである。これを多項式補間(Polynomial Interpolation)という。ここでは特に$x_i=i$という制限を加えて補間を行いやすくする。

準備としてまず降冪(Falling Power) $x^{\underline{n}}$を $x^{\underline{n}}=x(x-1)\cdots(x-n+1)$、特に $x^{\underline{0}}=1$と定義する。次に差分演算子(Difference Operator) $\Delta$を $\Delta f(x)=f(x+1)-f(x)$と定義する。そうすると、$n>1$ならば

$$\begin{eqnarray*} \Delta x^{\underline{n}}&=&(x+1)x\cdots(x-n+1)-x(x-1)\cdots(x-n+1) \\ &=&nx(x-1)\cdots(x-n+2) \end{eqnarray*}$$

となり、 $\Delta x^{\underline{1}}=1$なので、$n\ge1$に対して

$$\Delta x^{\underline{n}}=nx^{\underline{n-1}}$$ (3.36)

が成り立つ。さて、$x^n$は第2種Stirling数を用いて

$$x^n=\sum_{k=0}^{n}{n \brace k}x^{\underline{k}}$$ (3.37)

と表されるので、$w(x)$を

$$w(x)=\sum_{k=0}^{n-1}c_kx^{\underline{k}}$$ (3.38)

と書き換えたとき、$c_k$は全て負でない。ここで

$$\begin{eqnarray*} \Delta^0 w(x)&=&c_0+c_1x^{\underline{1}}+c_2x^{\underline{2}}+... ...nderline{1}}+\cdots \\ \Delta^3 w(x)&=&6c_3+\cdots \\ &\vdots& \end{eqnarray*}$$

であるから、$c_k$は

$$c_k=\frac{\Delta^k w(0)}{k!}$$ (3.39)

と簡単に求まる。よって、(3.11)において $x^{\underline{k}}$を全て展開すれば$w_k$を決定することができる。

3.4 Toom-Cook法

今まで述べてきたことを使って乗算を行うアルゴリズムをToom-Cook法(Toom-Cook Method)という。例として、 $1093=4\cdot16^2+4\cdot16+5$と $3511=13\cdot16^2+11\cdot16+7$の16進乗算を行ってみる。 $u(x)=4x^2+4x+5,v(x)=13x^2+11x+7,w(x)=u(x)v(x)$とおくと、

$$\begin{array}{r@{}c@{}rr@{}c@{}rr@{}c@{}rr@{}c@{}rr@{}c@{}r} ... ...&w(1)&=&403,&w(2)&=&2349,&w(3)&=&8321,&w(4)&=&22015 \end{array}$$

となる。$w(x)$の補間は次のような表を考えると分かりやすい。まず0段目に$w(i)$の値を並べる。

$$\begin{array}{c@{\quad:\quad}ccccccccc} 0&35&\phantom{368}&4... ...956}&2349&\phantom{5972}&8321&\phantom{13694}&22015 \end{array}$$

次に、0段目の隣り合う要素の差を1段目に並べる。

$$\begin{array}{c@{\quad:\quad}ccccccccc} 0&35& &403& &2349& &8321& &22015 \\ 1& &368& &1946& &5972& &13694& \end{array}$$

次に、1段目の隣り合う要素の差を2で割ったのものを2段目に並べる。

$$\begin{array}{c@{\quad:\quad}ccccccccc} 0&35& &403& &2349& &... ...1946& &5972& &13694& \\ 2& & &789& &2013& &3861& & \end{array}$$

$k-1$段目の隣り合う要素の差を$k$で割ったのものを$k$段目に並べる、という操作を繰り返すと、次の表が得られる。

$$\begin{array}{c@{\quad:\quad}ccccccccc} 0&35& &403& &2349& &... ... 3& & & & 408& & 616& & & \\ 4& & & & & 52& & & & \end{array}$$

このとき、$k$段目の左端の値が $x^{\underline{k}}$の係数となる。つまり、

$$w(x)=52x^{\underline{4}}+408x^{\underline{3}}+789x^{\underline{2}}+368x^{\underline{1}}+35x^{\underline{0}}$$

である。ここで $w(x)=(((52(x-3)+408)(x-2)+789)(x-1)+368)x+35$と書き直すと、$x^k$の係数は次のように計算できることが分かる。

$$\begin{array}{rrrrrr} &52& 408& & & \\ -& & 3\cdot52& & & ... ...cline{1-5} &52& 96&137\phantom{)}&83\phantom{)}&35 \end{array}$$

ここで求まった係数の繰上げを行うと、

$$\begin{eqnarray*} 35&=&2\cdot16+3 \\ 83+2&=&5\cdot16+5 \\ 137+5&=&8\cdot16+14 \\ 96+8&=&6\cdot16+8 \\ 52+6&=&3\cdot16+10 \end{eqnarray*}$$

となるから、

$$w(x)=3x^5+10x^4+8x^3+14x^2+5x+3$$

である。実際に値を代入すると、 $w(16)=3837523=1093\times3511$となるから、うまくいっている。