多倍長演算について

近藤健一

目次

1  基本的なアルゴリズム  2

1.1 多倍長演算とは . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 多倍長数の表現と加減乗算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 除算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.4 平方根 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.5 基数変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2  Karatsuba 系列のアルゴリズム  6

2.1 Karatsuba 乗算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2 Karatsuba 除算 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.3 Karatsuba 開平 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4 Karatsuba 基数変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3  Karatsuba 乗算の一般化  10

3.1 数と多項式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 連立方程式の解としての畳み込み . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.3 多項式補間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.4 Toom-Cook 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

4  高速 Fourier 変換を利用した乗算法  14

4.1 多項式の乗算の行列による表現 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 4.2 高速 Fourier 変換 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 4.3 Cooley-Tukey のアルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5  用語について  17

1 基本的なアルゴリズム

1.1 多倍長演算とは

CPU が直接扱える数の範囲には限りがあり、昔からのプログラミング言語ではその制限が残っていることが多い。 C 言語を例にとってみると、 int 型の変数は通常 CPU のレジスタのサイズとなるので、 32bit の CPU では 2^{32} = 4294967296 種類の数しか表現することができないことが多い。それに対して、 Lisp 系の言語や最近のスクリプト言語では任意の桁数の数を扱えるものがよくある。これらの言語は、限られた桁数の数のみによる任意の桁数の数の計算、つまり多倍長演算 (Multiple Precision Arithmetic) を行うプログラムを作って CPU の制限を解消しているのである。

1.2 多倍長数の表現と加減乗算

多倍長演算の対象となる数、つまり CPU が直接扱える範囲を超えるような数を多倍長数 (Multiple Precision Number) と呼ぶことにする。多倍長演算を行うにはまず多倍長数を表現しなければならないが、大きな数を日常的に 3141592 と一桁の数の並びで表すように、一桁ごとに変数を割り当てるのが最も単純な方法である。つまり

u_{6} = 3, u_{5} = 1, u_{4} = 4, …, u_{0} = 2

のようにすればよいのである。実装する際には 10 進数よりも 2^{32} 進数などで表現した方が、桁上げがビット演算で行えるのでよい。

多倍長数をこのように表現しておけば、加減乗算は筆算 (Schoolbook Method) を使うことで簡単に実現できるので説明は必要ないだろう。 n 桁同士の加減乗算に必要な時間は、それぞれ Θ(n), Θ(n), Θ(n^{2}) である。除算については次の節で詳しく見ることにする。

1.3 除算

実際にやってみれば分かるように、筆算では (n + 1) 桁と n 桁の除算が重要となってくる。そこで u, v

u = (u_{n}u_{n−1}u_{n−2} … u_{0})_{b}, v = (v_{n−1}v_{n−2} … v_{0})_{b}

に対して u = q_{0}v + r (q_{0} < b, r < v) となる場合の除算を考える。 q_{0} の満たす必要条件は

q_{0}v ≤ u ⟹ q_{0}v_{n−1}b^{n−1} ≤ u_{n}b^{n} + u_{n−1}b^{n−1} ⟺ q_{0} ≤ ⌊\frac{u_{n}b + u_{n−1}}{v_{n−1}} ⌋  (1.1)

となる。よって、 q_{0} ≤ b − 1 であることを考えて

q_{0} = min (⌊\frac{u_{n}b + u_{n−1}}{v_{n−1}} ⌋, b − 1) , r = u − qv  (1.2)

とおき、 r ≥ 0 となるまで

q_{0} := q_{0} − 1, r := r + v  (1.3)

を繰り返せば q_{0}, r が求まる。この繰り返しの回数は、

q′_{0} = min (⌊\frac{u_{n}b + u_{n−1}}{v_{n−1}} ⌋, b − 1)

とおけば q′_{0} − q_{0} で与えられる。ここで

q′_{0} ≤ \frac{u_{n}b + u_{n−1}}{v_{n−1}} ≤ \frac{u}{v_{n−1}b^{n−1}} , q_{0} > \frac{u}{v} − 1

であるから、

q′_{0} − q_{0} < \frac{u}{v_{n−1}b^{n−1}} − \frac{u}{v} + 1 = \frac{u}{v} · \frac{v − v_{n−1}b^{n−1}}{v_{n−1}b^{n−1}} + 1
≤ \frac{u}{v} · \frac{1}{v_{n−1}} + 1
< \frac{b}{v_{n−1}} + 1  (1.4)

となる。よって、 v_{n−1} ≥ b/2 ならば繰り返しは必ず 2 回以下となるのである。この条件を満たすためには適当な数 k を v に掛ければ良い。 u にも k を掛けて除算を行えば、商は q、余りは kr となるので、 kr を k で割れば求める結果が得られる。 b が 2 の冪の場合は k も 2 の冪となるので、シフト演算を用いることでより高速に行える。 2n 桁と n 桁の除算の実行時間は Θ(n^{2}) である。

