あざらしとペンギンの問題

主に漫画、数値計算、幾何計算、TCS、一鰭旅、水族館、鰭脚類のことを書きます。

Toom-Cook法:掛け算の数が減るカラクリ

今回は、Karatsuba法をダシにして、「高みの見物で構造が見えてくる」というお話をします。

最初の高速乗算アルゴリズムであるKaratsuba法は、それだけ見ると天下り的に上手いことやって式を作ったように思えます。あるいは実際そうだったのかもしれません。発明者のKaratsuba氏*1は主に数論の方で多くの研究をされてきたようなので、あるいは数論的な背景があったのかもしれません。この点に関しては、本人が既に亡くなっているので、真相は神のみぞ知るといったところです。

しかし、Karatsuba法の一般化であるToom-Cook法を見ると、掛け算の回数が減らせる仕組みが明らかになります。逆に言えば、その仕組みが判ったからこそ、一般化できたとも言えます。我々は一般化された形を見ることによって、特殊ケースの裏にある構造を容易に意識(ここが一番重要!)し、理解することができるのです。

イデア

前に言葉だけで触れたように、Toom-Cook法は数をもっと細かく分割して掛け算をする方法です。例えば、2つの n 桁の数 a, b に対して、両方を d 分割して行うものはToom-d と呼ばれます。この意味でKaratsuba法はToom-2になっています。より一般には、ab を異なる数に分割をすることも可能です。今回はKaratsuba法の考察が目的なので、より単純な一般化であるToom-d だけについて見ることにします。

基数 \ell は任意として入力 a, bn 桁の非負整数とします。簡単のため、nd で常に割り切れるとします。目的はもちろん、積 a \times b を計算することです。

分割

Toom-Cook法では、まず入力を d 分割して n/d 桁ずつにします。右(最下位)から左(最上位)に向かって aa_0, a_1, \dots, a_{d-1} に、bb_0, b_1, \dots, b_{d-1} に分割されるとすると、B = \ell^{n/d} として

a = a_0 + a_1 B + a_2 B^2 + \dots + a_{d-1} B^{d-1}
b = b_0 + b_1 B + b_2 B^2 + \dots + b_{d-1} B^{d-1}

となります。

そして、掛け算の結果は

R = a \times b = R_0 + R_1 B + R_2 B^2 + \dots + R_{2d-2} B^{2d-2}

と書けます。後はこの 2d-1 個の R_0, R_1, \dots, R_{2d-2} をどうやって計算するかが問題となります。

Toom-Cook法では、結果の分割の数と等しい回数の掛け算で R_0, R_1, \dots, R_{2d-2} の全てを計算します。つまり、この場合は 2d-1 回の掛け算でできるというわけです。

多項式化と係数の決定

Toom-Cook法で掛け算の回数を減らす鍵は、先の入力および結果の分割を多項式と見なすことにあります。すなわち、上の a, b, R の式において Bx に置き換え、

