円周率.jp > Toom-Cook 法

Toom-Cook 法

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

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

A Ad-1Rd-1Ad-2Rd-2 + … + A0 A(x) Ad-1xd-1Ad-2xd-2 + … + A0
B Bd-1Rd-1Bd-2Rd-2 + … + B0 B(x) Bd-1xd-1Bd-2xd-2 + … + B0
AB C C2d-2R2d-2C2d-3R2d-3 + … + C0 A(x)B(x) C2d-2x2d-2C2d-3x2d-3 + … + C0

手順

先に条件を述べておくと, 手順 3 での逆演算を容易にするため, 手順 1 で代入する値は絶対値が小さな整数を選ぶ. また,手順 3 において多倍長同士の乗算・ 除算は行われない方が実装が容易であり, (おそらく)演算も速いため, そのような範囲で分割数 d が制限される.

ここでは d=3 の例を提示しながらその計算手順を示す.

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

A(-2)= 4A2 − 2A1A0
A(-1)A2A1A0
A(0)A0
A(1)A2A1A0
A(2)= 4A2 + 2A1A0
B(-2)= 4B2 − 2B1B0
B(-1)B2B1B0
B(0)B0
B(1)B2B1B0
B(2)= 4B2 + 2B1B0

まず A(x),B(x) に x = −d+1,−d+2,…, d−2,d−1 を代入した (2d−1)×2 個の数を求める.

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

手順 2. 値での乗算

C(-2) A(-2) ×B(-2)
C(-1) A(-1) ×B(-1)
C(0) A(0) ×B(0)
C(1) A(1) ×B(1)
C(2) A(2) ×B(2)

x に同じ値 r を代入した A(r),B(r) 同士の積 C(r) = A(r) × B(r) を求める.

多倍長同士の乗算が複数回行われているが,演算量としては N/d 桁同士の乗算が 2d−1 回なので (2d−1)M(N/d) だけの演算量になる.

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

C4 =(C(-2)−4C(-1)+6C(0)−4C(1)+C(2))/24
C3 =(−C(-2)+4C(-1)−4C(1)+C(2))/12
C2 =(−C(-2)+16C(-1)−30C(0)+16C(1)+C(2))/24
C1 =(C(-2)−8C(-1)+8C(1)−C(2))/12
C0 C(0)

ここで,C は 2d−2 次の多項式で表現されるはずなので,AB と同様に値を代入した場合,(d=3 で例示すると)

C(-2)= 16C4 − 8C3 + 4C2 − 2C1C0
C(-1)C4C3C2C1C0
C(0)C0
C(1)C4C3C2C1C0
C(2)= 16C2 + 8C2 + 4C2 + 2C1C0

という形で構成されているはずなので, この連立方程式を解析的に解いた左図の式に値を代入すれば各係数が求められる. この逆演算の式は d によって決まるので, プログラムに組み込んでおくことで高速に実行することができる. さらに,同じ演算になっている部分の処理を共通化することで さらに高速な実装を組むことができる.

手順 4.

CC4R4C3R3C2R2C1RC0

全体を繋げて多倍長数としての正規化を行う. この時,分割された各値は約 2N/d 桁あり,R(繋げる際にずらす桁数) は N/d 桁である.

計算量

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

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)

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

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

M(n) = O(n logd (2d-1))

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

定番の改良

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

C (∞) = A (∞) × B (∞) = Ad-1 Bd-1

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

一応補足しておくと,正確には∞を代入しているのではなく,

A’(x) = \frac{1}{x%5e{d-1}}A(x), B’(x) = \frac{1}{x%5e{d-1}}B(x)
C’(x) = A’(x) × B’(x) = \frac{1}{x%5e{2d-2}}C(x)

として,x→∞ と飛ばした計算である.

定数例

掛け合わせる多倍長数の桁数が同程度になるとは限らない. そのため,R によっては分割数がバラバラになることもある. ここではその具体的な例をいくつか挙げる. この係数は分割数依存で決定されるためプログラムとしてはハードコーディングで行うことになる. そのため,0 との乗算については省略することで高速に処理できるようになる.

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

C2
C1
C0
001
-111
010
A(-1)×B(-1)
A(0)×B(0)
A(∞)×B(∞)
,  
A(-1) B(-1)
A(0) B(0)
A(∞) B(∞)
-11
01
10
A1 B1
A0 B0

3 分割 × 2 分割

C3
C2
C1
C0
\frac{1}{2}
0002
1-210
-101-2
0200
A(-1)×B(-1)
A(0)×B(0)
A(1)×B(1)
A(∞)×B(∞)
,  
A(-1)
A(0)
A(1)
A(∞)
1-11
001
111
100
A2
A1
A0
B(-1)
B(0)
B(1)
B(∞)
-11
01
11
10
B1
B0

3 分割 × 3 分割

C4
C3
C2
C1
C0
\frac{1}{6}
00006
-13-3112
03-63-6
1-632-12
00600
A(-2)×B(-2)
A(-1)×B(-1)
A(0)×B(0)
A(1)×B(1)
A(∞)×B(∞)
,  
A(-2)B(-2)
A(-1)B(1)
A(0)B(0)
A(1)B(1)
A(∞)B(∞)
4-21
1-11
001
111
100
A2B2
A1B1
A0B0

4 分割 × 2 分割

C4
C3
C2
C1
C0
\frac{1}{6}
00006
-13-3112
03-63-6
1-632-12
00600
A(-2)×B(-2)
A(-1)×B(-1)
A(0)×B(0)
A(1)×B(1)
A(∞)×B(∞)
,  
A(-2)
A(-1)
A(0)
A(1)
A(∞)
-84-21
-11-11
0001
1111
1000
A3
A2
A1
A0
B(-2)
B(-1)
B(0)
B(1)
B(∞)
-21
-11
01
11
10
B1
B0

4 分割 × 3 分割

C5
C4
C3
C2
C1
C0
\frac{1}{24}
0000024
1-46-410
-240-42-120
-116-3016-10
2-16016-296
0024000
A(-2)×B(-2)
A(-1)×B(-1)
A(0)×B(0)
A(1)×B(1)
A(2)×B(2)
A(∞)×B(∞)
,  
A(-2)
A(-1)
A(0)
A(1)
A(2)
A(∞)
-84-21
-11-11
0001
1111
8421
1000
A3
A2
A1
A0
B(-2)
B(-1)
B(0)
B(1)
B(2)
B(∞)
4-21
1-11
001
111
421
100
B2
B1
B0

4 分割 × 4 分割

C6
C5
C4
C3
C2
C1
C0
\frac{1}{360}
000000360
-4-10120180-406180
015-6090-6015-1800
2020-540-9001400-900
0-15240-450240-151440
-16-1024072080-6720
000360000
A(-1/2)×B(-1/2)
A(-2)×B(-2)
A(-1)×B(-1)
A(0)×B(0)
A(1)×B(1)
A(2)×B(2)
A(∞)×B(∞)
,  
A(-1/2)B(-1/2)
A(-2)B(-2)
A(-1)B(-1)
A(0)B(0)
A(1)B(1)
A(2)B(2)
A(∞)B(∞)
-12-48
-84-21
-11-11
0001
1111
8421
1000
A3B3
A2B2
A1B1
A0B0