離散フーリエ変換を用いた多倍長乗算の話

離散フーリエ変換のアルゴリズムを使って多倍長乗算ができるわけを書いた文章です。月刊電脳倶楽部 128 号 (1999 年 1 月号) の読み物横町のコーナーに収録されたものを HTML 化して加筆・修正しました。

はじめに
筆算の方法を用いた多倍長乗算
分割統治法を用いた多倍長乗算
離散フーリエ変換
高速フーリエ変換
FFT を用いた多倍長乗算
2 組の実数列の同時 FFT
2 組の実数列の同時 FFT を用いた多倍長乗算
他の多倍長演算のアプローチ
おわりに
参考文献
関連リンク
1. FFT
2. 多倍長計算

1. はじめに

CPU のレジスタに入り切らないような巨大あるいは高精度な数値データの演算を行うとき、四則演算をやるにもそのためのルーチンを用意しなければならない。^{注 1}

ここでは、「CPU のレジスタではちょっとはみ出す」という程度の数値データのことは考えない。何万桁、何百万桁という普通では目にしないような巨大な数値データを扱う場合に掛け算をどうすればよいかについて考える。^{注 2}

はっきり言ってそんな巨大な数が扱えたところで日常生活にこれといった恩恵があるわけではないが、敢えて常識を覆すことによって掛け算の本質が見えてくるかも知れない……などと論文調にこじつけてみたけれど、実は単なる好奇心だったりする。ただ、常識が覆されるのは事実だ。私は実際に巨大な数の計算を必要としているのではなくて、常識を覆えすことに快感を覚えるのである。^{注 3}

先に断っておくと、以下の解説を全部理解するにはそれなりの数学の知識が必要だと思われる。数式も出てくるので、数学が苦手な人は注意されたし。細かい定理や法則までは解説しないので、わからない点は各自で調べていただきたい。私だって、こんなことを学校で習ったわけではない。好奇心からその手の本を読み漁って覚えた程度なので、専門に勉強している人からのご指摘はお手柔らかに願いたい。^{注 4}

2. 筆算の方法を用いた多倍長乗算

筆算の方法で $N$ 桁の整数同士の乗算を行うには、だいたい $N^2$ に比例する時間が必要である。

$B$ 進法 ( $B$ は $2$ 以上の任意の整数) で表現したとき $N$ 桁になる整数 $X$ および $Y$ を、それぞれ次のように表現する。

$\begin{eqnarray} X & = & x_{N-1}B^{N-1}+x_{N-2}B^{N-2}+\cdots+x_2 B^2+x_1 B+x_0 \\ & = & \sum_{n=0}^{N-1}{x_n B^n} \\ Y & = & y_{N-1}B^{N-1}+y_{N-2}B^{N-2}+\cdots+y_2 B^2+y_1 B+y_0 \\ & = & \sum_{n=0}^{N-1}{y_n B^n} \end{eqnarray}$

ここで、 $x_n,y_n;n=0,1,\cdots,N-1$ は $B$ 進法で $1$ 桁の数とする。 $x_0$ は $B^0=1$ の位、 $x_1$ は $B^1=B$ の位、 $x_2$ は $B^2$ の位である。

$X$ と $Y$ の積を筆算の方法で求めると、次のようになる。

$\begin{eqnarray} X\times Y & = & \left\{x_{N-1}B^{N-1}+x_{N-2}B^{N-2}+\cdots+x_2 B^2+x_1 B+x_0\right\} \\ & \times & \left\{y_{N-1}B^{N-1}+y_{N-2}B^{N-2}+\cdots+y_2 B^2+y_1 B+y_0\right\} \\ & = & \sum_{n=0}^{N-1}{\left\{\sum_{k=0}^{N-1}{x_k B^k}\right\}y_n B^n}\\ & = & \sum_{n=0}^{N-1}{\sum_{k=0}^{N-1}{x_k y_n B^{k+n}}} \end{eqnarray}$

この結果、筆算の方法では $B$ 進法で $N$ 桁の整数同士の乗算結果を求めるために $B$ 進法で $1$ 桁の乗算が $N^2$ 回必要であることがわかる。筆算の方法を使った乗算のコストが桁数の $2$ 乗に比例するのはこのためである。

3. 分割統治法を用いた多倍長乗算

分割統治法を用いると、 $N$ 桁の整数同士の乗算を $N^{\log_2{3}}\approx N^{1.585}$ に比例する時間で行うことができる。これはディジタル法とも呼ばれる。^{注 5}

$B$ 進法で表現すると $2$ 桁になる数 $X$ および $Y$ があるとき、それぞれ下の位を $x_0$ と $y_0$ 、上の位を $x_1$ と $y_1$ とおく。すなわち、

$\begin{eqnarray} X & = & x_1 B+x_0 \\ Y & = & y_1 B+y_0 \end{eqnarray}$

とする。すると $X$ と $Y$ の積は

$\begin{eqnarray} X\times Y & = & (x_1 B+x_0)\times(y_1 B+y_0) \\ & = & x_1 y_1 B^2+(x_1 y_0+x_0 y_1)B+x_0 y_0\hspace{3em}\text{(1)} \end{eqnarray}$

となり、この式の計算には $B$ 進法で $1$ 桁の数同士の乗算が $x_1 y_1,x_1 y_0,x_0 y_1,x_0 y_0$ の $4$ 回必要になると思われる。

ところが、実は $B$ 進法で $1$ 桁の乗算を $3$ 回で済ませることができる。そのためには、 $\text{(1)}$ の右辺で乗算を $2$ 回使っている第 $2$ 項の係数を次のように変形すればよい。