1.4 平方根

平方根の計算方法を調べるために、 u = (u_{2n−1}u_{2n−2} … u_{0})_{b} の平方根 s = (s_{n−1}s_{n−2} … s_{0})_{b} を求めることを考える。まず、 s_{n−1} は

(s_{n−1}00 … 0)_{b}^{2} ≤ u ⟺ s_{n−1}^{2} ≤ u_{2n−1}b + u_{2n−2}

を満たす最大の整数であるから、

s_{n−1} = ⌊\sqrt{u_{2n−1}b + u_{2n−2}}⌋  (1.5)

と計算できる。ここで誤差 (余り) r_{n−1} を

r_{n−1} = (u_{2n−1}b + u_{2n−2}) − s_{n−1}^{2}  (1.6)

とおく。次に、 s_{n−2} は

(s_{n−1}s_{n−2}0 … 0)_{b}^{2} ≤ u
⟺ (s_{n−1}s_{n−2})_{b}^{2} ≤ (u_{2n−1}u_{2n−2}u_{2n−3}u_{2n−4})_{b}
⟺ (s_{n−1}s_{n−2})_{b}^{2} − (s_{n−1}0)_{b}^{2} ≤ (u_{2n−1}u_{2n−2}u_{2n−3}u_{2n−4})_{b} − (s_{n−1}0)_{b}^{2}
⟺ (2s_{n−1}b + s_{n−2})s_{n−2} ≤ r_{n−1}b^{2} + u_{2n−3}b + u_{2n−4}

を満たす最大の整数であり、誤差 r_{n−2} は右辺から左辺を引いた

r_{n−2} = r_{n−1}b^{2} + u_{2n−3}b + u_{2n−4} − (2s_{n−1}b + s_{n−2})s_{n−2}

となる。同様にして、一般の s_{k} は

(2(s_{n−1}s_{n−2} … s_{k+1})_{b}b + s_{k})s_{k} ≤ r_{k+1}b^{2} + u_{2k+1}b + u_{2k}  (1.7)

を満たす最大の整数で、誤差 r_{k} は

r_{k} = r_{k+1}b^{2} + u_{2k+1}b + u_{2k} − (2(s_{n−1}s_{n−2} … s_{k+1})_{b}b + s_{k})s_{k}  (1.8)

となることが分かる。 s_{k}, r_{k} を求める方法を考えよう。 (1.7) から

(2(s_{n−1}s_{n−2} … s_{k+1})_{b}b + s_{k})s_{k} ≤ r_{k+1}b^{2} + u_{2k+1}b + u_{2k}
⟹ 2(s_{n−1}s_{n−2} … s_{k+1})_{b}bs_{k} ≤ r_{k+1}b^{2} + u_{2k+1}b + u_{2k}
⟺ 2(s_{n−1}s_{n−2} … s_{k+1})_{b}s_{k} ≤ r_{k+1}b + u_{2k+1}
⟺ s_{k} ≤ \frac{r_{k+1}b + u_{2k+1}}{2(s_{n−1}s_{n−2} … s_{k+1})_{b}}  (1.9)

となる。最後の除算の商を s_{k}、余りを r′_{k} とおけば、 r_{k} は

r_{k} = r_{k+1}b^{2} + u_{2k+1}b + u_{2k} − (2(s_{n−1}s_{n−2} … s_{k+1})_{b}b + s_{k})s_{k}
= r′_{k}b + u_{2k} − s_{k}^{2}  (1.10)

となる。 s_{k} が真値よりも大きければ r_{k} は負となるので、 r_{k} の値が正となるまで

r_{k} := r_{k} + 2s_{k} − 1, s_{k} := s_{k} − 1  (1.11)

を繰り返せば良い。

さて、除算のときと同じように繰り返しの回数を検討しよう。 k = n − 1 の場合は単純に平方根をとるだけであるから、 k ≤ n − 2 の場合を考える。

s′_{k} = ⌊\frac{r_{k+1}b + u_{2k+1}}{2(s_{n−1}s_{n−2} … s_{k+1})_{b}} ⌋

とおけば、繰り返しの回数は s′_{k} − s_{k} で与えられる。

s′_{k} ≤ \frac{r_{k+1}b + u_{2k+1}}{2(s_{n−1}s_{n−2} … s_{k+1})_{b}} ≤ \frac{r_{k+1}b^{2} + u_{2k+1}b + u_{2k}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b}} ,
s_{k} > \frac{r_{k+1}b^{2} + u_{2k+1}b + u_{2k}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b} + s_{k} + 1} − 1

であるから、

