Skip to content →

乗算:Karatsuba, Toom-Cook, FFT

多項式同士の乗算を高速に計算するアルゴリズムである、Karatsuba乗算、Toom-Cook法、FFT を紹介します(正確には FFT を途中で使うアルゴリズム)。

解くのは次の問題です:

問題:
\(\mathbb{Z}/p\mathbb{Z}\) 係数(\(p\) : 素数)の多項式 $$f(x) = \sum_{i=0}^n a_ix^i,$$$$g(x) = \sum_{i=0}^n b_ix^i$$ が与えられます。二つの積 \(f(x)g(x)\) を求めてください。

Karatsuba乗算は \(O(n^{\log_2(3)})\), Toom-Cook法は \(O(n^{1+\varepsilon})\), FFTは \(O(n\log(n))\) で \(f(x)g(x)\) が求まります。

Karatsuba乗算

\(f(x)=a_0+a_1x,g(x)=b_1+b_1x\) の積を計算するには、素朴には \(a_0b_0,a_0b_1, a_1b_0,a_1b_1\) の4つの積を計算する必要があります。しかし、次のように変形すると、$$\begin{align} &(a_0+a_1x)(b_0+b_1x) \\ =& a_0b_0 + (a_0b_1 + a_1b_0)x + a_1b_2x^2 \\ =& a_0b_0 + \{(a_0 + a_1)(b_0 + b_1) – a_0b_0- a_1b_1\} x + a_1b_1x^2\end{align}$$\(a_0b_0,(a_0+a_1)(b_0+b_1),a_1b_1\) の3つで済みます。\(f(x),g(x)\) が \(n\) 次式のときは \( n/2 \) 次式に分けて再帰的に掛け算します。この計算方法をKaratsuba乗算といいます。\(n\) 次式の掛け算の計算量を \(T(n)\) とすると、$$T(n)=3T(n/2)+O(n)$$です。\(3T(n/2)\) が3つの乗算、\(O(n)\) が1つの加算です。\(n^c = 3(n/2)^c \) を解くと \(c=\log_2(3) > 1\) なので、$$T(n) = O(n^{\log_2(3)}) \subset O(n^{1.59})$$ を得ます。愚直 \(O(n^2)\) より高速になりました! Karatsubaは、Kolmogorovが「乗算は \(\Theta(n^2)\) が最適だろう」と予想したのを聞いて、一週間でこのアルゴリズムを思いついたと振り返っています。

Karatsuba乗算は2次式の多点評価、各点積、多項式補間で理解できます。
\(h(x)=c_0+c_1x+c_2x^2\) は多点評価 \(h(0),h(1),h(\infty)\) から $$h(x)=h(0)+\{h(1)-h(0)-h(\infty)\}x+h(\infty)x^2$$と多項式補間できます(\(h(\infty)\) は \(x^2\) の係数)。\(h(x)=f(x)g(x)\) とすると、\(x=0,1,\infty\) での多点評価は、\(f(0),f(1),f(\infty),g(0),g(1),g(\infty)\) を個別に求めてから、それぞれ掛ければよいです(\(f(\infty), g(\infty)\) は \(x\) の係数)。これは、Karatsuba乗算そのものです。

ところで、”2”次式の多点評価、各点積、多項式補間が最適なのでしょうか。3次式、4次式、5次式、…etcだと、どうなるのでしょうか(→Toom-Cook法, FFT)。また、評価点は \(x=0,1,\infty\) が最良でしょうか(→FFT)。

演習問題

添え字を bitwise or / and / xor で和を取った時の畳み込みにKaratsuba乗算を適用して高速化してください。つまり、\(\sum_{i|j=k} a_ib_j x^k, \sum_{i \& j=k} a_ib_j x^k, \sum_{i \oplus j=k} a_ib_j x^k\)を求めて下さい。