$\begin{eqnarray} x_1 y_0+x_0 y_1 & = & -(-x_1 y_0-x_0 y_1) \\ & = & x_1 y_1+x_0 y_0-(x_1 y_1-x_1 y_0-x_0 y_1+x_0 y_0) \\ & = & x_1 y_1+x_0 y_0-(x_1-x_0)(y_1-y_0)\hspace{3em}\text{(2)} \end{eqnarray}$

これでは乗算が $2$ 回から $3$ 回に増えてしまったように見えるが、実際には $x_1 y_1$ と $x_0 y_0$ は $\text{(1)}$ の第 $1$ 項と第 $3$ 項にあるので、実質的に乗算は $1$ 回で済む。また、 $x_1-x_0$ と $y_1-y_0$ は負になってしまう場合があるが、符号の管理に必要なコストは乗算のコストと比較すれば無視できる。

$\text{(2)}$ を $\text{(1)}$ に代入すると、

$\begin{eqnarray} X\times Y & = & (x_1 B+x_0)\times(y_1 B+y_0) \\ & = & x_1 y_1 B^2+(x_1 y_0+x_0 y_1)B+x_0 y_0 \\ & = & \underbrace{x_1 y_1}_{\text{(a)}}B^2+\left\{(\underbrace{x_1 y_1}_{\text{(a)}}+\underbrace{x_0 y_0}_{\text{(b)}})-\underbrace{(x_1 - x_0)(y_1-y_0)}_{\text{(c)}}\right\}B+\underbrace{x_0 y_0}_{\text{(b)}} \end{eqnarray}$

となる。乗算は $\text{(a),(b),(c)}$ で示した $3$ 回だけで済んでいる。

すなわち、 $B$ 進法で表現したとき $2$ 桁になる整数同士の乗算は、 $B$ 進法で $1$ 桁の整数同士の乗算を $3$ 回と、それよりも十分に少ない時間で済む加減算を用いて行うことができる。つまり、分割統治法を用いると、桁数を $2$ 倍にしたときに乗算のコストが $3$ 倍になるということだ。

$N=2^k$ と仮定したとき、 $B$ 進法で $N$ 桁の整数同士の乗算を分割統治法で行うと、 $B$ 進法で $1$ 桁の整数同士の乗算が $3^k$ 回必要となる。 $k=\log_2{N}$ より、この方法での $B$ 進法で $N$ 桁の整数同士の乗算のコストは、 $3^k=3^{\log_2{N}}=N^{log_2{3}}\approx N^{1.585}$ 程度で済むことになる。実際には $N$ が $2^k$ で表現できるとは限らないので、 $2^k$ で表現できないときは $N^{1.585}$ よりも少し余計にかかる。

なお、結果の桁数が整数の場合の半分でよい多倍精度浮動小数点数の乗算には $2$ 分割よりも $4$ 分割のほうが効率がよく、整数の場合の $\frac{5}{6}$ の時間で済む。^{注 6}

4. 離散フーリエ変換

$N$ 個の複素数からなる数列 $x_n;n=0,1,\cdots,N-1$ の離散フーリエ変換 (DFT; Discrete Fourier Transform) を $X_m;m=0,1,\cdots,N-1$ とする。 $x_n$ は $X_m$ の逆離散フーリエ変換である。

$\begin{eqnarray} X_m & = & \sum_{n=0}^{N-1}\omega_{N}^{-m n}x_n & \hspace{3em} & m=0,1,\cdots,N-1 \\ x_n & = & \frac{1}{N}\sum_{m=0}^{N-1}\omega_{N}^{m n}X_m & & n=0,1,\cdots,N-1 \end{eqnarray}$

ここで $\omega_{N}$ は $1$ の $N$ 乗根の $1$ つであり、次のように定義される。

$\begin{eqnarray} \omega_{N} & = & e^{2\pi\mathrm{i}/N} = & \cos{\frac{2\pi}{N}}+\mathrm{i}\cdot\sin{\frac{2\pi}{N}} \end{eqnarray}$

( $e$ は自然対数の底 (ネイピア数)、 $\mathrm{i}$ は虚数単位)

$\omega_{N}$ は回転子とも呼ばれ、次の性質を持つ。

$\begin{eqnarray} \omega_{N}^N & = & 1 \end{eqnarray}$

( $\omega_{N}^n;n=0,1,\cdots,N-1$ は、複素平面の単位円に内接し、点 $1$ を頂点の $1$ つとする、正 $N$ 角形の各頂点の座標を示す)

$\begin{eqnarray} \omega_{N}^{-n} & = & \overline{\omega_{N}^{n}} \end{eqnarray}$

( $\overline{z}$ は $z$ の共役複素数)

$N$ が偶数のとき

$\begin{eqnarray} \omega_{N}^2 & = & \omega_{N/2} \\ \omega_{N}^{N/2} & = & -1 \end{eqnarray}$

離散フーリエ変換は式の通りだと $X_m;m=0,1,\cdots,N-1$ をすべて求めるのに乗算を $N^2$ 回必要とするが、これを $N\log{N}$ に比例する程度の時間で済ませるアルゴリズムが存在する。それが高速フーリエ変換 (FFT; Fast Fourier Transform) である。

5. 高速フーリエ変換

5.1. Cooley-Tukey の時間間引き型 FFT

時間のインデックスを間引くので時間間引き型 FFT と呼ばれる。

$N$ を偶数に限る。

