Toom-Cook 乗算法

概説

多倍長整数同士の乗算において、被乗数・乗数を ある程度大きな桁数(基数$x$) で分割した $l$次・$m$次多項式 とみなし、多項式同士の乗算を行う。（改めて多項式乗算の積に基数$x$を代入し、繰り上がりの正規化を行えば多倍長整数の積が求められる）

この時、多項式係数同士の乗算回数を $(l+m+1)$ に抑える方法を検討する。(基数$x$が充分に大きければ、伴って発生する加減算・定数乗除算の計算量は少ないとみなす)

$\{A_i\},\{B_j\},\{C_k\}$ をそれぞれ被乗数・乗数・多項式乗算の積、の係数列とする。

$$\left( A_i = \begin{cases} \mbox{(any)} & (0 \le i \le l) \\ 0 & (i \lt 0, l \lt i) \end{cases},\ B_j = \begin{cases} \mbox{(any)} & (0 \le j \le m) \\ 0 & (j \lt 0, m \lt j) \end{cases},\ C_k = \sum_{i=0}^{k} A_{i}B_{k-i} \right)$$

$$\begin{align*} &\ \quad(A_lx^l+A_{l-1}x^{l-1}+\cdots +A_1x+A_0)(B_mx^m+B_{m-1}x^{m-1}+\cdots +B_1x+B_0)\\ &=A_lB_mx^{l+m}+(A_lB_{m-1}+A_{l-1}B_m)x^{l+m-1}+\cdots +(A_1B_0+A_0B_1)x+A_0B_0\\ &=C_{l+m}x^{l+m}+C_{l+m-1}x^{l+m-1}+\cdots +C_1x+C_0\\ &=\Biggl(\sum_{i=0}^lA_ix^i\Biggr)\Biggl(\sum_{j=0}^mB_jx^j\Biggr)=\sum_{k=0}^{l+m}\Biggl(\sum_{i+j=k}A_iB_j\Biggr)x^k=\sum_{k=0}^{l+m}\Biggl(\sum_{i=0}^{k}A_{i}B_{k-i}\Biggr)x^k=\sum_{k=0}^{l+m}C_kx^k \end{align*}$$

$(\sum_{i=0}^lA_ix^i),(\sum_{j=0}^mB_jx^j),(\sum_{k=0}^{l+m}C_kx^k)$ に整数の組 $(n,d)$ を用いて $x=n/d$ と代入し、それぞれ $d^l,d^m,d^{l+m}$ を乗じた値 $\mathcal{A}(n/d),\mathcal{B}(n/d),\mathcal{C}(n/d)$ を考えると（ $n/d$ は整数比を示す為の便宜上の表記、実質的には2変数の関数? $1/0$ は $\infty$ と表記 ）

$$\mathcal{A}(n/d) = d^l \sum_{i=0}^l A_i (n/d)^i = \sum_{i=0}^l A_i n^i d^{l-i} = A_l n^l + A_{l-1} n^{l-1} d + \cdots + A_i n^i d^{l-i} +\cdots+ A_1 n d^{l-1} + A_0 d^l$$

$$\mathcal{B}(n/d) = d^m \sum_{j=0}^m B_j (n/d)^j = \sum_{j=0}^m B_j n^j d^{m-j} = B_m n^m + B_{m-1} n^{m-1} d + \cdots + B_j n^j d^{m-j} + \cdots + B_1 n d^{m-1} + B_0 d^m$$

$$\mathcal{C}(n/d) = d^{l+m} \sum_{k=0}^{l+m} C_k (n/d)^k = \sum_{k=0}^{l+m} C_k n^k d^{l+m-k} = \left( d^l \sum_{i=0}^l A_i (n/d)^i \right) \left( d^m \sum_{j=0}^m B_j (n/d)^j \right) = \mathcal{A}(n/d) \times \mathcal{B}(n/d)$$

$(l+m+1)$ 個の異なる整数比の組 $\{(n_0,d_0),(n_1,d_1),\dots,(n_{l+m},d_{l+m})\}$ $[\forall p,q \in \{0,1,\dots,l+m\} ; p \neq q \to n_pd_q \neq n_qd_p]$ を用意し、下記のように逆行列を求める。