s′_{k} − s_{k} < \frac{r_{k+1}b^{2} + u_{2k+1}b + u_{2k}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b}} − \frac{r_{k+1}b^{2} + u_{2k+1}b + u_{2k}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b} + s_{k} + 1} + 1
= \frac{r_{k+1}b^{2} + u_{2k+1}b + u_{2k}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b} + s_{k} + 1} · \frac{s_{k} + 1}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b}} + 1
< \frac{b^{2}}{2(s_{n−1}s_{n−2} … s_{k+1}0)_{b}} + 1
≤ \frac{b^{k+3−n}}{2s_{n−1}} + 1  (1.12)

となる。この式から、繰り返しは k < n−2 の場合常に 1 回以下であり、また s_{n−1} ≥ b/2 ならば k = n−2 の場合も 1 回以下となることが分かる。よって s_{n−1} ≥ b/2 とするために、 u に適当な数 k^{2} を掛けて u_{2n−1} ≥ b/4 とすると、平方根が ks、余りが k^{2}r となるので、それぞれ k, k^{2} で割って補正すれば良い。実行時間は Θ(n^{2}) である。

1.5 基数変換

基数を 10 の冪以外にしている場合には、数を表示する際に基数変換が必要となることが多い。自然数 n を基数 b で表すには、 n を b で割った商と余りを求め、商を新たに n と置く、という操作を繰り

返せば良い。出てきた余りの列は、 n を b 進表記したときの下位の桁から順に並べたものになる。

次に整数部が 0 の正の実数 r の基数変換を考える。今度は r を b 倍して整数部と小数部に分け、小数部を新たに r と置く、という操作を繰り返せば良い。出てきた整数部の列は、 r の小数部を上位の桁から順に並べたものになる。

一般の実数は整数部と小数部に分け、それぞれを上に述べた方法で基数変換すれば良い。実行時間は Θ(n^{2}) である。

2 Karatsuba 系列のアルゴリズム

2.1 Karatsuba 乗算

2n 桁の数 u, v の乗算を行うことを考える。 u, v を基数 b の冪 x = b^{n} を用いて

u = u_{1}x + u_{0}, v = v_{1}x + v_{0}  (2.1)

のように上半分と下半分に分ければ、 x を基数とした乗算と考えることができる。このとき単純に w = uv を求めると、

w = u_{1}v_{1}x^{2} + (u_{0}v_{1} + u_{1}v_{0})x + u_{0}v_{0}

のように n × n^{*1}が 4 回必要である。しかし

u_{0}v_{1} + u_{1}v_{0} = u_{1}v_{1} + u_{0}v_{0} − (u_{1} − u_{0})(v_{1} − v_{0})  (2.2)

という関係を使えば n × n は 3 回で済むことが分かる。ここでさらに u_{1}, u_{0}, v_{1}, v_{0} の上位の桁を 0 で埋めて桁数を偶数にするか、あるいは始めから u, v の桁数を 2 の冪にしておけば、 3 回の n × n に同じアルゴリズムを再帰的に適用することができる。この乗算アルゴリズムを Karatsuba 乗算 (Karatsuba Multiplication) という。

このアルゴリズムで n × n に必要な時間を T (n) とする。 n × n には n/2 × n/2 が 3 回と n/2 ± n/2 が数回必要であるから、

T (n) = 3T (n/2) + cn/2
= 3^{lg n}T (1) + cn ((3/2)^{lg n} − 1)
= n^{lg 3}(T (1) + c) − cn
= Θ(n^{lg 3})  (lg 3 ≈ 1.585)  (2.3)

となり、 Θ(n^{2}) よりも少ない時間で乗算が実行できることが分かる。

2.2 Karatsuba 除算

4n 桁の数 u と 2n 桁の数 v (u < vb^{2n}) の除算を行い、 2n 桁の商 w と 2n 桁の余り r を求めることを考える。 u, v, w, r を基数 b の冪 x = b^{n} を用いて

u = u_{3}x^{3} + u_{2}x^{2} + u_{1}x + u_{0}, v = v_{1}x + v_{0}, w = w_{1}x + w_{0}  (2.4)

のように表せば、 x を基数とした除算と考えることができる。筆算のアルゴリズムから、 w_{1} を求めるためには除算 (u_{3}x + u_{2})/v_{1} を行わなければならないが、これはちょうど 4n ÷ 2n の半分の大きさであるから、 Karatsuba 乗算と同じように桁数を調節することで再帰的に求めることができる。後はこのようにして求めた商と余りを w_{1}, r′_{1} として、

r_{1} = r′_{1}x + u_{1} − v_{0}w_{1}  (2.5)

とおき、 r_{1} ≥ 0 となるまで

w_{1} := w_{1} − 1, r_{1} := r_{1} + v  (2.6)

*1 n 桁同士の乗算をこのように略記する。その他の演算も同様。