時間のインデックスを偶数と奇数に分ける。すなわち、 $n=0,1,\cdots,N-1$ を $n=2k;k=0,1,\cdots,\frac{N}{2}-1$ と $n=2k+1;k=0,1,\cdots,\frac{N}{2}-1$ に分ける。

周波数のインデックスを前半と後半に分ける。すなわち、 $m=0,1,\cdots,N-1$ を $m=j;j=0,1,\cdots,\frac{N}{2}-1$ と $m=\frac{N}{2}+j;j=0,1,\cdots,\frac{N}{2}-1$ に分ける。

$x_n$ のインデックスが偶数の要素を $p_k$ 、インデックスが奇数の要素を $q_k$ とする。

$\begin{eqnarray} p_k & = & x_{2k} & \hspace{3em} & k=0,1,\cdots,\frac{N}{2}-1 \\ q_k & = & x_{2k+1} & & k=0,1,\cdots,\frac{N}{2}-1 \end{eqnarray}$

$p_k,q_k$ のフーリエ変換を $P_j,Q_j$ とする。

$\begin{eqnarray} P_j & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}p_k} & \hspace{3em} & j=0,1,\cdots,\frac{N}{2}-1 \\ Q_j & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}q_k} & & j=0,1,\cdots,\frac{N}{2}-1 \end{eqnarray}$

周波数のインデックスが前半のとき

$\begin{eqnarray} X_j & = & \sum_{n=0}^{N-1}{\omega_{N}^{-j n}x_n} \\ & = & \underbrace{\sum_{k=0}^{N/2-1}{\omega_{N}^{-j(2k)}x_{2k}}}_{\text{偶数}}+ \underbrace{\sum_{k=0}^{N/2-1}{\omega_{N}^{-j(2k+1)}x_{2k+1}}}_{\text{奇数}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-j(2k)}p_k}+ \sum_{k=0}^{N/2-1}{\omega_{N}^{-j(2k)-j}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-j(2k)}p_k}+ \sum_{k=0}^{N/2-1}{\omega_{N}^{-j}\omega_{N}^{-j(2k)}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{2(-j k)}p_k}+ \omega_{N}^{-j}\sum_{k=0}^{N/2-1}{\omega_{N}^{2(-j k)}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\left(\omega_{N}^2\right)^{-j k}p_k}+ \omega_{N}^{-j}\sum_{k=0}^{N/2-1}{\left(\omega_{N}^2\right)^{-j k}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}p_k}+ \omega_{N}^{-j}\sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}q_k} \\ & = & P_j+\omega_{N}^{-j}Q_j \end{eqnarray}$

周波数のインデックスが後半のとき

$\begin{eqnarray} X_{N/2+j} & = & \sum_{n=0}^{N-1}{\omega_{N}^{-(N/2+j)n}x_n} \\ & = & \underbrace{\sum_{k=0}^{N/2-1}{\omega_{N}^{-(N/2+j)(2k)}x_{2k}}}_{\text{偶数}} +\underbrace{\sum_{k=0}^{N/2-1}{\omega_{N}^{-(N/2+j)(2k+1)}x_{2k+1}}}_{\text{奇数}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{N(-k)-j(2k)}p_k} +\sum_{k=0}^{N/2-1}{\omega_{N}^{N(-k)-j(2k)-N/2-j}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{N(-k)}\omega_{N}^{-j(2k)}p_k} +\sum_{k=0}^{N/2-1}{\omega_{N}^{N(-k)}\omega_{N}^{-j(2k)}\omega_{N}^{-N/2}\omega_{N}^{-j}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\left(\omega_{N}^N\right)^{-k}\omega_{N}^{2(-j k)}p_k} +\sum_{k=0}^{N/2-1}{\left(\omega_{N}^N\right)^{-k}\omega_{N}^{2(-j k)}\overline{\omega_{N}^{N/2}}\omega_{N}^{-j}q_k} \\ & = & \sum_{k=0}^{N/2-1}{1^{-k}\left(\omega_{N}^2\right)^{-j k}p_k} +\sum_{k=0}^{N/2-1}{1^{-k}\left(\omega_{N}^2\right)^{-j k}(\overline{-1})\omega_{N}^{-j}q_k} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}p_k} -\omega_{N}^{-j}\sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}q_k} \\ & = & P_j-\omega_{N}^{-j}Q_j \end{eqnarray}$

最初にインデックスのビット反転を行う。インデックスが前半の要素をインデックスが偶数の領域に、インデックスが後半の要素をインデックスが奇数の領域に移動させることで、入力と出力を同じ配列に割り当てることができる。

$\omega_{N}^{-j}$ は三角関数の値をあらかじめテーブルに展開しておいてテーブル参照で求める。 $1$ 回の操作で数列の長さが半分になる代わりに数列の数が $2$ 倍に増えるので、 $1$ 回あたりテーブル参照と加減算と乗算を合わせて $N$ に比例する程度の時間がかかる。それをおよそ $\log_{2}{N}$ 回繰り返すことになるので、高速フーリエ変換の所要時間は全体で $N\log_{2}{N}$ に比例する程度になる。

5.2. Sande-Tukey の周波数間引き型 FFT

周波数のインデックスを間引くので周波数間引き型 FFT と呼ばれる。

$N$ を偶数に限る。

時間のインデックスを前半と後半に分ける。すなわち、 $n=0,1,\cdots,N-1$ を $n=k;k=0,1,\cdots,\frac{N}{2}-1$ と $n=\frac{N}{2}+k;k=0,1,\cdots,\frac{N}{2}-1$ に分ける。

