円周率.jp (http://xn--w6q13e505b.jp/method/toomcook.html)

Toom-Cook 法や FFT を使った多倍長乗算アルゴリズムでは、 基本的に多倍長整数を $k$ ワードずつの長さの短い多倍長整数に分割し計算する。 この時、分割された状態の多倍長数列をただの数値列ではなく 一度多項式の係数として捉えなおす。

$$\begin{eqnarray} && A = A_{d-1}R^{(d-1)k}+A_{d-2}R^{(d-2)k}+\cdots+A_1R^k+A_0\\ &\Longrightarrow& A(x) = A_{d-1}x^{d-1}+A_{d-2}x^{d-2}+\cdots+A_1x+A_0 \end{eqnarray}$$

すると元の数値とは

$$A(R^k) = A$$

という関係が成り立つ。

この多項式の形で考えると少し別の計算路が見えてくる。 ちなみにこの時 $A=A(R^k)$、$B=B(R^k)$ が成り立つ必要があるので $k$ は共通化しておかなければならない。

結果、その多項式での積 $C(x)=A(x)B(x)$ を求めることができれば $C(R^k)$ を計算(正規化)することで本来求めたかった積 $C=A\times B$ を求めることができる、という流れである。

Karatsuba 法では長桁の数を 2 分割して重複する乗算部分を作ることで全体の演算量を削減した。 Toom-Cook 法では上記のような流れでより多くの多倍長整数へと分割し、 多項式の形でアルゴリズムを構築する。

Toom-Cook 法では多項式 $A(x)$、$B(x)$ の次数をそれぞれ $d_A-1$ 次、$d_B-1$ 次とするよう設定する。 計算中、分割後の多倍長整数のワード数 $k$ は条件を満たすよう求められる。 すると $C(x)$ は $d_A+d_B-2$ 次(定数値)の多項式ということになるので $x$ に具体的な数値を代入した $C(x_i)=A(x_i)B(x_i)$ を $d_A+d_B-1$ 組以上求めることができれば $C(x)$ を復元することができる。 そこで $\{x_i\}$ として絶対値の小さな整数を用い、各 $C(x_i)$ を計算してから $C(x)$、$C$ を復元する。 以下、具体例を用いつつその手順を紹介する。

ここでは $d=d_A=d_B=3$ とした場合の例を提示しながら計算手順を示す。

手順 1. 多項式→値の変換

まず $A(x)$ と $B(x)$ の $x$ に $\{-d+1, -d+2, \cdots, d-1\}$ を代入した値を求める。

$$\left\{ \begin{eqnarray} A(-2) &=& 4A_2 - 2A_1 + A_0\\ A(-1) &=& A_2 - A_1 + A_0\\ A(0) &=& A_0\\ A(1) &=& A_2 + A_1 + A_0\\ A(2) &=& 4A_2 + 2A_1 + A_0 \end{eqnarray} \right.$$

演算としては「単精度数と多倍長数の乗算」と「多倍長数同士の加減算」 のみで構成されるので、計算量は $2(2d-1)(d-1)A(n/d)$ 程度になる。 実際には $x=0$ で計算が実質的に必要無かったり、$x=\pm1$ で単精度数との乗算が必要無かったり、 $\pm x$ 代入での単精度数との乗算が共通化できたりするなどの演算量削減が可能であるが、 ここではそれを無視して計算量を求めている。

手順 2. 値での乗算

次に、$x$ に同じ値を代入した $A(x)$ と $B(x)$ の積 $C(x)=A(x)B(x)$ を求める。

$$\left\{ \begin{eqnarray} C(-2) &=& A(-2) \times B(-2)\\ C(-1) &=& A(-1) \times B(-1)\\ C(0) &=& A(0) \times B(0)\\ C(1) &=& A(1) \times B(1)\\ C(2) &=& A(2) \times B(2) \end{eqnarray} \right.$$

多倍長整数同士の乗算が複数回行われているが、 $n/d$ ワード同士の乗算が$2d-1$ 回なので 全体の演算量は $(2d-1)M(n/d)$ になる。

手順 3. 畳込み乗算への逆変換

ここで $C(x)$ の多項式表現を知っている前提で $\{x_i\}$ を代入する計算を考えると

$$\begin{pmatrix}C(-2)\\C(-1)\\C(0)\\C(1)\\C(2)\end{pmatrix} = \begin{pmatrix} 1 & -2 & 4 & -8 & 16\\ 1 & -1 & 1 & -1 & 1\\ 1 & 0 & 0 & 0 & 0\\ 1 & 1 & 1 & 1 & 1\\ 1 & 2 & 4 & 8 & 16\\ \end{pmatrix} \begin{pmatrix}C_0\\C_1\\C_2\\C_3\\C_4\end{pmatrix}$$

という形になる。 しかし実際には逆に代入後の値 $C(x_i)$ が分かっていて係数 $C_i$ を求めたいので、この連立方程式を解くことになる。

$$\begin{eqnarray} \begin{pmatrix}C_0\\C_1\\C_2\\C_3\\C_4\end{pmatrix} &=& \begin{pmatrix} 1 & -2 & 4 & -8 & 16\\ 1 & -1 & 1 & -1 & 1\\ 1 & 0 & 0 & 0 & 0\\ 1 & 1 & 1 & 1 & 1\\ 1 & 2 & 4 & 8 & 16\\ \end{pmatrix}^{-1} \begin{pmatrix}C(-2)\\C(-1)\\C(0)\\C(1)\\C(2)\end{pmatrix} \\ &=& \dfrac{1}{24} \begin{pmatrix} 0 & 0 & 24 & 0 & 0\\ 2 & -16 & 0 & 16 & -2\\ -1 & 16 & -30 & 16 & -1\\ -2 & 4 & 0 & -4 & 2\\ 1 & -4 & 6 & -4 & 1 \end{pmatrix} \begin{pmatrix}C(-2)\\C(-1)\\C(0)\\C(1)\\C(2)\end{pmatrix} \\ \end{eqnarray}$$

この逆演算の式は $d$ によって決まるので プログラムに組み込まれることが前提になっている。 また、$x$ の値は十分に小さなものを設定したので 演算に必要な数値も単精度で収まることが期待される。

手順 4.

$$C = C(R^k) = C_4R^{4k} + C_3R^{3k} + C_2R^{2k} + C_1R^k + C_0$$

全体を繋げて多倍長数としての正規化を行う。 この時、分割された各値 $\{C_i\}$ は約 $2n/d$ ワードであり 繋げる際にずらすワード数は $k\simeq n/d$ なので、 $(2d-2)A(n/d)$ の演算を要する。

Toom-Cook 法の計算量を求める。 上記計算手順に書いた必要計算量を足し合わせると

$$\begin{eqnarray} M(n) &=& 2(2d-1)(d-1)A(n/d) + (2d-1)M(n/d) + (2d-2)(d-1)A(2n/d) + (2d-2)A(n/d)\\ &=& 2(4d-3)(d-1)A(n/d) + (2d-1)M(n/d) \end{eqnarray}$$

となる。ここで加減算のコストは $A(n)=O(n)$ であることを利用している。

多倍長同士の演算については、加減算よりも乗算のコストが非常に高い、 つまり $A(n) \ll M(n)$ という前提に立つと、

$$M(n) = O(n^{\log_d(2d-1)})$$

となる。$d=2$ の時に Karatsuba 法と同じ計算量になり、 $d$ が大きくなれば $O(n^{1+\varepsilon})$ へ漸近的に近付く。 しかし、$d$ が大きな値になると $A(n) \ll M(n)$ という前提が成り立たない、 $A(n/d)$ の係数が大きい、上記手順内の係数が単精度で収まらない、 等の理由でコストが増えるため、$d$ を大きくすれば速い というわけではないことに注意していただきたい。

分割数を変更することで多倍長乗算の回数も変動するが、 単精度数との乗算もなるべく少なく、小さい数との演算に制限したい。 そこで上の説明では多倍長整数を 1 変数 $x$ の多項式

$$A(x) = \sum_{k=0}^{d_A-1} A_kx^k$$

に変換したが、ここでは 2 変数 $(x,y)$ の斉次多項式

$$A(x, y) = \sum_{k=0}^{d-1} A_kx^ky^{d-1-k}$$

に変換する。 代入する値が 2 次元に広がることで同程度に絶対値が小さな値の候補が増える。 ただし代入する数の組は代入後の数式が一次独立になるようにしなければならない。 分かりやすい例としては $(x,y)=(1,0)$ を代入することで

$$C(1,0) = A(1,0) \times B(1,0) = A_{d-1} B_{d-1}$$

と、$n/d$ ワードの乗算 1 回で 1 つの数値を得ることができる。

実際の演算では掛け合わせる多倍長数の長さが同程度になるとは限らないため、 分割数によっては分割後の一部の "数" が 0 になることもあり、 常に同じ分割数を利用していると無駄な計算をするかもしれない。 また、計算に用いた係数は分割数にしか依存しないためプログラムでは ハードコーディングで行うことになる。 そのためいくつかの分割数のパターンを実装することで 0 との乗算を省略したり 代入する数を小さくして演算を簡単にすることができる。 ここではその具体的な例をいくつか挙げる。

2 分割 × 2 分割 (Karatsuba 法相当)

$$\begin{pmatrix} A(1,0) & B(1,0)\\ A(0,1) & B(0,1)\\ A(-1,1) & B(-1,1) \end{pmatrix} = \begin{pmatrix} 0 & 1\\ 1 & 0\\ 1 & -1\\ \end{pmatrix} \begin{pmatrix} A_0 & B_0\\ A_1 & B_1 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2 \end{pmatrix} = \begin{pmatrix} 0 & 1 & 0\\ 1 & 1 & -1\\ 1 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1) \end{pmatrix}$$

3 分割 × 2 分割

$$\begin{pmatrix} A(1,0)\\ A(0,1)\\ A(-1,1)\\ A(1,1) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 1 & -1 & 1\\ 1 & 1 & 1\\ \end{pmatrix} \begin{pmatrix} A_0\\ A_1\\ A_2 \end{pmatrix}$$ $$\begin{pmatrix} B(1,0)\\ B(0,1)\\ B(-1,1)\\ B(1,1) \end{pmatrix} = \begin{pmatrix} 0 & 1\\ 1 & 0\\ 1 & -1\\ 1 & 1\\ \end{pmatrix} \begin{pmatrix} B_0\\ B_1 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2\\ C_3 \end{pmatrix} = \frac{1}{2} \begin{pmatrix} 0 & 2 & 0 & 0\\ -2 & 0 & -1 & 1\\ 0 & -2 & 1 & 1\\ 2 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1)\\ C(1,1) \end{pmatrix}$$

3 分割 × 3 分割

$$\begin{pmatrix} A(1,0) & B(1,0)\\ A(0,1) & B(0,1)\\ A(-1,1) & B(-1,1)\\ A(1,1) & B(1,1)\\ A(-2,1) & B(-2,1) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 1 & -1 & 1\\ 1 & 1 & 1\\ 1 & -2 & 4\\ \end{pmatrix} \begin{pmatrix} A_0 & B_0\\ A_1 & B_1\\ A_2 & B_2 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2\\ C_3\\ C_4 \end{pmatrix} = \frac{1}{6} \begin{pmatrix} 0 & 6 & 0 & 0 & 0\\ -12 & 3 & -6 & 2 & 1\\ -6 & -6 & 3 & 3 & 0\\ 12 & -3 & 3 & 1 & -1\\ 6 & 0 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1)\\ C(1,1)\\ C(-2,1) \end{pmatrix}$$