を繰り返せば w_{1} が求まる。この繰り返しを 2 回以下にするには、筆算のときと同様に v_{1} ≥ x/2 となるように、したがって b 進表示における v の最上位の桁が b/2 以上となるように調節すれば良い。 w_{0} についても同様にして求めることができる。 r_{1}/v_{1} の商 w_{0} と余り r′_{0} を再帰的に求め、

r_{0} = r′_{0}x + u_{0} − v_{0}w_{0}  (2.7)

とおき、 r_{0} ≥ 0 となるまで r_{0} := r_{0} + v, w_{0} := w_{0} − 1 を繰り返せば良い。ここで求まった r_{0} がそのまま余り r となる。この除算アルゴリズムを Karatsuba 除算 (Karatsuba Division) という。

このアルゴリズムで 2n ÷ n に必要な時間を T (n) とする。 2n ÷ n には n ÷ n/2 が 2 回と n/2 × n/2 が 2 回、 n + n 等が数回必要であるから、 n × n に必要な時間を M(n) として、

T (n) = 2T (n/2) + 2M(n/2) + cn
= nT (1) + ∑_{k=1}^{lg n} 2^{k}M(n/2^{k}) + cn lg n
= ∑_{k=1}^{lg n} 2^{k}M(n/2^{k}) + Θ(n lg n)  (2.8)

となる。具体的に M(n) = cn^{2} とすれば

T (n) = ∑_{k=1}^{lg n} cn^{2}/2^{k} + Θ(n lg n)
= cn^{2}(1 − 1/n) + Θ(n lg n)
= M(n) + Θ(n lg n)  (2.9)

となり、 M(n) = cn^{lg 3} とすれば

T (n) = ∑_{k=1}^{lg n} cn^{lg 3}(2/3)^{k} + Θ(n lg n)
= 2cn^{lg 3}(1 − n/n^{lg 3}) + Θ(n lg n)
= 2M(n) + Θ(n lg n)  (2.10)

となる。第 4 章で述べるアルゴリズムを使えば M(n) = cn lg n となるので、

T (n) = ∑_{k=1}^{lg n} cn lg(n/2^{k}) + Θ(n lg n)
= cn lg^{2} n − cn lg n(1 + lg n)/2 + Θ(n lg n)
= (M(n) lg n)/2 + Θ(n lg n)  (2.11)

となる。

2.3 Karatsuba 開平

4n 桁の数 u を開平して 2n 桁の平方根 s と余り r を求めることを考える。 u, s を基数 b の冪 x = b^{n} を用いて

u = u_{3}x^{3} + u_{2}x^{2} + u_{1}x + u_{0}
s = s_{1}x + s_{0}

のように表せば、 x を基数とした開平と考えることができる。筆算のアルゴリズムより、まず u_{3}x + u_{2} の平方根と余りを求めなければならないが、これはちょうど u の半分のサイズであるから再帰的に求めることができる。このようにして求めた平方根を s_{1}、余りを r_{1} とおく。次に s_{0} を求めるために、 (r_{1}x + u_{1})/2s_{1} の商を s_{0}、余りを r′_{0}、そして r_{0} = r′_{0}x + u_{0} − s_{0}^{2} とおく。 r_{0} ≥ 0 となるまで

r_{0} := r_{0} + 2s_{0} − 1, s_{0} := s_{0} − 1  (2.12)

を繰り返せば s_{0} が求まる。特に u_{3} ≥ b/4 の場合には、 1 回条件分岐をするだけとなる。この開平アルゴリズムを Karatsuba 開平 (Karatsuba Square Root) という。

このアルゴリズムで 2n 桁の数を開平するのに必要な時間を T (n) とする。 4n 桁の開平には sqrt(2n) が 1 回、 2n ÷ n が 1 回、 n × n が 1 回、 2n − 2n が 1 回、 n + n 等が数回必要である。 2n ÷ n に必要な時間を D(n)、 n × n に必要な時間を M(n) とすれば、

T (n) = T (n/2) + D(n/2) + M(n/2) + cn/2
= T (1) + ∑_{k=1}^{lg n}(D(n/2^{k}) + M(n/2^{k}) + cn/2^{k})
= ∑_{k=1}^{lg n}(D(n/2^{k}) + M(n/2^{k})) + Θ(n)  (2.13)

となる。ここで Karatsuba 乗除算を使えば

T (n) = ∑_{k=1}^{lg n}(3M(n/2^{k}) + Θ(n/2^{k} lg(n/2^{k}))) + Θ(n)
= (9M(n/2))/2 + Θ(n lg n)
= 3M(n)/2 + Θ(n lg n)  (2.14)

となる。

2.4 Karatsuba 基数変換