周波数のインデックスを偶数と奇数に分ける。すなわち、 $m=0,1,\cdots,N-1$ を $m=2j;j=0,1,\cdots,\frac{N}{2}-1$ と $m=2j+1;j=0,1,\cdots,\frac{N}{2}-1$ に分ける。

周波数のインデックスが偶数のとき

$\begin{eqnarray} X_{2j} & = & \sum_{n=0}^{N-1}{\omega_{N}^{-(2j)n}x_n} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\underbrace{\omega_{N}^{-(2j)k}x_k}_{\text{前半}}+\underbrace{\omega_{N}^{-(2j)(N/2+k)}x_{N/2+k}}_{\text{後半}}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\omega_{N}^{-(2j)k}x_k+\omega_{N}^{N(-j)-(2j)k}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\omega_{N}^{-(2j)k}x_k+\omega_{N}^{N(-j)}\omega_{N}^{-(2j)k}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-(2j)k}\left\{x_k+\omega_{N}^{N(-j)}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{2(-j k)}\left\{x_k+\left(\omega_{N}^{N}\right)^{-j}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\left(\omega_{N}^2\right)^{-j k}\left(x_k+1^{-j}x_{N/2+k}\right)} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}\left(x_k+x_{N/2+k}\right)} \\ \end{eqnarray}$

これは長さが半分の数列 $x_k+x_{N/2+k};k=0,1,\cdots,\frac{N}{2}-1$ のフーリエ変換に等しい。

周波数のインデックスが奇数のとき

$\begin{eqnarray} X_{2j+1} & = & \sum_{n=0}^{N-1}{\omega_{N}^{-(2j+1)n}x_n} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\underbrace{\omega_{N}^{-(2j+1)k}x_k}_{\text{前半}}+\underbrace{\omega_{N}^{-(2j+1)(N/2+k)}x_{N/2+k}}_{\text{後半}}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\omega_{N}^{-(2j+1)k}x_k+\omega_{N}^{N(-j)-N/2-(2j+1)k}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\left\{\omega_{N}^{-(2j+1)k}x_k+\omega_{N}^{N(-j)}\omega_{N}^{-N/2}\omega_{N}^{-(2j+1)k}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-(2j+1)k}\left\{x_k+\left(\omega_{N}^{N}\right)^{-j}\overline{\omega_{N}^{N/2}}x_{N/2+k}\right\}} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-(2j)k-k}\left(x_k+1^{-j}(\overline{-1})x_{N/2+k}\right)} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{-(2j)k}\omega_{N}^{-k}\left(x_k-x_{N/2+k}\right)} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N}^{2(-j k)}\omega_{N}^{-k}(x_k-x_{N/2+k})} \\ & = & \sum_{k=0}^{N/2-1}{\left(\omega_{N}^2\right)^{-j k}\omega_{N}^{-k}\left(x_k-x_{N/2+k}\right)} \\ & = & \sum_{k=0}^{N/2-1}{\omega_{N/2}^{-j k}\omega_{N}^{-k}\left(x_k-x_{N/2+k}\right)} \\ \end{eqnarray}$

これは長さが半分の数列 $\omega_{N}^{-k}\left(x_k-x_{N/2+k}\right);k=0,1,\cdots,\frac{N}{2}-1$ のフーリエ変換に等しい。

$X_{2k}$ を $x_k$ の領域に、 $X_{2k+1}$ を $x_{N/2+k}$ の領域に格納することで、入力の $x_j$ と出力の $X_m$ を同じ配列に割り当てることができる。要素の順序が変わってしまうので、最後にインデックスのビット反転を行う。

$\omega_{N}^{-k}$ は三角関数の値をあらかじめテーブルに展開しておいてテーブル参照で求める。 $1$ 回の操作で数列の長さが半分になる代わりに数列の数が $2$ 倍に増えるので、 $1$ 回あたりテーブル参照と加減算と乗算を合わせて $N$ に比例する程度の時間がかかる。それをおよそ $\log_{2}{N}$ 回繰り返すことになるので、高速フーリエ変換の所要時間は全体で $N\log_{2}{N}$ に比例する程度になる。

5.3. 基数を一般化する

準備

数列 $X_m;m=0,1,\cdots,N-1$ を数列 $x_n;n=0,1,\cdots,N-1$ のフーリエ変換とする。

$\begin{eqnarray} X_m & = & \sum_{n=0}^{N-1}{\omega_{N}^{-m n}x_n} & \hspace{3em} & m=0,1,\cdots,N-1 \end{eqnarray}$

話を簡単にするため、ここでは要素数 $N$ を基数 $K$ の累乗に限る。基数 $K$ は $2$ 以上の小さい自然数であり、 $2$ または $4$ が使われることが多い。

$\begin{eqnarray} N & = & K^d \end{eqnarray}$

分割

周波数のインデックス $m$ を $K$ で割った商を $j$ 、余りを $u$ とする。また、時間のインデックス $n$ を $\frac{N}{K}$ で割った商を $v$ 、余りを $k$ とする。

$\begin{eqnarray} m & = & K j+u & \hspace{3em} & j=0,1,\cdots,\frac{N}{K}-1, u=0,1,\cdots,K-1 \\ n & = & \frac{N}{K}v+k & \hspace{3em} & v=0,1,\cdots,K-1, k=0,1,\cdots,\frac{N}{K}-1 \end{eqnarray}$

$X_m$ の $m$ と $n$ を $K j+u$ と $\frac{N}{K}v+k$ で置き換える。