a(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_{d-1} x^{d-1}
b(x) = b_0 + b_1 x + b_2 x^2 + \dots + b_{d-1} x^{d-1}
R(x) = a(x) \times b(x) = R_0 + R_1 x + R_2 x^2 + \dots + R_{2d-2} x^{2d-2}

とするのです。最初の時点で R(x) の係数は未知です。そこで、a, b の分割を用いて R(x) の係数を決めてやるのです。最終的な結果は R = R(B) として得られます。

ここで、係数 R_k は普通の多項式の掛け算でも計算することができ、すべて整数値です。ただし、そのように計算すると、すべての係数を得るのに d^2 回の掛け算が必要になってしまいます。これを減らすための一般的な方法を見ることが今回の話の主題です。

Toom-Cook法の基本的なアイデアは、R(x) の係数 R_0, R_1, \dots, R_{2d-2} を決めるために、それらと同じ 2d-1 個の R(x) の値を求めて、連立方程式を解いてやるというものです。これ自体は高校数学でよく見かける問題です。例えば、2次と判っている多項式 R(x) について、

R(0) = 3
R(1) = 5
R(2) = 3

であるとき、R(x) を決定せよという問題が与えられたとします。これを解くには、R(x) = R_0 + R_1 x + R_2 x^2 とおいて、左辺に各 x の値を代入した式、右辺にその値を入れた

R(0) = R_0 = 3
R(1) = R_0 + R_1 + R_2 = 5
R(2) = R_0 + 2R_1 + 4R_2 = 3

連立方程式として立てます。これを解けば、R(x) の係数が

R_0 = 3, R_1 = 4, R_2 = -2

と求まりますね。Toom-Cook法ではこのような計算を行います。

一般には、R(x) の次数を t-1 として、食わせる値を x_1, x_2, \dots, x_t とすると、線形方程式系の一般論より、次の行列(以下、「変数行列」と呼ぶことにします。)

f:id:azapen6:20140412030934p:plain

が可逆ならば、解はいつでもちょうど一つあります。また、適当な方法でそれを計算することができます。加えて、変数行列はいわゆるヴァンデルモンド行列と呼ばれるものになっていて、x_1, x_2, \dots, x_t がすべて異なれば可逆です。よって、t 個の異なる点で R(x) の値を求めることができたなら、それらから係数 R_0, R_1, \dots, R_{t-1} を求めることができます。

元の問題に戻りましょう。上の考察より、R(x) の値を 2d-1 個の点 x_1, x_2, \dots, x_{2d-1} で計算できれば、求めたい係数 R_0, R_1, \dots, R_{2d-2} が決まります。

以上を図に表すと次のようになります。

f:id:azapen6:20160412164211p:plain

同じ四角にまとめたものは同じ種類のもので、左から順に、数、多項式多項式の評価値となっています。上二つ(掛ける数の処理)と下(結果の算出)の推移が逆変換の関係になっていることが見て取れるでしょう。

掛け算の数が減るカラクリ

Toom-Cook法では、主要(掛ける2数の桁数が多い、上の図では右の四角の掛け算)な掛け算を指定した点 x = x_k に対する

R(x_k) = a(x_k) \times b(x_k)

を計算するときにだけ行います。これを x_1, x_2, \dots, x_{2d-1} のすべてに対して行うため、ちょうど 2d-1 回の掛け算が行われるのです。仮に R(x)多項式として先に計算しようとすると、最初に言った通り d^2 回の掛け算が必要になります。しかし、R(x) の係数を決めるという発想のもとで、必要十分な数の点で上を評価してやれば、掛け算を 2d-1 回に減らせるというわけです。これこそがToom-Cook法並びにKaratsuba法で掛け算の回数が減るカラクリです。

ここで、主要でない掛け算というのは、一方の数が入力に依存しない(+小さい)定数であるものを言います。例えば x=x_k に対して a_i x_k^i を計算するときの掛け算のことです。この場合では、x_k^i はプログラムに予め書き込まれる数です。このような掛け算については、筆算によって高々 O(n) の計算量(足し算、引き算と同じオーダー)で実行できるため、主要でない掛け算として足し算、引き算と同じ扱いをします。

後は、主要な掛け算を再帰的に行うことで、全体の計算量も減らせるというわけです。

アルゴリズム

以上の考察を元にアルゴリズムを書き下すと、次のようになります。x_1, x_2, \dots, x_{2d-1} はパラメータとして最初に決めます。

入力:\ell 進法 n 桁で書かれた非負整数 a, b
出力:\ell 進法で書かれた a \times b

  1. 入力を分割して a_0, a_1, \dots, a_{d-1}, b_0, b_1, \dots, b_{d-1} を得る。
  2. k=1, 2, \dots, 2d-1 について3-4を行う。
  3. a(x_k), b(x_k) を計算する。
  4. R(x_k) = a(x_k) \times b(x_k)再帰的に計算する。
  5. 連立方程式 R_0 + R_1 x_k + \dots + R_{2d-2} x_k^{2d-2} = R(x_k) \quad (k = 1, 2, \dots, 2d-1) を解いて R(x) の係数を確定する。
  6. B=\ell^{n/d} として R(B) を計算し、出力する。

計算量

n 桁の入力に対する計算量 T(n) を評価しましょう。再帰の各ステップで行われる主な計算を抜き出すと、

  • a(x_k), b(x_k) の計算×2d-1
  • 4の掛け算×2d-1
  • 5の連立方程式を解く計算
  • R(B) の計算

まず、4の主要な掛け算について見ます。これは再帰によって行うので、入力となる a(x_k), b(x_k) の桁数の最大値を L と置くと、一つの計算あたりの計算量が T(L) で、それを 2d-1 回行うので、合計で (2d-1) T(L) の寄与になります。L は大雑把には n/d と見なせる程度のものですが、実際には x_k の値によって変わります。これが少し面倒なところです。x_k には 2d-1 個の異なるものを持ってくればいいので、仮に x_k=0, \pm 1, \pm 2, \dots, \pm (d-1) としましょう。このとき、|a(x_k)| \leq d (d-1) B となります。ここで、|a_k| < B を用いました。b(x_k) についても同様です。よって、L \leq \lceil \log_{\ell} (d (d-1) B) \rceil \leq \lceil n/d + 2 \log_{\ell} d \rceil となります。これより、4に要する計算量は高々 (2d-1) T(\lceil n/d + 2 \log_{\ell} d \rceil) です。

以下、残りを全部合わせても O(n) の計算量で実行できることを示します。

a(x_k) の計算については、a(x_k) = a_0 + a_1 x_k + \dots + a_{d-1} x_k^{d-1} をそのまま行うだけですが、各演算は足し算(引き算)、主要でない掛け算のみで、それぞれが O(n) の計算量で実行できるので、一つの a(x_k) の計算は O(d) \times O(n) = O(dn) の計算量で実行できます。b(x_k) についても同様です。それらを 2d-1 回行うので、全体の計算量は O(d) \times O(dn) = O(d^2 n) となり、d が定数なので O(d^2 n) = O(n) となります。

R(B) については、筆算の最後における足し算と同じく、n 桁の数をシフトしながら順次加えていくだけなので、計算量はやはり O(n) です。

5の連立方程式を解いて係数を求める計算についても、d を定数とすれば O(n) で実行可能であることが言えます。例えば、変数行列の逆行列を事前(プログラムを組むとき)に計算しておいて掛けるという実装にすると、実行の際にはそれぞれ d^2 回の足し算(引き算)と、主要でない掛け算だけで行えます。ただし、逆行列には分数、すなわち割り算が含まれる可能性があるので注意を要します。割り算についてはまだ扱っていませんが、割る数が定数ならば、筆算によって O(n) の計算量で実行できます。よって、係数を求める部分全体の計算量は O(d^2) \times O(n) = O(d^2 n) = O(n) となります。

以上から T(n) の漸化式は

T(n) = (2d-1) T(\lceil n/d + 2 \log_{\ell} d \rceil) + O(n) \quad (n > 2 \log_{\ell} d  + 1)
T(n) \leq 5 n^2 \quad (n \leq 2 \log_{\ell} d  + 1)

と書けます。ただし、初項だけでは値が定まらないので、最初の方は筆算の計算量で抑えておくことにします。つまり、そのようにプログラムを組んだときの計算量を考えているということです。実装上も桁数が少ないところでは筆算でやる方が速いです。

厳密な評価は見送るとして、上の漸化式を大雑把に評価すると、次の漸化式で定められる数列 \tilde{T}(n)T(n) と同程度のオーダーになると考えられます。

\tilde{T}(n) = (2d-1) \tilde{T}(n/d) + O(n)
\tilde{T}(1) = 1

この数列は所謂「マスター定理」より

\tilde{T}(n) = O(n^{\log_d (2d-1)})

というオーダーを持ちます。これより

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

であると予想されます。

具体例

以下では、具体例によってKaratsuba法との対応を確認します。

Toom-Cook法を実装するには、まず、x_k の値を決めてやることが必要です。それには、計算が簡単になるものを選ぶのがよいでしょう。便宜上、x_k = \top を代入すると最高次の項の係数が得られるとします。(概念的なもので計算は何もしません。)逆に、x_k = 0 を代入すると最低次の項が得られますね。これら2つを最初に置いて、残りを決めるのが通例です。x_k = \top を置いた場合、変数行列の一つの行は 0 0 \dots 0 1 になりますが、これによって可逆性が崩れることはありません。また、R(\top) = a(\top) \times b(\top) がそのまま成り立ちます。

x_k を決めた後は、R(x) の係数を求めるための逆行列を事前に計算しておくと便利です。ただし、実装上は必ずしも逆行列をそのまま掛ける必要はありません。

Toom-2

d = 2 の場合は 2d - 1 = 3 なので、x_1 = 0, x_3 = \top 以外にもう一つ用意します。ここでは x_2 = -1 としましょう。

すると、変数行列は

f:id:azapen6:20140412031006p:plain

となります。これの逆行列を計算すると、

f:id:azapen6:20140412031016p:plain

となります。(元の X と一致したのはただの偶然です。)これに

R(0) = a(0) \times b(0) = a_0 \times b_0
R(1) = a(1) \times b(1) = (a_0 - a_1) \times (b_0 - b_1)
R(\top) = a(\top) \times b(\top) = a_1 \times b_1

を縦にならべたベクトルに掛けて R(x) の係数を出すと、

R_0 = R(0) = a_0 \times b_0
R_1 = R(0) + R(\top) - R(1) = R_0 + R_2 + (a_1 - a_0) \times (b_0 - b_1)
R_2 = R(\top) = a_1 \times b_1

となります。

これらの式をよく見ると、Karatsuba法そのものであることが解ります。結局、Karatsuba法はToom-2において x_1 = 0, x_2 = -1, x_3 = \top と置いた場合だったのです。

Toom-3

d = 3 の場合は 2d - 1 = 5 なので、x_k を5個用意します。ここでは x_1 = 0, x_2 = -1, x_3 = 1, x_4 = 2, x_5 = \top と置きましょう。このとき、変数行列は

f:id:azapen6:20140412031230p:plain

となります。逆行列を計算すると、

f:id:azapen6:20140412031049p:plain

となります。これを

R(0) = a(0) \times b(0) = a_0 \times b_0
R(-1) = a(-1) \times b(-1) = (a_0 - a_1 + a_2) \times (b_0 - b_1 + b_2)
R(1) = a(1) \times b(1) =  (a_0 + a_1 + a_2) \times (b_0 + b_1 + b_2)
R(2) = a(2) \times b(2) =  (a_0 + 2 a_1 + 4 a_2) \times (b_0 + 2 b_1 + 4 b_2)
R(\top) = a(\top) \times b(\top) = a_2 \times b_2

を縦にならべたベクトルに掛けて R(x) の係数を出すと、便宜上一部を6倍で計算して

R_0 = R(0)
6 R_1 = -3 R(0) - 2 R(-1) + 6 R(1) - R(2) + 12  R(\top)
6 R_2 = -6 R(0) +3 R(-1) + 3 R(1) - 6  R(\top)
6 R_3 = -3 R(0) - R(-1) - 3 R(1) + R(2) - 12  R(\top)
R_4 = R(\top)

となります。各 R_k は整数なので、R_1, R_2, R_3 の値は上で求めたものを6で割れば厳密に出せます。

まとめ

今回の話では、Karatsuba法で掛け算の数が減らせる仕組みを見るために、それを一般化したToom-Cook法について見ました。そして、Karatsuba法が特別な場合であることを具体例として示しました。

このように、一般化によって背後に隠れていた構造を明らかにすることは、理論における一般化のモチベーションの一つです。この意味で、Toom-Cook法はKaratsuba法の意味のある一般化になっています。逆に、無意味な一般化というのは、ただパラメータ的に拡張を行っただけで、見逃していた構造が浮かび上がるでもなく、他に応用があるわけでもなく、元から特に情報量が増えないというものです。当然ながら、無意味な一般化は価値を認められません。

あなたがもし何かを一般化できるかもしれないと思ったならば、それが意味のある一般化になり得るかどうか、最初に考えて、また少し進んだら振り返って、今一度意義を確認して進めてください。

*1:露:Анато́лий Алексе́евич Карацу́ба、1937-2008