例えば 3141592 という 10 進表現は、 3, 141, 592 のように 3 桁ごとに区切れば 10^{3} = 1000 進表現としてみることができる。このように基数を b から b^{n} に変換する、あるいはその逆は単に (b 進表示における) 桁の区切りを変えるだけで良く、隣の桁に対する繰上げや繰り下げが起こることはない。このことを利用したのがここで述べるアルゴリズムである。

2n_{0} 桁の b_{0} 進自然数 u を基数変換し、 b_{1} 進表現を得たいとする。まず b_{0}^{n_{0}} ≤ b_{1}^{n_{1}} となるような最小の n_{1} を求め、 x = b_{1}^{n_{1}} とする。このとき、 u の b_{1} 進表示は 2n_{1} 桁以下である。次に u の x 進表示を得るために u/x を行い、その商を q、余りを r とする。そうすると u は x 進表示では q, r の 2 桁で表されることになるが、先に述べたことから q, r はそのまま b_{1} 進表示における上位 n_{1} 桁と下位 n_{1} 桁になっている。後は q, r に対して同じアルゴリズムを再帰的に適用すれば u の b_{1} 進表示が得られる。この基数変換アルゴリズムを Karatsuba 基数変換 (Karatsuba Radix Conversion) という。実際に必要とされる基数は限られているであろうから、予め b_{1}^{2^{n}} 等を計算しておくとよい。

このアルゴリズムで b_{0} 進 n 桁の数を b_{1} 進表現に変換するのに必要な時間を T (n) とする。 b_{1} の冪を予め求めておけば単に除算を繰り返すだけであるから、

T (n) = 2T (n/2) + D(n/2)

= ∑_{k=1}^{lg n} D(n/2^{k})  (2.15)

となる。具体的に Karatsuba 乗算を用いて D(n) = 2M(n) + Θ(n lg n) とすれば、

T (n) = ∑_{k=1}^{lg n}(2M(n/2^{k}) + Θ(n/2^{k} lg(n/2^{k})))
= M(n) + Θ(n lg n)  (2.16)

となる。

3 Karatsuba 乗算の一般化

3.1 数と多項式

b 進 n 桁の数 u = (u_{n−1}u_{n−2} … u_{0})_{b} の値は

u_{n−1}b^{n−1} + u_{n−2}b^{n−2} + ⋯ + u_{0}

であるが、これは多項式

u(x) = u_{n−1}x^{n−1} + u_{n−2}x^{u−2} + ⋯ + u_{0}  (3.1)

において、 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 を多項式で行うと、

u(x) = x^{2} + 2x + 8
v(x) = x^{2} + 8x + 6
w(x) = u(x) + v(x) = 2x^{2} + 10x + 14

となり、 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_{2}x^{2} + u_{1}x + u_{0}, v(x) = v_{2}x^{2} + v_{1}x + v_{0}

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

w(x) = u_{2}v_{2}x^{4} + (u_{1}v_{2} + u_{2}v_{1})x^{3} + (u_{0}v_{2} + u_{1}v_{1} + u_{2}v_{0})x^{2} + (u_{0}v_{1} + u_{1}v_{0})x + u_{0}v_{0}

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

u(x) = ∑_{n=0}^{∞} u_{n}x^{n}, v(x) = ∑_{n=0}^{∞} v_{n}x^{n}  (3.2)

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

w(x) = ∑_{n=0}^{∞} w_{n}x^{n}  (3.3)

のように表されているとき、 w_{n} は

w_{n} = ∑_{k=0}^{n} u_{k}v_{n−k}  (3.4)

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

u_{0}v_{1} + u_{1}v_{0} = u_{1}v_{1} + u_{0}v_{0} − (u_{1} − u_{0})(v_{1} − v_{0})

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

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

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

u(x) = u_{1}x + u_{0}
v(x) = v_{1}x + v_{0}
w(x) = w_{2}x^{2} + w_{1}x + w_{0}

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

w(x_{i}) = u(x_{i})v(x_{i})  (i = 0, 1, 2)  (3.5)

を解けばよいことになる。具体的に x_{0} = 0, x_{1} = 1, x_{2} = −1 とすれば

w_{0} = u_{0}v_{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})

となって

w_{0} = u_{0}v_{0}
w_{1} = u_{1}v_{1} + u_{0}v_{0} − (u_{1} − u_{0})(v_{1} − v_{0})  (3.6)
w_{2} = u_{1}v_{1}

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

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

u(x) = u_{2}x^{2} + u_{1}x + u_{0}
v(x) = v_{2}x^{2} + v_{1}x + v_{0}  (3.7)
w(x) = w_{4}x^{4} + w_{3}x^{3} + w_{2}x^{2} + w_{1}x + w_{0}