$\begin{eqnarray} X_m & = & \sum_{n=0}^{N-1}{\omega_{N}^{-m n}x_n} \\ X_{K j+u} & = & \sum_{k=0}^{N/K-1}{\sum_{v=0}^{K-1}{\omega_{N}^{-(K j+u)((N/K)v+k)}x_{(N/K)v+k}}} \\ & = & \sum_{k=0}^{N/K-1}{\sum_{v=0}^{K-1}{\omega_{N}^{-N j v-K j k-(N/K)u v-u k}x_{(N/K)v+k}}} \\ & = & \sum_{k=0}^{N/K-1}{\omega_{N}^{-K j k}\omega_{N}^{-u k}\sum_{v=0}^{K-1}{\omega_{N}^{-(N/K)u v}x_{(N/K)v+k}}} \\ & = & \sum_{k=0}^{N/K-1}{\omega_{N/K}^{-j k}\left\{\omega_{N}^{-u k}\sum_{v=0}^{K-1}{\omega_{K}^{-u v}x_{(N/K)v+k}}\right\}} \end{eqnarray}$

ここで

$\begin{eqnarray} x^{\prime}_{(N/K)u+k} & = & \omega_{N}^{-u k}\sum_{v=0}^{K-1}{\omega_{K}^{-u v}x_{(N/K)v+k}} \end{eqnarray}$

と置くと、

$\begin{eqnarray} X_{K j+u} & = & \sum_{k=0}^{N/K-1}{\omega_{N/K}^{-j k}x^{\prime}_{(N/K)u+k}} \end{eqnarray}$

となって、 $u$ を固定して考えると $X_{K j+u}$ は $x^{\prime}_{(N/K)u+k}$ のフーリエ変換に等しい。すなわち、周波数のインデックス $m$ を基数 $K$ で割った余り $u$ で分類して商 $j$ を新たな周波数のインデックスとすることで、要素数 $N$ の数列 $x_n$ のフーリエ変換 $1$ 回を要素数 $\frac{N}{K}$ の数列 $x^{\prime}_k$ のフーリエ変換 $K$ 回に置き換えることができる。

三角関数の計算

通常は同じ $N$ で何度もフーリエ変換を行うので、 $\omega_{N}^{k}$ の値をあらかじめ $N$ 通りすべて計算しておく。 $\mathrm{sin}$ 関数の値を $\frac{1}{4}$ 周期分計算して $\frac{5}{4}$ 周期分テーブルに展開しておけばよい。

配列変数の再利用

$x_{(N/K)v+k}$ から $x^{\prime}_{(N/K)u+k}$ を求めるとき、 $x^{\prime}_{(N/K)u+k}$ を $x_{(N/K)u+k}$ に格納することで $x_n$ と $x^{\prime}_k$ を同じ配列変数に割り当てることができる。 $u$ と $v$ は $0,1,\cdots K-1$ の $K$ 通りしかないので、 $k$ 毎に $K$ 個の $x_{(N/K)v+k}$ をすべて読み出してから $x^{\prime}_{(N/K)u+k}$ を計算して $x_{(N/K)u+k}$ に書き戻せばよい。

データの並べ替え

$x_n$ のフーリエ変換を $x^{\prime}_k$ のフーリエ変換に置き換えることを $d$ 回繰り返すと要素数が $1$ になる。要素数が $1$ の数列はフーリエ変換を行っても値が変化しないので、この時点で $x_n$ が入っていた配列変数には $X_m$ の値が並んでいることになる。ただし、そのままでは $X_{K j+u}$ を $x_{(N/K)u+j}$ から取り出さなければならないので、最後にデータを並べ替える。インデックスを $d$ 桁の $K$ 進数で書いたとき、 $K j+u$ の $u$ は右端、 $\frac{N}{K}u+j$ の $u$ は左端にある。右端の数字を左端に移動させることを桁数を減らしながら $d$ 回繰り返すことになるので、 $X_m$ は $x_n$ が入っていた配列変数のインデックスを $d$ 桁の $K$ 進数で書いて数字を左右逆順にした位置に格納されている。

所要時間

大雑把に見積もると、分割にかかる時間は $N$ に比例する程度であり、それを $d=\log_{K}{N}$ 回繰り返すので、この方法のフーリエ変換の所要時間は $N\log{N}$ に比例する程度になる。

基数 $2$ の FFT

$x^{\prime}_{(N/K)u+k}$ の式に $K=2$ を代入し、 $u=0,1$ と $v=0,1$ をそれぞれ展開する。

$\begin{eqnarray} x^{\prime}_{k} & = & x_{k}+x_{N/2+k} \\ x^{\prime}_{N/2+k} & = & \omega_{N}^{-k}(x_{k}-x_{N/2+k}) \end{eqnarray}$

最後にインデックスが $2$ 進数で左右逆順になるようにデータを並べ替える。

基数 $4$ の FFT

$x^{\prime}_{(N/K)u+k}$ の式に $K=4$ を代入し、 $u=0,1,2,3$ と $v=0,1,2,3$ をそれぞれ展開する。

ここでは基数 $2$ の FFT と共通のコードを使用できるように最後にインデックスが $4$ 進数で左右逆順になるように並べ替える代わりに $2$ 進数で左右逆順になるように並べ替えることにする。そのために $x^{\prime}_{N/2+k}$ と $x^{\prime}_{N/4+k}$ を入れ替えておく ( $4$ 進数の $1$ と $2$ は $2$ 進数で $01$ と $10$ である)。