$$\begin{bmatrix}\mathcal{A}(n_{l+m}/d_{l+m})\\\vdots\\\mathcal{A}(n_1/d_1)\\\mathcal{A}(n_0/d_0)\end{bmatrix} =\begin{bmatrix}(n_{l+m})^l&(n_{l+m})^{l-1}(d_{l+m})&\cdots&(n_{l+m})(d_{l+m})^{l-1}&(d_{l+m})^l\\\vdots&\vdots&\ddots&\vdots&\vdots\\(n_1)^l&(n_1)^{l-1}(d_1)&\cdots&(n_1)(d_1)^{l-1}&(d_1)^l\\(n_0)^l&(n_0)^{l-1}(d_0)&\cdots&(n_0)(d_0)^{l-1}&(d_0)^l\end{bmatrix} \begin{bmatrix}A_{l}\\\vdots\\A_1\\A_0\end{bmatrix}\tag{1}$$

$$\begin{bmatrix}\mathcal{B}(n_{l+m}/d_{l+m})\\\vdots\\\mathcal{B}(n_1/d_1)\\\mathcal{B}(n_0/d_0)\end{bmatrix} =\begin{bmatrix}(n_{l+m})^m&(n_{l+m})^{m-1}(d_{l+m})&\cdots&(n_{l+m})(d_{l+m})^{m-1}&(d_{l+m})^m\\\vdots&\vdots&\ddots&\vdots&\vdots\\(n_1)^m&(n_1)^{m-1}(d_1)&\cdots&(n_1)(d_1)^{m-1}&(d_1)^m\\(n_0)^m&(n_0)^{m-1}(d_0)&\cdots&(n_0)(d_0)^{m-1}&(d_0)^m\end{bmatrix} \begin{bmatrix}B_{m}\\\vdots\\B_1\\B_0\end{bmatrix}\tag{2}$$

$$\begin{bmatrix}\mathcal{A}(n_{l+m}/d_{l+m})\times\mathcal{B}(n_{l+m}/d_{l+m})\\\vdots\\\mathcal{A}(n_1/d_1)\times\mathcal{B}(n_1/d_1)\\\mathcal{A}(n_0/d_0)\times\mathcal{B}(n_0/d_0)\end{bmatrix} =\begin{bmatrix}(n_{l+m})^{l+m}&(n_{l+m})^{l+m-1}(d_{l+m})&\cdots&(n_{l+m})(d_{l+m})^{l+m-1}&(d_{l+m})^{l+m}\\\vdots&\vdots&\ddots&\vdots&\vdots\\(n_1)^{l+m}&(n_1)^{l+m-1}(d_1)&\cdots&(n_1)(d_1)^{l+m-1}&(d_1)^{l+m}\\(n_0)^{l+m}&(n_0)^{l+m-1}(d_0)&\cdots&(n_0)(d_0)^{l+m-1}&(d_0)^{l+m}\end{bmatrix} \begin{bmatrix}C_{l+m}\\\vdots\\C_1\\C_0\end{bmatrix}$$

$$\begin{bmatrix}C_{l+m}\\\vdots\\C_1\\C_0\end{bmatrix} =\left(\begin{bmatrix}(n_{l+m})^{l+m}&(n_{l+m})^{l+m-1}(d_{l+m})&\cdots&(n_{l+m})(d_{l+m})^{l+m-1}&(d_{l+m})^{l+m}\\\vdots&\vdots&\ddots&\vdots&\vdots\\(n_1)^{l+m}&(n_1)^{l+m-1}(d_1)&\cdots&(n_1)(d_1)^{l+m-1}&(d_1)^{l+m}\\(n_0)^{l+m}&(n_0)^{l+m-1}(d_0)&\cdots&(n_0)(d_0)^{l+m-1}&(d_0)^{l+m}\end{bmatrix}^{-1}\right) \begin{bmatrix}\mathcal{A}(n_{l+m}/d_{l+m})\times\mathcal{B}(n_{l+m}/d_{l+m})\\\vdots\\\mathcal{A}(n_1/d_1)\times\mathcal{B}(n_1/d_1)\\\mathcal{A}(n_0/d_0)\times\mathcal{B}(n_0/d_0)\end{bmatrix}\tag{3}$$