解答:\(f = a + bx , g = c + dx\) のとき
bitwise or 畳み込み : \(ac+((a+b)(c+d)-ac)x\)
bitwise and 畳み込み : \((a+b)(c+d)-bd+bdx\)
bitwise xor 畳み込み : \((ac+bd)(1-x)+(a+b)(c+d)x\)

bitwise or / and 畳み込みはゼータ変換・メビウス変換により高速化、bitwise xor 畳み込みは FFT により高速化できます。

熨斗袋さんが問題提起していたことですが、2 が非可逆の場合、後述のToom-Cook法で高速化できるのか気になっています。未解決です。

Toom-Cook法

Karatsuba法は2次式の多点評価、各点積、多項式補間でした。次数を上げて同じことをしたのがToom-Cook法です。

\(2d\) 次式で同じことをします。\(x=x_0,x_1,\ldots,x_d\) での多点評価は次の行列を、多項式の係数からなるベクトルに掛けることに相当します。
$$\begin{pmatrix}
x_0^0 & x_0^1 & x_0^2 & \cdots & x_0^{2d} \\
x_1^0 & x_1^1 & x_1^2 & \cdots & x_1^{2d} \\
x_2^0 & x_2^1 & x_2^2 & \cdots & x_2^{2d} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
x_{2d}^0 & x_{2d}^1 & x_{2d}^{2} & \cdots & x_{2d}^{2d}
\end{pmatrix} \begin{pmatrix}
h_0 \\ h_1 \\ h_2 \\ \vdots \\ h_{2d}
\end{pmatrix} = \begin{pmatrix}h(x_0) \\ h(x_1) \\ h(x_2) \\ \vdots \\ h(x_{2d})\end{pmatrix}$$この行列はヴァンデルモンド行列と呼ばれ、\(2d+1\) 点が全て異なるならば逆行列が存在します。多項式補間は逆行列を掛けることに相当します。\(d\) を定数とすると、行列を掛けることや逆行列を求めることは定数時間でできます。多点評価、各点積、多項式補間を \(2d\) 次式ですることにします。

\(n\) 次式 \(f(x)\) と \(g(x)\) の積を求めるアルゴリズムは次のようになります。

  1. \(n\) 次式 \(f(x)\) を \(d+1\) 個の \(n/(d+1)\) 次以下の式 \(f(x)=f_0(x)+f_1(x)x^{\lceil n/(d+1) \rceil}+f_2(x)x^{2 \lceil n/(d+1) \rceil}+\cdots+f_{d}(x)x^{d\lceil n/(d+1) \rceil}\) に分けます。
  2. \(n\) 次式 \(g(x)\) を \(d+1\) 個の \(n/(d+1)\) 次以下の式 \(g(x)=g_0(x)+g_1(x)x^{\lceil n/(d+1) \rceil}+g_2(x)x^{2\lceil n/(d+1) \rceil}+\cdots+g_{d}(x)x^{d\lceil n/(d+1) \rceil}\) に分けます。
  3. 先の行列を用いて、\(f(x)\) と \(g(x)\) の積を \(2d\) 次式の多点評価、各点積、多項式補間に帰着します。各点積では、\(2d+1\) 個の高々 \(n/(d+1)\) 次の多項式の積が計算されています。再帰的にこの各点の多項式の積を求めます。

\(n\) 次式の乗算をこのアルゴリズムに従って行うときの計算量を \(T(n)\) とすると$$T(n)=(2d+1)T(n/(d+1))+O(n)$$が成り立ち、解くと \(T(n)=O(n^{\log_{d+1}(2d+1)})\) です。\(d\) について、\(\log_{d+1}(2d+1)\) は単調減少で、\(d\to\infty\) のとき、\(\log_{d+1}(2d+1)\to 1\) です。つまり、任意の \(\varepsilon>0\) について、\(O(n^{1+\varepsilon})\) で多項式乗算できるアルゴリズムが手に入りました。

FFT