$\begin{eqnarray} x^{\prime}_{k} & = & x_{k}+x_{N/4+k}+x_{N/2+k}+x_{3 N/4+k} \\ x^{\prime}_{N/2+k} & = & \omega_{N}^{-k}(x_{k}-\mathrm{i}\,x_{N/4+k}-x_{N/2+k}+\mathrm{i}\,x_{3 N/4+k}) \\ x^{\prime}_{N/4+k} & = & \omega_{N}^{-2 k}(x_{k}-x_{N/4+k}+x_{N/2+k}-x_{3 N/4+k}) \\ x^{\prime}_{3 N/4+k} & = & \omega_{N}^{-3 k}(x_{k}+\mathrm{i}\,x_{N/4+k}-x_{N/2+k}-\mathrm{i}\,x_{3 N/4+k}) \end{eqnarray}$

最後にインデックスが $2$ 進数で左右逆順になるようにデータを並べ替える。

6. FFT を用いた多倍長乗算

$N$ 個の複素数からなる数列 $x_n,y_n;n=0,1,\cdots,N-1$ の離散フーリエ変換を $X_m,Y_m;m=0,1,\cdots,N-1$ とする。すなわち、

$\begin{eqnarray} X_m & = & \sum_{n=0}^{N-1}{\omega_{N}^{-m n}x_n}\hspace{3em}m=0,1,\cdots,N-1 \\ Y_m & = & \sum_{n=0}^{N-1}{\omega_{N}^{-m n}y_n}\hspace{3em}m=0,1,\cdots,N-1 \end{eqnarray}$

とする。

ここで、数列 $Z_m$ を次のように定義する。

$\begin{eqnarray} Z_m & = & X_m\times Y_m\hspace{3em}m=0,1,\cdots,N-1 \end{eqnarray}$

$Z_m$ を $x_n$ と $y_n$ で表現すると、

$\begin{eqnarray} Z_m & = & X_m\times Y_m \\ & = & \left\{\sum_{n=0}^{N-1}{\omega_{N}^{-m n}x_n}\right\}\times\left\{\sum_{n=0}^{N-1}{\omega_{N}^{-m n}y_n}\right\} \\ & = & \left\{\omega_{N}^{-m 0}x_0+\omega_{N}^{-m 1}x_1+\cdots+\omega_{N}^{-m(N-1)}x_{N-1}\right\} \\ & \times & \left\{\omega_{N}^{-m 0}y_0+\omega_{N}^{-m 1}y_1+\cdots+\omega_{N}^{-m(N-1)}y_{N-1}\right\} \\ & = & \omega_{N}^{-m 0}\left\{x_0 y_0+x_1 y_{N-1}+\cdots+x_{N-1}y_1\right\} \\ & + & \omega_{N}^{-m 1}\left\{x_0 y_1+x_1 y_0+x_2 y{N-1}+\cdots+x_{N-1}y_2\right\} \\ & + & \omega_{N}^{-m 2}\left\{x_0 y_2+x_1 y_1+x_2 y_0+x_3 y{N-1}+\cdots+x_{N-1}y_3\right\} \\ & \vdots & \\ & + & \omega_{N}^{-m(N-2)}\left\{x_0 y_{N-2}+x_1 y_{N-3}+\cdots+x_{N-1}y_{N-1}\right\} \\ & + & \omega_{N}^{-m(N-1)}\left\{x_0 y_{N-1}+x_1 y_{N-2}+\cdots+x_{N-1}y_0\right\} \\ & = & \sum_{n=0}^{N-1}{\omega_{N}^{-m n}\sum_{k=0}^{N-1}{x_k y_{(n-k)\pmod{N}}}} \end{eqnarray}$

となる。最後の式に注目すると、これは

$\begin{eqnarray} z_n & = & \sum_{k=0}^{N-1}{x_k y_{(n-k)\pmod{N}}}\hspace{3em}n=0,1,\cdots,N-1 \end{eqnarray}$

で定義された数列 $z_n$ のフーリエ変換に等しい。すなわち、数列 $Z_m;m=0,1,\cdots,N-1$ は数列 $z_n;n=0,1,\cdots,N-1$ のフーリエ変換である。

ここで $N$ を偶数に限定して $N=2M$ とおき、次の条件を加える。

$\begin{eqnarray} x_n & = & y_n & = & 0\hspace{3em}n=M,M+1,\cdots,2M-1 \end{eqnarray}$

つまり、 $x_n$ と $y_n$ の“後半”をすべて $0$ にする。すると、 $z_n$ は次のようになる。

$\begin{array}{rclclccclcl} z_0 & = & x_0 y_0 & & & & & & & & \\ z_1 & = & x_0 y_1 & + & x_1 y_0 & & & & & & \\ & & & \vdots & & & & & & & \\ z_{M-2} & = & x_0 y_{M-2} & + & x_1 y_{M-3} & + & \cdots & + & x_{M-2} y_0 & & \\ z_{M-1} & = & x_0 y_{M-1} & + & x_1 y_{M-2} & + & \cdots & + & x_{M-2} y_1 & + & x_{M-1} y_0 \\ z_M & = & & & x_1 y_{M-1} & + & \cdots & + & x_{M-2} y_2 & + & x_{M-1} y_1 \\ & & & & & & & & & \vdots & \\ z_{2M-3} & = & & & & & & & x_{M-2} y_{M-1} & + & x_{M-1} y_{M-2} \\ z_{2M-2} & = & & & & & & & & & x_{M-1} y_{M-1} \\ z_{2M-1} & = & 0 & & & & & & & & \end{array}$