を考え、 w(x_{i}) = u(x_{i})v(x_{i}) (i = 0, 1, …, 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
= Θ(n^{log_{3} 5})  (log_{3} 5 ≈ 1.465)  (3.8)

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

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

3.3 多項式補間

問題は (x_{i}, y_{i}) (i = 0, …, n − 1) という点を通る多項式 w(x) = ∑_{k=0}^{n−1} w_{k}x^{k} を決定することである。これを多項式補間 (Polynomial Interpolation) という。ここでは特に x_{i} = i という制限を加えて補間を行いやすくする。

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

Δx\underline{n} = (x + 1)x ⋯ (x − n + 1) − x(x − 1) ⋯ (x − n + 1)
= nx(x − 1) ⋯ (x − n + 2)

となり、 Δx\underline{1} = 1 なので、 n ≥ 1 に対して

Δx\underline{n} = nx\underline{n−1}  (3.9)

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

x^{n} = ∑_{k=0}^{n} \left\{{n \atop k}\right\}x\underline{k}  (3.10)

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

w(x) = ∑_{k=0}^{n−1} c_{k}x\underline{k}  (3.11)

と書き換えたとき、 c_{k} は全て負でない。ここで

Δ^{0}w(x) = c_{0} + c_{1}x\underline{1} + c_{2}x\underline{2} + c_{3}x\underline{3} ⋯
Δ^{1}w(x) = c_{1} + 2c_{2}x\underline{1} + 3c_{3}x\underline{2} + ⋯
Δ^{2}w(x) = 2c_{2} + 6c_{3}x\underline{1} + ⋯
Δ^{3}w(x) = 6c_{3} + ⋯
⋮

であるから、 c_{k} は

c_{k} = \frac{Δ^{k}w(0)}{k!}  (3.12)

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

3.4 Toom-Cook 法

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

\begin{aligned}u(0) & = & 5, & u(1) & = & 13, & u(2) & = & 29, & u(3) & = & 53, & u(4) & = & 85 \\ v(0) & = & 7, & v(1) & = & 31, & v(2) & = & 81, & v(3) & = & 157, & v(4) & = & 259 \\ w(0) & = & 35, & w(1) & = & 403, & w(2) & = & 2349, & w(3) & = & 8321, & w(4) & = & 22015 \end{aligned}

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

0 : 35  403  2349  8321  22015

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

0 : 35  403  2349  8321  22015 1 :  368  1946  5972  13694

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

0 : 35  403  2349  8321  22015
1 :  368  1946  5972  13694
2 :  789  2013  3861

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

0 : 35  403  2349  8321  22015
1 :  368  1946  5972  13694
2 :  789  2013  3861
3 :  408  616
4 :  52

このとき、 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} の係数は次のように計算できることが分かる。

52 408
−  3 · 52
52 252 789
−  2( 52 252)
52 148 285 368
−  1( 52 148 285)
52  96 137 83 35

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

35 = 2 · 16 + 3
83 + 2 = 5 · 16 + 5
137 + 5 = 8 · 16 + 14
96 + 8 = 6 · 16 + 8
52 + 6 = 3 · 16 + 10

となるから、

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

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

4 高速 Fourier 変換を利用した乗算法

4.1 多項式の乗算の行列による表現

先の章で述べたことを行列を用いて整理してみよう。多項式

u(x) = u_{0} + u_{1}x + ⋯ + u_{n−1}x^{n−1}  (4.1)

に n 個の値 x_{0}, x_{1}, …, x_{n−1} を代入すれば、 n 個の値 y_{0}, y_{1}, …, y_{n−1} が得られる。行列を用いて表現すれば

( y_{0} ; y_{1} ; ⋮ ; y_{n−1} ) = ( \begin{matrix}1 & x_{0} & … & x_{0}^{n−1} \\ 1 & x_{1} & … & x_{1}^{n−1} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{n−1} & … & x_{n−1}^{n−1} \end{matrix} ) ( u_{0} ; u_{1} ; ⋮ ; u_{n−1} )  (4.2)

となる。 x_{k}, y_{k} が分かっているならば、上の式に (もし存在するならば) 逆行列を掛けることで u_{k} が求まる (補間)。逆行列が存在するように x_{k} を選びさえすれば、 u_{k} の代わりに y_{k} を使って多項式を表すことができるのである。問題の行列は Vandermonde 行列 (Vandermonde Matrix) と呼ばれているもので、

V (x_{0}, x_{1}, …, x_{n−1}) = ( \begin{matrix}1 & x_{0} & … & x_{0}^{n−1} \\ 1 & x_{1} & … & x_{1}^{n−1} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & x_{n−1} & … & x_{n−1}^{n−1} \end{matrix} )  (4.3)

と表され、行列式は

detV (x_{0}, x_{1}, …, x_{n−1}) = ∏_{i>j}(x_{i} − x_{j})  (4.4)

で与えられるので、 x_{k} が全て異なれば逆行列が存在する。

