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

Toom-Cook 法

Karatsuba 法では長桁の数を 2 分割して重複する乗算部分を作ることで全体の演算量を削減した. Toom-Cook 法では長桁の数を $d$ 分割してこれと同様の演算を行う.

$d$ 分割した多倍長数を $d-1$ 次多項式として捉えると, 2 多項式 $A(x)$,$B(x)$ の積は $2d-2$ 次の多項式になる. Toom-Cook 法を始めとした多くのアルゴリズムでは $R$ が変数の多項式 $x$ であると考えて計算を行い, 多項式の積に基底 $R$ を代入して多倍長に戻す, という考え方で計算する.

\[ \left\{ \begin{eqnarray} A &=& a_{d-1}R^{d-1} + a_{d-2}R^{d-2} + \cdots + a_0\\ B &=& b_{d-1}R^{d-1} + b_{d-2}R^{d-2} + \cdots + b_0\\ A\cdot B = C &=& c_{2d-2}R^{2d-2} + c_{2d-3}R^{2d-3} + \cdots + c_0\\ \end{eqnarray}\right. \] \[ \Longleftrightarrow \left\{ \begin{eqnarray} A(x) &=& a_{d-1}x^{d-1} + a_{d-2}x^{d-2} + \cdots + a_0\\ B(x) &=& b_{d-1}x^{d-1} + b_{d-2}x^{d-2} + \cdots + b_0\\ A(x)\cdot B(x) = C(x) &=& c_{2d-2}x^{2d-2} + c_{2d-3}x^{2d-3} + \cdots + c_0\\ \end{eqnarray}\right. \]

手順

ここでは $d=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=1$ で単精度数との乗算が必要無かったり, $\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)$ の多項式表現を知っている前提で $A(x)$、$B(x)$ と同じ値を代入したとすると

\[ \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} \]

という形で構成されているはずである。 一方でその具体的な値は手順 2 で分かっているので、 この連立方程式を解くことで $C(x)$ の係数が分かる。

\[ \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$ によって決まるので, プログラムに組み込んでおくことで高速に実行することができる. さらに,同じ演算になっている部分の処理を共通化することで さらに高速な実装を組むことができる.

手順 4.

\[ C = C(R) = c_4R^4 + c_3R_3 + c_2R^2 + c_1R + c_0 \]

全体を繋げて多倍長数としての正規化を行う. この時,分割された各値 $\{c_i\}$ は約 $2N/d$ 桁あり, $R$ (繋げる際にずらす桁数) は $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$ を大きくすれば速い というわけではないことに注意していただきたい.

定番の改良

分割数を変更することで多倍長乗算の回数を減少させるが, 多倍長加減算の回数も実際的な計算速度には影響を与える. GMP などの実装では $x=d-1$ を代入する計算代わりに(形式的に) $x=\infty$ を代入し,

\[ C(\infty) = A(\infty) \times B(\infty) = a_{d-1} b_{d-1} \]

で計算することで加減算の数を減らしている. 実際の計算の過程においては何を代入しているか はあまり重要ではなく,手順 3 における係数変換行列, そしてその逆行列 (例示している変換式) の構築が容易かどうかだけの問題なので, $\infty$ を代入しても問題無い.

定数例

掛け合わせる多倍長数の桁数が同程度になるとは限らないため, 分割数によっては分割後の一部の "数" が 0 になることもある. また、計算に用いた係数は分割数にしか依存しないため、プログラムでは ハードコーディングで行うことになる. そのため,0 との乗算については省略するルーチンを作ったり、 仮代入する数を小さくすることで演算を簡単にするようになる.

上の説明では 1 変数 $x$ の多項式

\[ A(x) = \sum_{k=0}^{d-1} a_kx^k \]

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

\[ A(x, y) = \sum_{k=0}^{d-1} a_kx^ky^{d-1-k} \]

として考えることでより効率的に式を組み立てることができる。 ここではその具体的な例をいくつか挙げる.

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} \]