これはまさに $M$ 桁の多倍長乗算そのものである。言い換えると、 $B$ 進法で $M$ 桁の数値 $X$ と $Y$ の $B^n$ の位を数をそれぞれ $x_n,y_n;n=0,1,\cdots,M-1$ としたとき、

$\begin{eqnarray} X & = & \sum_{n=0}^{M-1}{x_n B^n} \\ Y & = & \sum_{n=0}^{M-1}{y_n B^n} \end{eqnarray}$

であり、 $X$ と $Y$ の積は、

$\begin{eqnarray} X\times Y & = & \sum_{n=0}^{2M-1}{z_n B^n} \end{eqnarray}$

である。

まとめると、 $B$ 進法で $M$ 桁の数 $X$ と $Y$ が与えられたとき、 $X$ と $Y$ をそれぞれ $B$ 進法で $1$ 桁の数を $M$ 個並べた数列と見なし、上位の側に $0$ を $M$ 個追加した $2M$ 個の数からなる数列を作る。 $X$ と $Y$ それぞれの数列を離散フーリエ変換した結果の同じ位置の数同士を掛け合わせ、逆フーリエ変換すると、 $X$ と $Y$ の積を $B$ 進法で表した数列が得られるのだ。

離散フーリエ変換は FFT (高速フーリエ変換) のアルゴリズムを用いるとだいたい $N\log{N}$ に比例する程度のコストで求めることができるので、FFT を用いると多倍長乗算を桁数 $N$ に対してだいたい $N\log{N}$ に比例する程度の時間で求めることができる。

なお、 $z_n$ は高々 $M B^2$ までの数だが、 $B$ 進法の $1$ 桁には収まっていないので、最後に大きすぎる分を上の位に繰り上げる必要がある。しかしそれは桁数 $N$ に比例する程度の時間で済む。

7. 2 組の実数列の同時 FFT

$2$ 組の実数列 $x_n,y_n;n=0,1,\cdots,N-1$ が与えられたとき、それぞれの離散フーリエ変換 $X_m,Y_m;m=0,1,\cdots,N-1$ を以下の手順で求めることができる。

複素数列 $t_n;n=0,1,\cdots,N-1$ を次のようにおく。

$\begin{eqnarray} t_n & = & x_n+\mathrm{i}\cdot y_n \end{eqnarray}$

( $x_n$ と $y_n$ はいずれも実数なので、 $t_n$ の実数部は $x_n$ 、虚数部は $y_n$ である)

FFT を用いて $t_n$ の離散フーリエ変換 $T_m;m=0,1,\cdots,N-1$ を求める。

$T_m$ から、 $x_n$ と $y_n$ の離散フーリエ変換 $X_m,Y_m$ を次の式で求める。

$\begin{eqnarray} X_m & = & \frac{\mathrm{Re}(T_m)+\mathrm{Re}(T_{N-m})}{2}+\mathrm{i}\cdot\frac{\mathrm{Im}(T_m)-\mathrm{Im}(T_{N-m})}{2} \\ Y_m & = & \frac{\mathrm{Im}(T_m)+\mathrm{Im}(T_{N-m})}{2}-\mathrm{i}\cdot\frac{\mathrm{Re}(T_m)-\mathrm{Re}(T_{N-m})}{2} \end{eqnarray}$

( $m=0$ のとき、 $T_{N-m}=T_N=T_0$ )

この方法を用いると、 $2$ 組の実数列の離散フーリエ変換を、 $1$ 組の複素数列の離散フーリエ変換とほぼ同じ時間で求めることができる。^[2]

8. 2 組の実数列の同時 FFT を用いた多倍長乗算

FFT を用いた多倍長乗算では、引数となる多倍長データを数列と見なしてそれぞれの離散フーリエ変換を行い、同じ番号の要素同士を掛け合わせて逆フーリエ変換を行う。このように引数の数列 $2$ つに対して離散フーリエ変換を行うが、それらの引数は実数列なので、 $2$ 組の実数列の同時 FFT を用いることができる。

引数の数列を $x_n,y_n;n=0,1,\cdots,N-1$ とおく。 $2$ 組の実数列の同時 FFT なので、

$\begin{eqnarray} t_n & = & x_n+\mathrm{i}\cdot y_n \end{eqnarray}$

とおき、FFT を用いて $t_n$ の離散フーリエ変換 $T_m$ を求める。

以降、式が冗長になるので次のような文字で置き換えて表記する。

$\begin{eqnarray} p & := & \mathrm{Re}(T_m) \\ q & := & \mathrm{Im}(T_m) \\ r & := & \mathrm{Re}(T_{N-m}) \\ s & := & \mathrm{Im}(T_{N-m}) \end{eqnarray}$

( $?=\mathrm{Re}(?)+\mathrm{i}\cdot\mathrm{Im}(?)$ 。 $\mathrm{Re}(?)$ は $?$ の実数部、 $\mathrm{Im}(?)$ は $?$ の虚数部)

$2$ 組の実数列の同時 FFT の関係から、 $X_m,Y_m;m=0,1,\cdots,N-1$ は次のようになる。

$\begin{eqnarray} X_m & = & \frac{\mathrm{Re}(T_m)+\mathrm{Re}(T_{N-m})}{2}+\mathrm{i}\cdot\frac{\mathrm{Im}(T_m)-\mathrm{Im}(T_{N-m})}{2} \\ & = & \frac{p+r}{2}+\mathrm{i}\cdot\frac{q-s}{2} \\ Y_m & = & \frac{\mathrm{Im}(T_m)+\mathrm{Im}(T_{N-m})}{2}-\mathrm{i}\cdot\frac{\mathrm{Re}(T_m)-\mathrm{Re}(T_{N-m})}{2} \\ & = & \frac{q+s}{2}-\mathrm{i}\cdot\frac{p-r}{2} \end{eqnarray}$