FFTでは多点評価する点をうまく選ぶことで、ヴァンデルモンド行列とその逆行列を高速で掛け算できます。

\(0\) で高次の項を埋めて、オーダーを変えずに掛け算する多項式の次数が \(2^m\) で表せると仮定できます。\(1\) の原始 \(2^{m+1}\) 乗根を \(\zeta\) と置いて、多点評価する \(n=2^{m+1}\) 点を \(x=\zeta^0,\zeta^1,\ldots,\zeta^{n-1}\) とします。

\(\sum_{i=0}^{n-1} \zeta^{ai} \zeta^{-bi} = n\delta_{a,b}\) なので、ヴァンデルモンド行列 \((\zeta^{ij})\) の逆行列は \((\frac{1}{n}(\zeta^{-1})^{ij})\) です。 多項式補間をするには、\(\zeta\) を \(\zeta^{-1}\) で置き換えて多点評価をしてから、\(n\) で割ればよいです。従って、ヴァンデルモンド行列を高速に掛けられることを示せば十分です。

因数定理より点 \(a\) での評価は \(f(a) = (f(x) \bmod (x-a))\) と書けます。つまり、\(i=0,1,\ldots,n-1\) に対し $$f(\zeta^i) = f(x) \bmod (x-\zeta^i)$$ です。mod を取る多項式を $$x^n-1=(x-\zeta^0)(x-\zeta^1)\cdots(x-\zeta^{n-1})$$ から半分ずつ因数を落としていくことで、高速化します。

中国剰余定理に従って、
\begin{align}
k[x] / (x^n-1) &\cong k[x] / (x^{n/2}-1) \times k[x] / ((\zeta^{-1} x)^{n/2} – 1) \\
& \cong k[x] / (x^{n/2}-1) \times k[\zeta x] / (x^{n/2} – 1)
\end{align}
という分解を繰り返せばよいです。通常、多項式の mod は掛け算以上に計算量が掛かりますが、FFTでは単純な形をしているので簡単です。\(f \bmod (x^{n/2}-1)\) は \(x^{n/2}=1\) をただ代入すればよいです。以上を再帰的に繰り返すと、計算量は \(O(\sum_{i} 2^i (n/2^i)) = O(n\log(n))\) となります。

\(k[x^2] / (x^n-1) \cong k[x] / (x^{n/2}-1)\) と \(f(x) = f_e(x^2) + xf_o(x^2)\) を用いてもよいです。

ところで、ヴァンデルモンド行列 \((\zeta^{ij})\) の逆行列は先の \(O(n\log(n))\) のアルゴリズムを一つずつ逆行することでも得られます。\(A=f \bmod (x^{2^k}-a^{2^{k}})\) と \(B=f \bmod (x^{2^k}-(\zeta^{n/2^{k+1}}a)^{2^k})\) から \(C=f \bmod (x^{2^{k+1}}-a^{2^{k+1}})\) を復元すればよいです。ぐっとにらむと $$C=\frac{1}{2}\{A+B+(a^{-1}x)^{2^k}(A-B)\}$$ です。逆行列は \((\frac{1}{n} \zeta^{-ij})\) の形でしたから、$$C=A+B+(ax)^{2^k}(A-B)$$ とすると、\((\zeta^{ij})\) を掛ける別のアルゴリズムが得られます。この二つのアルゴリズムの手順は異なります。先のアルゴリズムが周波数間引き、こちらのアルゴリズムが時間間引きと呼ばれています。

参考文献

  1. Karatsuba, Anatolii Alexeevich. “The complexity of computations.” Proceedings of the Steklov Institute of Mathematics-Interperiodica Translation 211 (1995): 169-183.
  2. Fiduccia, Charles M. “Polynomial evaluation via the division algorithm the fast Fourier transform revisited.” Proceedings of the fourth annual ACM symposium on Theory of computing. 1972.

Published in データ構造

Comments

コメントを残す

メールアドレスが公開されることはありません。

CAPTCHA