多項式の乗算 w(x) = u(x)v(x) を行列を用いて大雑把に述べると、 u(x), v(x), w(x) の係数を並べた列ベクトルをそれぞれ u, v, w として、

1. V (x_{0}, x_{1}, …, x_{n−1}) を u, v に掛け、 u′, v′ を得る。

2. u′, v′ を要素ごとに掛け、 w′ を得る。

3. V (x_{0}, x_{1}, …, x_{n−1})^{−1} を w′ に掛け、 w を得る。

となる。大雑把というのは必要な次元を考慮していないからである。上で述べたままだと、 w(x) が n − 1 次、つまり u(x), v(x) の上位の係数合わせて n − 1 個が 0 でないと上手くいかない。これは w(x) = u(x)v(x) の次数が u(x), v(x) の次数よりも高くなるからで、 Toom-Cook 法で連立方程式を立てる際に、 u(x), v(x) が 2 次なら w(x) = u(x)v(x) は 4 次となり、従って x に 5 個の異なる値を代入しなければならなかったのと同じである。

Toom-Cook 法は 2. の要素ごとの乗算を再帰的に行うアルゴリズムである。従って 1. や 3. における V (x_{0}, x_{1}, …, x_{n−1}) の乗算をいかに速く行うかが問題となるが、 u, v の分割数を大きくするとこの部分が複雑になり、位数も余り下がらなくなるので、大した効果は期待できなくなる。ここで発想を変えて、 1. や 3. で再帰を行うのが次に述べる高速 Fourier 変換である。

4.2 高速 Fourier 変換

Vandermonde 行列の要素を

x_{k} = ω_{k} = e^{2kπi/n}  (4.5)

としたときに得られる u′ を u(x) の離散 Fourier 変換 (Discrete Fourier Transform, DFT) と呼び、 DFT_{n}(u) で表す。特に n を 2 の冪に限った場合、高速 Fourier 変換 (Fast Fourier Transform, FFT) というアルゴリズムによって Θ(n lg n) の時間で DFT_{n}(u) が求まるのである。

FFT の基本的な考え方は簡単である。

u(x) = u_{0} + u_{1}x + ⋯ + u_{n−1}x^{n−1}  (4.6)

とおくと、 DFT_{n}(u) を得るためには u(ω_{k}) を評価しなければならないが、もし n = 2m ならば u_{k} を k の偶奇で分けた 2 つの m − 1 次多項式

u_{0}(x) = u_{0} + u_{2}x^{2} + ⋯ + u_{2m−2}u^{m−1}  (4.7)
u_{1}(x) = u_{1} + u_{3}x^{2} + ⋯ + u_{2m−1}u^{m−1}  (4.8)

を用いることによって

u(ω_{k}) = u_{0}(ω_{k}^{2}) + ω_{k}u_{1}(ω_{k}^{2})
= u_{0}(ω_{2k}) + ω_{k}u_{1}(ω_{2k})  (4.9)

のように計算することができる。

ここで k = 0, 1, …, 2m − 1 について u(ω_{k}) の値を求めることを考えると、 ω_{k} は 2m 個の値をとるのに対して、 ω_{2k} は m 個の値しかとらず、しかもそれは複素 m 乗根であることが分かる。つまり DFT_{n}(u) を得る問題は DFT_{m}(u_{0}), DFT_{m}(u_{1}) を得る問題に帰着されるのである。具体的に n = 4 の場合を考えると、

u(ω_{0}) = u_{0}(ω_{0}) + ω_{0}u_{1}(ω_{0})
u(ω_{1}) = u_{0}(ω_{2}) + ω_{1}u_{1}(ω_{2})
u(ω_{2}) = u_{0}(ω_{0}) − ω_{0}u_{1}(ω_{0})
u(ω_{3}) = u_{0}(ω_{2}) − ω_{1}u_{1}(ω_{2})

のように計算でき、計算量がほぼ半分に減るのである。当然 u_{0}(x), u_{1}(x) にも同じ方法を適用できれば計算量を減らすことができるので、 n が 2 の冪であれば最も効率が良いが、これは 0 を係数として加えるだけで実現できる。結局このアルゴリズムで実際に乗算が行われる部分は ω_{m}u_{1}(ω_{2m}) のみとなり、 DFT_{n}(u) を得るのに必要な時間 T (n) は

T (n) = 2T (n/2) + cn
= nT (1) + cn lg n
= Θ(n lg n)  (4.10)

となる。このアルゴリズムを再帰的 FFT(Recursive FFT) という。後は逆変換 (Inverse DFT, DFT^{−1}) が高速に実行できれば良いが、

V (ω_{0}, ω_{1}, …, ω_{n−1})^{−1} = \frac{1}{n} ( \begin{matrix}1 & 1 & … & 1 \\ 1 & ω_{1}^{−1} & … & ω_{1}^{−(n−1)} \\ ⋮ & ⋮ & ... & ⋮ \\ 1 & ω_{n−1}^{−1} & … & ω_{n−1}^{−(n−1)} \end{matrix} )  (4.11)

