乗算：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)$ が求まります。

1 Karatsuba乗算
- 1.1 演習問題
2 Toom-Cook法
3 FFT
4 参考文献

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$ の３つで済みます。$f(x),g(x)$ が $n$ 次式のときは $ n/2 $ 次式に分けて再帰的に掛け算します。この計算方法をKaratsuba乗算といいます。$n$ 次式の掛け算の計算量を $T(n)$ とすると、$$T(n)=3T(n/2)+O(n)$$です。$3T(n/2)$ が３つの乗算、$O(n)$ が１つの加算です。$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乗算は２次式の多点評価、各点積、多項式補間で理解できます。
$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乗算そのものです。

ところで、”２”次式の多点評価、各点積、多項式補間が最適なのでしょうか。３次式、４次式、５次式、…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法は２次式の多点評価、各点積、多項式補間でした。次数を上げて同じことをしたのが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)$ の積を求めるアルゴリズムは次のようになります。

$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}$ に分けます。
$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}$ に分けます。
先の行列を用いて、$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})$ を掛ける別のアルゴリズムが得られます。この二つのアルゴリズムの手順は異なります。先のアルゴリズムが周波数間引き、こちらのアルゴリズムが時間間引きと呼ばれています。

参考文献

Karatsuba, Anatolii Alexeevich. “The complexity of computations.” Proceedings of the Steklov Institute of Mathematics-Interperiodica Translation 211 (1995): 169-183.
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.