よって、 $\{(n_0,d_0),(n_1,d_1),\dots,(n_{l+m},d_{l+m})\}$ の組から逆行列をあらかじめ求めておき、式$(1),(2),(3)$の計算をそれぞれ行うと、$\{A_i\},\{B_j\}$から$\{C_k\}$を求められる事になる。

参照

• http://円周率.jp/method/toomcook.html 註:4分割×4分割の係数例に誤記あり
• http://円周率.jp/method/mult.html 乗算アルゴリズムの概説
• https://en.wikipedia.org/wiki/Toom–Cook_multiplication Toom-Cook乗算法
• https://ja.wikipedia.org/wiki/行列 行列の記法・定義、行列の積
• https://ja.wikipedia.org/wiki/正則行列 逆行列
• https://ja.wikipedia.org/wiki/カラツバ法 Karatsuba法
• https://ja.wikipedia.org/wiki/畳み込み 多項式乗算と線形畳み込みの関係
• https://ja.wikipedia.org/wiki/ランダウの記号 O-記法
• https://ja.wikipedia.org/wiki/高速フーリエ変換 FFT
• https://en.wikipedia.org/wiki/Schönhage–Strassen_algorithm FFTを用いた高速乗算法
• http://www.bodrato.it/software/toom.html GMPへのToom-Cook実装
• https://gmplib.org/manual/Multiplication-Algorithms.html GMPで用いられる乗算法の解説
• https://gmplib.org/devel/thres/ GMP開発版 算法切り替えしきい値
• ./data/toom_mat.json Toom-Cook係数例データ(JSON形式)

再帰的適用による漸近的な計算コスト

おおむね、上の物ほど小さな桁数での計算コストが小さく、下の物ほど大きな桁数での計算コストが小さくなる。

• Schoolbook(筆算): $O(n^2)$ (桁数1000倍で約1000000倍相当の時間計算量)
• Toom-2(Karatsuba): $O(n^{\log(3)/\log(2)}) \simeq O(n^{1.585})$ (桁数1000倍で約56871倍相当の時間計算量)
• Toom-3: $O(n^{\log(5)/\log(3)}) \simeq O(n^{1.465})$ (桁数1000倍で約24827倍相当の時間計算量)
• Toom-4: $O(n^{\log(7)/\log(4)}) \simeq O(n^{1.404})$ (桁数1000倍で約16257倍相当の時間計算量)
• Toom-6: $O(n^{\log(11)/\log(6)}) \simeq O(n^{1.338})$ (桁数1000倍で約10348倍相当の時間計算量)
• Toom-8: $O(n^{\log(15)/\log(8)}) \simeq O(n^{1.302})$ (桁数1000倍で約8070倍相当の時間計算量)
• Toom-10: $O(n^{\log(19)/\log(10)}) \simeq O(n^{1.279})$ (桁数1000倍で約6859倍相当の時間計算量)
• Toom-12: $O(n^{\log(23)/\log(12)}) \simeq O(n^{1.262})$ (桁数1000倍で約6102倍相当の時間計算量)
• Toom-16: $O(n^{\log(31)/\log(16)}) \simeq O(n^{1.239})$ (桁数1000倍で約5196倍相当の時間計算量)
• Toom-20: $O(n^{\log(39)/\log(20)}) \simeq O(n^{1.223})$ (桁数1000倍で約4664倍相当の時間計算量)
• Toom-24: $O(n^{\log(47)/\log(24)}) \simeq O(n^{1.211})$ (桁数1000倍で約4310倍相当の時間計算量)
• Toom-36: $O(n^{\log(71)/\log(36)}) \simeq O(n^{1.190})$ (桁数1000倍で約3703倍相当の時間計算量)
• Toom-44: $O(n^{\log(87)/\log(44)}) \simeq O(n^{1.180})$ (桁数1000倍で約3471倍相当の時間計算量)
• Toom-56: $O(n^{\log(111)/\log(56)}) \simeq O(n^{1.170})$ (桁数1000倍で約3235倍相当の時間計算量)
• Toom-64: $O(n^{\log(127)/\log(64)}) \simeq O(n^{1.165})$ (桁数1000倍で約3121倍相当の時間計算量)
• FFT(Schönhage–Strassen): $O(n \log n \log\log n)$ (十分に大きな桁数[10進数で3万桁〜15万桁程度以上?]ではToom-Cook法より少ない時間計算量に)

係数例