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

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

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

$$u(x)=u_0+u_1x+ \cdots +u_{n-1}x^{n-1}$$ (4.40)

に$n$個の値 $x_0,x_1,\ldots,x_{n-1}$を代入すれば、$n$個の値 $y_0,y_1, \ldots ,y_{n-1}$が得られる。行列を用いて表現すれば

$$\left( \begin{array}{@{ }c@{ }} y_0 \ y_1 \ \vdots \ y_... ...{ }c@{ }} u_0 \ u_1 \ \vdots \ u_{n-1} \end{array}\right)$$ (4.41)

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

$$V(x_0,x_1, \ldots ,x_{n-1})= \left( \begin{array}{@{ }cccc@... ...s \\ 1 & x_{n-1} & \ldots & x_{n-1}^{n-1} \end{array}\right)$$ (4.42)

と表され、行列式は

$$\mathrm{det}V(x_0,x_1,\ldots,x_{n-1})=\prod_{i>j}(x_i-x_j)$$ (4.43)

で与えられるので、$x_k$が全て異なれば逆行列が存在する。

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

1. $V(x_0,x_1,\ldots,x_{n-1})$ を $\mathbf{u},\mathbf{v}$ に掛け、 $\mathbf{u'},\mathbf{v'}$ を得る。
2. $\mathbf{u'},\mathbf{v'}$ を要素ごとに掛け、$\mathbf{w'}$ を得る。
3. $V(x_0,x_1,\ldots,x_{n-1})^{-1}$ を $\mathbf{w'}$ に掛け、$\mathbf{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,\ldots,x_{n-1})$の乗算をいかに速く行うかが問題となるが、$u,v$の分割数を大きくするとこの部分が複雑になり、位数も余り下がらなくなるので、大した効果は期待できなくなる。ここで発想を変えて、1.や3.で再帰を行うのが次に述べる高速Fourier変換である。

4.2 高速Fourier変換

Vandermonde行列の要素を

$$x_k=\omega_k=e^{2k\pi i/n}$$ (4.44)

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

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

$$u(x)=u_0+u_1x+ \cdots +u_{n-1}x^{n-1}$$ (4.45)

とおくと、 $\mathrm{DFT}_n(u)$を得るためには$u(\omega_k)$を評価しなければならないが、もし$n=2m$ならば$u_k$を$k$の偶奇で分けた2つの$m-1$次多項式

$$u_0(x) = u_0+u_2x^2+ \cdots +u_{2m-2}u^{m-1}$$ (4.46)

$$u_1(x) = u_1+u_3x^2+ \cdots +u_{2m-1}u^{m-1}$$ (4.47)

を用いることによって

$$u(\omega_k) = u_0(\omega_k^2)+\omega_ku_1(\omega_k^2)$$

$$= u_0(\omega_{2k})+\omega_ku_1(\omega_{2k})$$ (4.48)

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

ここで $k=0,1, \dots ,2m-1$について$u(\omega_k)$の値を求めることを考えると、$\omega_k$は$2m$個の値をとるのに対して、$\omega_{2k}$は$m$個の値しかとらず、しかもそれは複素$m$乗根であることが分かる。つまり $\mathrm{DFT}_n(u)$を得る問題は $\mathrm{DFT}_m(u_0),\mathrm{DFT}_m(u_1)$を得る問題に帰着されるのである。具体的に$n=4$の場合を考えると、

$$\begin{eqnarray*} u(\omega_0)&=&u_0(\omega_0)+\omega_0u_1(\omega_0) \\ u(\omega... ...(\omega_0) \\ u(\omega_3)&=&u_0(\omega_2)-\omega_1u_1(\omega_2) \end{eqnarray*}$$

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

$$T(n) = 2T(n/2)+cn$$

$$= nT(1)+cn \lg n$$

$$= \Theta(n \lg n)$$ (4.49)

となる。このアルゴリズムを再帰的FFT(Recursive FFT)という。

後は逆変換(Inverse DFT, $\mathrm{DFT}^{-1}$)が高速に実行できれば良いが、

$$V(\omega_0,\omega_1,\ldots,\omega_{n-1})^{-1}= \frac{1}{n} \... ...{n-1}^{-1} & \ldots & \omega_{n-1}^{-(n-1)} \end{array}\right)$$ (4.50)

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

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

1. FFTを用いて $\mathbf{u'}=\mathrm{DFT}_{n}(u),\mathbf{v'}=\mathrm{DFT}_{n}(v)$ を求める。必要な時間は $\Theta(n \lg n)$。
2. $\mathbf{u'},\mathbf{v'}$ の要素同士を掛けて $\mathbf{w'}$ を求める。必要な時間は $\Theta(n)$。
3. FFTを用いて $w=\mathrm{DFT}_n^{-1}(\mathbf{w'})$ を求める。必要な時間は $\Theta(n \lg n)$。

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

4.3 Cooley-Tukeyのアルゴリズム

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

$$\begin{array}{c} (0,1,2,3,4,5,6,7) \\ (0,2,4,6),(1,3,5,7) \\... ...,6),(1,5),(3,7) \\ (0),(4),(2),(6),(1),(5),(3),(7) \end{array}$$

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