であるあから、 FFT を少し変更するだけで逆変換が行え、従って必要な時間も Θ(n lg n) となる。

以上のことを考えると、多項式の乗算は次のようになる。下準備として、 w(x) = u(x)v(x) の次数より小さくない最小の 2 の冪を n とし、 u(x), v(x) の次数が n となるように 0 を高次の係数として加えておく。そうすると、

1. FFT を用いて u′ = DFT_{n}(u), v′ = DFT_{n}(v) を求める。必要な時間は Θ(n lg n)。

2. u′, v′ の要素同士を掛けて w′ を求める。必要な時間は Θ(n)。

3. FFT を用いて w = DFT_{n}^{−1}(w′) を求める。必要な時間は Θ(n lg n)。

という計算を行うことによって乗算が Θ(n lg n) の時間で実行できることになる。

4.3 Cooley-Tukey のアルゴリズム

再帰を行わず、繰り返しのみで FFT が実行できれば速度が向上するはずである。そこで u_{k} がどのように振り分けられるのかを考える。例えば n = 8 の場合、添字だけを書けば

(0, 1, 2, 3, 4, 5, 6, 7)

(0, 2, 4, 6), (1, 3, 5, 7)
(0, 4), (2, 6), (1, 5), (3, 7)
(0), (4), (2), (6), (1), (5), (3), (7)

となる。少し考えると、 k の 2 進表記の上位と下位を入れ替えた値が u_{k} の位置となることが分かる。このようにして u_{k} を予め並べ替えておけば、繰り返しだけで簡単に FFT が実行できる。この方法を Cooley-Tukey のアルゴリズム (Cooley-Tukey Algorithm) という。

5 用語について

この文書で用いている用語についての注意などを述べる。

多倍長演算を意味する語は日本語では恐らくこれだけであるが、英語では multiple precision の他に multiprecision、 arbitrary precision、 multidigit といった形容詞や、 arithmetic の代わりに computation も使うようである。多倍長数という日本語は余り使われていないようであり、英語でも multiple precision number よりも単に bignum とするようである。筆算も日本語ではこれしかないと思われるが、英語では他に schoolboy’s method などが使われる。

Karatsuba 乗算、 Karatsuba 除算、 Karatsuba 開平、 Karatsuba 基数変換は日本語の文献には見られないが、英語との対応や分かりやすさを考慮して使用した。ただし Karatsuba radix conversion は、名前が見当たらなかったので筆者が勝手に付けたものである。 Karatsuba 乗算は、 Karatsuba-Ofman 法 (Karatsuba-Ofman method)、ディジタル法 (digital method)、 divide and conquer method、 fractal multiplication などの呼び名がある。ただし、このアルゴリズムは Karatsuba にのみクレジットされているようであり、またここで紹介したものはオリジナルとは少し違うようであるが、論文が手に入らないので詳しいことは分からない。

降冪は、日本語では一般に階乗冪、下降階乗冪などと呼ばれるが、冗長に感じるので筆者は勝手に降冪をこの意味にも使っている。英語では falling power、 falling factorial power の両方が使われるようである。また Pochhammer の記号 (Pochhammer’s symbol) (x)_{n} を降冪として使うこともあるが、種々の「冪ではないが、冪のようなもの」の定義に使われているので紛らわしい。

Cooley-Tukey のアルゴリズムは、繰り返し FFT(iterative FFT)、 Cooley-Tukey 法 (Cooley-Tukey method)、 Cooley-Tukey FFT などの呼び名がある。

参考文献

[1] Donald E. Knuth, 中川圭介 (訳), 『準数値算法/算術演算』 , サイエンス社, 1986.

[2] Torbjörn Granlund, ”GNU MP: The GNU Multiple Precision Arithmetic Library” Edition 4.1, 2002. http://swox.com/gmp/

[3] Thomas H. Cormen, Charles E. Leiserson and Ronald L. Rivest, ”Introduction to Algorithms”, MIT Press, 1990.

[4] Christoph Burnikel and Joachim Ziegler, ”Fast Recursive Division”, Research Report MPI-I-98-1-022, Max-Planck-Institute für Informatik, 1998. http://data.mpi-sb.mpg.de/internet/reports.nsf/NumberView/1998-1-022

[5] Paul Zimmermann, ”Karatsuba Square Root”, Research Report RR-3805, INRIA, 1999. http://www.inria.fr/rrrt/rr-3805.html

[6] D. J. Bernstein, ”Multidigit multiplication for mathematicians”, Advances in Applied Mathematics, to appear, 2001. http://cr.yp.to/papers/m3.ps