4 分割 × 2 分割

$$\begin{pmatrix} A(1,0)\\ A(0,1)\\ A(-1,1)\\ A(1,1)\\ A(-2,1) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & -1 & 1 & -1\\ 1 & 1 & 1 & 1\\ 1 & -2 & 4 & -8\\ \end{pmatrix} \begin{pmatrix} A_0\\ A_1\\ A_2\\ A_3 \end{pmatrix}$$ $$\begin{pmatrix} B(1,0)\\ B(0,1)\\ B(-1,1)\\ B(1,1)\\ B(-2,1) \end{pmatrix} = \begin{pmatrix} 0 & 1\\ 1 & 0\\ 1 & -1\\ 1 & 1\\ 1 & -2\\ \end{pmatrix} \begin{pmatrix} B_0\\ B_1 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2\\ C_3\\ C_4 \end{pmatrix} = \frac{1}{6} \begin{pmatrix} 0 & 6 & 0 & 0 & 0\\ -12 & 3 & -6 & 2 & 1\\ -6 & -6 & 3 & 3 & 0\\ 12 & -3 & 3 & 1 & -1\\ 6 & 0 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1)\\ C(1,1)\\ C(-2,1) \end{pmatrix}$$

4 分割 × 3 分割

$$\begin{pmatrix} A(1,0)\\ A(0,1)\\ A(-1,1)\\ A(1,1)\\ A(-2,1)\\ A(2,1) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & -1 & 1 & -1\\ 1 & 1 & 1 & 1\\ 1 & -2 & 4 & -8\\ 1 & 2 & 4 & 8\\ \end{pmatrix} \begin{pmatrix} A_0\\ A_1\\ A_2\\ A_3 \end{pmatrix}$$ $$\begin{pmatrix} B(1,0)\\ B(0,1)\\ B(-1,1)\\ B(1,1)\\ B(-2,1)\\ B(2,1) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1\\ 1 & 0 & 0\\ 1 & -1 & 1\\ 1 & 1 & 1\\ 1 & -2 & 4\\ 1 & 2 & 4\\ \end{pmatrix} \begin{pmatrix} B_0\\ B_1\\ B_2 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2\\ C_3\\ C_4\\ C_5 \end{pmatrix} = \frac{1}{24} \begin{pmatrix} 0 & 24 & 0 & 0 & 0 & 0\\ 96 & 0 & -16 & 16 & 2 & -2\\ 0 & -30 & 16 & 16 & -1 & -1\\ -120 & 0 & 4 & -4 & -2 & 2\\ 0 & 6 & -4 & -4 & 1 & 1\\ 24 & 0 & 0 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1)\\ C(1,1)\\ C(-2,1)\\ C(2,1) \end{pmatrix}$$