$Z_m=X_m\times Y_m$ を $p,q,r,s$ を使って展開すると、

$\begin{eqnarray} Z_m & = & X_m\times Y_m \\ & = & \left\{\frac{p+r}{2}+\mathrm{i}\cdot\frac{q-s}{2}\right\}\times\left\{\frac{q+s}{2}-\mathrm{i}\cdot\frac{p-r}{2}\right\} \\ & = & \frac{\left\{(p+r)+\mathrm{i}\cdot(q-s)\right\}\times\left\{(q+s)-\mathrm{i}\cdot(p-r)\right\}}{4} \\ & = & \frac{(p+r)(q+s)+(q-s)(p-r)}{4}+\mathrm{i}\cdot\frac{-(p+r)(p-r)+(q-s)(q+s)}{4} \\ & = & \frac{p q+p s+q r+r s+p q-q r-p s+r s}{4}+\mathrm{i}\cdot\frac{-p^2+r^2+q^2-s^2}{4} \\ & = & \frac{p q+r s}{2}+\mathrm{i}\cdot\frac{(r^2+q^2)-(p^2+s^2)}{4} \end{eqnarray}$

このように実数部が打ち消し合って簡単な式になる。表記を $p,q,r,s$ から $T_m$ に戻すと、

$\begin{eqnarray} Z_m & = & \frac{\mathrm{Re}(T_m)\mathrm{Im}(T_m)+\mathrm{Re}(T_{N-m})\mathrm{Im}(T_{N-m})}{2} \\ & + & \mathrm{i}\cdot\frac{\mathrm{Re}(T_{N-m})^2+\mathrm{Im}(T_m)^2-\mathrm{Re}(T_m)^2-\mathrm{Im}(T_{N-m})^2}{4} \end{eqnarray}$

となる。

$m=0$ の場合について考えると、

$\begin{eqnarray} Z_0 & = & \mathrm{Re}(T_0)\mathrm{Im}(T_0) \end{eqnarray}$

である。また、 $N$ が偶数のとき、 $m=\frac{N}{2}$ では、

$\begin{eqnarray} Z_{N/2} & = & \mathrm{Re}(T_{N/2})\mathrm{Im}(T_{N/2}) \end{eqnarray}$

である。

この式を使うと、 $x_n,y_n$ という $2$ 組の実数列の離散フーリエ変換を同時に行うだけでなく、それぞれの離散フーリエ変換 $X_m,Y_m$ を算出せずに直接 $Z_m=X_m\times Y_m$ を求めることになるので、FFT を用いた多倍長乗算の計算量を減らすことができる。^{注 7}

9. 他の多倍長演算のアプローチ

詳しいことは次の機会に改めて書きたいが、多倍長除算について少しだけ触れておこう。

何万桁というオーダーになったとき、多倍長の数値データ同士の除算をまともにやろうとするのは間違いである。多倍長の数値データ同士の除算は乗算よりも遥かに大きな計算量を必要とするので、なるべく除算をしなくて済むようにするべきだ。

しかし、どうしても除算が必要なときは、「Newton 法の反復で除数の逆数を求めてから被除数に掛ける」という手段を使う。^[1]

例えば、定数 $A$ の逆数を求めるには、適切な方法で初期値 $x_0$ を決定し、次の反復を行うとよい。

$\begin{eqnarray} x_{n+1} & = & 2 x_n-A x_n^2 \end{eqnarray}$

この反復は $\frac{1}{A}$ に $2$ 次の収束を示す ( $1$ 回の反復で有効桁数がほぼ $2$ 倍になる)。

多倍精度実数の平方根を求める場合にも Newton 法が有効だ。

当然、Newton 法の過程で必要になる多倍長の整数同士の乗算には FFT を使う。

10. おわりに

「X68k で円周率を計算しました」なんて言うと、理論をよく理解もせずにただ速いだけのパソコンを使っている……というよりパソコンに使われているような人に「そんな遅いパソコンで円周率を計算して何が楽しいのか」なんて言われて頭にくることがあるので、もしも円周率の計算に関して書くとしたら、プログラムを示すことよりも、可能な限り自分で理解した上で解説文を書くことに重点を置きたいと思う。

それにしても AGM (算術幾何平均) の理論 ^[7] は難しいのぉ。私にゃいまだによくわからん。

11. 参考文献

D.E.Knuth, 中川圭介訳 , KNUTH The Art of Computer Programming 第 4 分冊準数値算法／算術演算 , サイエンス社 , ISBN4-7819-0426-2
五十嵐善英 , 整数の乗算法 , 情報処理 1983 Apr. Vol.24 No.4 pp.358-362
小池慎一 , C による科学技術計算 , CQ 出版社 , ISBN4-7898-3032-2
奥村晴彦 , C 言語による最新アルゴリズム辞典 , 技術評論社 , ISBN4-87408-414-1
H.P. スウ , 佐藤平八訳 , 工学基礎演習シリーズ 1 フーリエ解析 , 森北出版 , ISBN4-627-93010-0
丹明彦 , そこに π があるから , Oh!X 1988.8 pp.51-60
JONATHAN M. BORWEIN/PETER B. BORWEIN, Pi and the AGM, A Wiley-Interscience Publication, ISBN0-471-83138-7