4 分割 × 4 分割

$$\begin{pmatrix} A(1,0) & B(1,0)\\ A(0,1) & B(0,1)\\ A(-1,1) & B(-1,1)\\ A(1,1) & B(1,1)\\ A(-2,1) & B(-2,1)\\ A(2,1) & B(2,1)\\ A(-1,2) & B(-1,2) \end{pmatrix} = \begin{pmatrix} 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\\ 1 & -1 & 1 & -1\\ 1 & 1 & 1 & 1\\ 1 & -2 & 4 & -8\\ 1 & 2 & 4 & 8\\ 8 & -4 & 2 & -1\\ \end{pmatrix} \begin{pmatrix} A_0 & B_0\\ A_1 & B_1\\ A_2 & B_2\\ A_3 & B_3 \end{pmatrix}$$ $$\begin{pmatrix} C_0\\ C_1\\ C_2\\ C_3\\ C_4\\ C_5\\ C_6 \end{pmatrix} = \frac{1}{360} \begin{pmatrix} 0 & 360 & 0 & 0 & 0 & 0 & 0\\ 720 & 720 & 240 & 80 & -10 & -6 & -16\\ 1440 & -450 & 240 & 240 & -15 & -15 & 0\\ -900 & -900 & -540 & 140 & 20 & 0 & 20\\ -1800 & 90 & -60 & -60 & 15 & 15 & 0\\ 180 & 180 & 120 & -40 & -10 & 6 & -4\\ 360 & 0 & 0 & 0 & 0 & 0 & 0\\ \end{pmatrix} \begin{pmatrix} C(1,0)\\ C(0,1)\\ C(-1,1)\\ C(1,1)\\ C(-2,1)\\ C(2,1)\\ C(-1,2) \end{pmatrix}$$