━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 離散フーリエ変換を用いた多倍長乗算の話 鎌田　誠 ──────────────────────────────────── 　CPU のレジスタに入り切らないような巨大あるいは高精度な数値データの演算 を行うとき、四則演算をやるにもそのためのルーチンを用意しなければならない。(*1) 　ここでは、「CPU のレジスタではちょっとはみ出す」という程度の数値データ のことは考えない。何万桁、何百万桁という普通では目にしないような巨大な数 値データを扱う場合に掛け算をどうすればよいかについて考える。(*2) 　はっきり言ってそんな巨大な数が扱えたところで日常生活にこれといった恩恵 があるわけではないが、敢えて常識を覆すことによって掛け算の本質が見えてく るかも知れない……などと論文調にこじつけてみたけれど、実は単なる好奇心だ ったりする。ただ、常識が覆されるのは事実だ。私は実際に巨大な数の計算を必 要としているのではなくて、常識を覆えすことに快感を覚えるのである。(*3) 　先に断っておくと、以下の解説を全部理解するにはそれなりの数学の知識が必 要だと思われる。数式も出てくるので、数学が苦手な人は注意されたし。細かい 定理や法則までは解説しないので、わからない点は各自で調べていただきたい。 私だって、こんなことを学校で習ったわけではない。好奇心からその手の本を読 み漁って覚えた程度なので、専門に勉強している人からのご指摘はお手柔らかに 願いたい。(*4) (*1)そのような多倍長の数値データの四則演算を行うライブラ リとしては、電脳倶楽部でも Vol.26 と Vol.28 に LONG.C と いうプログラムが掲載されている。が、LONG.C では今回解説 するような高度なアルゴリズムは使われていない。 (*2)こんなの使うの、円周率の計算くらいだよね。 (*3)勿論、アルゴリズムの検証のために実際に計算してみるこ とは重要だ。 (*4) FFT で掛け算をやろうという時点で既に FFT の使い方を 間違っているという話もある。 ━─────────────────────────────────── ●筆算の方法を用いた多倍長乗算 　筆算の方法で N 桁の整数同士の乗算を行うには、だいたい N^2 に比例する時 間が必要である。 　B 進法（B は 2 以上の任意の整数）で表現したとき N 桁になる整数 X およ び Y を、それぞれ次のように表現する。 N-1 X = Σ x(n)・B^n n=0 N-1 Y = Σ y(n)・B^n n=0 ここで、x(n),y(n);n=0,1,…,N-1 は B 進法で 1 桁の数とする。 　X と Y の積を筆算の方法で求めると、次のようになる。 N-1 N-1 X*Y = Σ { Σ x(k)・B^k }*y(n)・B^n n=0 k=0 N-1 N-1 = Σ Σ x(k)*y(n)・B^(k+n) n=0 k=0 　この結果、筆算の方法では B 進法で N 桁の整数同士の乗算結果を求めるため に B 進法で 1 桁の乗算が N^2 回必要であることがわかる。筆算の方法を使っ た乗算のコストが桁数の 2 乗に比例するのはこのためである。 ━─────────────────────────────────── ●分割統治法を用いた多倍長乗算 　分割統治法を用いると、N 桁の整数同士の乗算を N^log2(3) ≒ N^1.585 に比 例する時間で行うことができる。これはディジタル法とも呼ばれる。(*5) 　B 進法で表現すると 2 桁になる数 X および Y があるとき、それぞれ下の位 を x0 と y0、上の位を x1 と y1 とおく。すなわち、 X = x1・B + x0 Y = y1・B + y0 とする。すると X と Y の積は X*Y = (x1・B + x0)*(y1・B + y0) = x1*y1・B^2 + (x1*y0 + x0*y1)・B + x0*y0 (1) となり、この式の計算には B 進法で 1 桁の数同士の乗算が x1*y1、x1*y0、 x0*y1、x0*y0 の 4 回必要になると思われる。 　ところが、実は B 進法で 1 桁の乗算を 3 回で済ませることができる。その ためには、(1) の右辺で乗算を 2 回使っている第 2 項の係数を次のように変形 すればよい。 x1*y0 + x0*y1 = -(- x1*y0 - x0*y1) = x1*y1 + x0*y0 - (x1*y1 - x1*y0 - x0*y1 + x0*y0) = x1*y1 + x0*y0 - (x1 - x0)*(y1 - y0) (2) 　これでは乗算が 2 回から 3 回に増えてしまったように見えるが、実際には x1*y1 と x0*y0 は (1) の第 1 項と第 3 項にあるので、実質的に乗算は 1 回 で済む。また、x1 - x0 と y1 - y0 は負になってしまう場合があるが、符号の 管理に必要なコストは乗算のコストと比較すれば無視できる。 　(2) を (1) に代入すると、 X*Y = (x1・B + x0)*(y1・B + y0) = x1*y1・B^2 + (x1*y0 + x0*y1)・B + x0*y0 = x1*y1・B^2 + {x1*y1 + x0*y0 - (x1 - x0)*(y1 - y0)}・B + x0*y0 (a) (a) (b) (c) (b) となる。乗算は (a),(b),(c) で示した 3 回だけで済んでいる。 　すなわち、B 進法で表現したとき 2 桁になる整数同士の乗算は、B 進法で 1 桁の整数同士の乗算を 3 回と、それよりも十分に少ない加減算を用いて行うこ とができる。つまり、分割統治法を用いると、桁数を 2 倍にしたときに乗算の コストが 3 倍になるということだ。 　N = 2^k と仮定したとき、B 進法で N 桁の整数同士の乗算を分割統治法で行 うと、B 進法で 1 桁の整数同士の乗算が 3^k 回必要となる。k = log2(N) より、 この方法での B 進法で N 桁の整数同士の乗算のコストは、3^k = 3^log2(N) = N^log2(3) ≒ N^1.585 程度で済むことになる。実際には N が 2^k で表現でき るとは限らないので、2^k で表現できないときは N^1.585 よりも少し余計にか かる。 　なお、結果の桁数が整数の場合の半分でよい多倍精度浮動小数点数の乗算には 2 分割よりも 4 分割のほうが効率がよく、整数の場合の 5/6 の時間で済む。 (*6) (*5)この方法は 1962 年に A.Karatsuba と Y.Ofman によって 発見された。簡単な方法であるにも関わらず 1962 年以前に報 告されていなかったというのは意外だ。 (*6)この事実は私が学生だった頃に自分で発見したものです。 ━─────────────────────────────────── ●離散フーリエ変換 　離散フーリエ変換というのは、本来、波形のサンプリングデータから周波数成 分を得るために使われる理論だ。ここではその詳細についての解説は省略し、離 散フーリエ変換の式を示すに留める。 　N 個の複素数からなる数列 x(n);n=0,1,…,N-1 の離散フーリエ変換を X(m);m=0,1,…,N-1 とする。すなわち、 N-1 X(m) = Σ ω^(-m*n)*x(n) m=0,1,…,N-1 n=0 （一般的に離散フーリエ変換には係数 1/N が付くが、ここでは省略する） と定義する。 　ここで、ω は 1 の N 乗根の 1 つであり、次のように定義される。 ω = e^(2πi/N) = cos(2π/N) + i・sin(2π/N) （e は自然対数の底、i は虚数単位） 　ω は回転子とも呼ばれ、次の性質を持つ。 ω^N = 1 （ω^n;n=0,1,…,N-1 は、複素平面の単位円に内接し、点 1 を頂点の 　1 つとする、正 N 角形の各頂点の座標を示す） 　離散フーリエ変換は式の通りだと X(m);m=0,1,…,N-1 をすべて求めるのに乗 算を N^2 回必要とするが、これを N*log(N) に比例する程度の時間で済ませる アルゴリズムが存在する。それが FFT（高速フーリエ変換）である。FFT につい ては電気や音声の専門書を見れば載っているので、各自で勉強していただきたい。 ━─────────────────────────────────── ● FFT を用いた多倍長乗算 　N 個の複素数からなる数列 x(n),y(n);n=0,1,…,N-1 の離散フーリエ変換を X(m),Y(m);m=0,1,…,N-1 とする。すなわち、 N-1 X(m) = Σ ω^(-m*n)*x(n) m=0,1,…,N-1 n=0 N-1 Y(m) = Σ ω^(-m*n)*y(n) m=0,1,…,N-1 n=0 とする。 　ここで、数列 Z(m) を次のように定義する。 Z(m) = X(m)*Y(m) m=0,1,…,N-1 　Z(m) を x(n) と y(n) で表現すると、 Z(m) = X(m)*Y(m) N-1 N-1 = { Σ ω^(-m*n)*x(n) }*{ Σ ω^(-m*n)*y(n) } n=0 n=0 = { ω^(-m*0)*x(0) + ω^(-m*1)*x(1) + … + ω^(-m*(N-1))*x(N-1) } * { ω^(-m*0)*y(0) + ω^(-m*1)*y(1) + … + ω^(-m*(N-1))*y(N-1) } = ω^(-m*0)*{ x(0)*y(0) + x(1)*y(N-1) + … + x(N-1)*y(1) } + ω^(-m*1)*{ x(0)*y(1) + x(1)*y(0) + x(2)*y(N-1) + … + x(N-1)*y(2) } + ω^(-m*2)*{ x(0)*y(2) + x(1)*y(1) + x(2)*y(0) + x(3)*y(N-1) + … + x(N-1)*y(3) } : + ω^(-m*(N-2))*{ x(0)*y(N-2) + x(1)*y(N-3) + … + x(N-1)*y(N-1) } + ω^(-m*(N-1))*{ x(0)*y(N-1) + x(1)*y(N-2) + … + x(N-1)*y(0) } N-1 N-1 = Σ ω^(-m*n)*{ Σ x(k)*y((n-k)mod N) } n=0 k=0 となる。最後の式に注目すると、これは N-1 z(n) = Σ x(k)*y((n-k)mod N) n=0,1,…,N-1 k=0 で定義された数列 z(n) のフーリエ変換に等しい。すなわち、数列 Z(m);m=0,1,…,N-1 は数列 z(n);n=0,1,…,N-1 のフーリエ変換である。 　ここで N を偶数に限定して N=2M とおき、次の条件を加える。 x(n) = y(n) = 0 n=M,M+1,…,2M-1 つまり、x(n) と y(n) の“後半”をすべて 0 にする。すると、z(n) は次のよ うになる。 z(0) = x(0)*y(0) z(1) = x(0)*y(1) + x(1)*y(0) : z(M-2) = x(0)*y(M-2) + x(1)*y(M-3) + … + x(M-2)*y(0) z(M-1) = x(0)*y(M-1) + x(1)*y(M-2) + … + x(M-2)*y(1) + x(M-1)*y(0) z(M) = x(1)*y(M-1) + … + x(M-2)*y(2) + x(M-1)*y(1) : 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 　これはまさに M 桁の多倍長乗算そのものである。言い換えると、B 進法で M 桁の数値 X と Y の B^n の位を数をそれぞれ x(n),y(n);n=0,1,…,M-1 とした とき、 M-1 X = Σ x(n)・B^n n=0 M-1 Y = Σ y(n)・B^n n=0 であり、X と Y の積は、 2M-1 X*Y = Σ z(n)・B^n n=0 である。 　まとめると、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 に比例する程度の時間で済む。 ━─────────────────────────────────── ● 2 組の実数列の同時 FFT 　2 組の実数列 x(n),y(n);n=0,1,…,N-1 が与えられたとき、それぞれの離散フ ーリエ変換 X(m),Y(m);m=0,1,…,N-1 を以下の手順で求めることができる。 　複素数列 t(n);n=0,1,…,N-1 を次のようにおく。 t(n) = x(n) + i・y(n) （x(n) と y(n) はいずれも実数なので、t(n) の実数部は x(n)、虚数 部は y(n) である） 　FFT を用いて t(n) の離散フーリエ変換 T(m);m=0,1,…,N-1 を求める。 　T(m) から、x(n) と y(n) の離散フーリエ変換 X(m),Y(m) を次の式で求める。 X(m) = {Re(T(m)) + Re(T(N-m))}/2 + i・{Im(T(m)) - Im(T(N-m))}/2 Y(m) = {Im(T(m)) + Im(T(N-m))}/2 - i・{Re(T(m)) - Re(T(N-m))}/2 （m=0 のとき、T(N-m) = T(N) = T(0)） 　この方法を用いると、2 組の実数列の離散フーリエ変換を、1 組の複素数列の 離散フーリエ変換とほぼ同じ時間で求めることができる[2]。 ━─────────────────────────────────── ● 2 組の実数列の同時 FFT を用いた多倍長乗算 　FFT を用いた多倍長乗算では、引数となる多倍長データを数列と見なしてそれ ぞれの離散フーリエ変換を行い、同じ番号の要素同士を掛け合わせて逆フーリエ 変換を行う。このように引数の数列 2 つに対して離散フーリエ変換を行うが、 それらの引数は実数列なので、2 組の実数列の同時 FFT を用いることができる。 　引数の数列を x(n),y(n);n=0,1,…,N-1 とおく。2 組の実数列の同時 FFT な ので、 t(n) = x(n) + i・y(n) とおき、FFT を用いて t(n) の離散フーリエ変換 T(m) を求める。 　以降、式が冗長になるので次のような文字で置き換えて表記する。 p := Re(T(m)) q := Im(T(m)) r := Re(T(N-m)) s := Im(T(N-m)) （?=Re(?)+i・Im(?)。Re(?) は ? の実数部、Im(?) は ? の虚数部） 　2 組の実数列の同時 FFT の関係から、X(m),Y(m);m=0,1,…,N-1 は次のように なる。 X(m) = {Re(T(m)) + Re(T(N-m))}/2 + i・{Im(T(m)) - Im(T(N-m))}/2 = (p + r)/2 + i・(q - s)/2 Y(m) = {Im(T(m)) + Im(T(N-m))}/2 - i・{Re(T(m)) - Re(T(N-m))}/2 = (q + s)/2 - i・(p - r)/2 　Z(m) = X(m)*Y(m) を p,q,r,s を使って展開すると、 Z(m) = X(m)*Y(m) = {(p + r)/2 + i・(q - s)/2}*{(q + s)/2 - i・(p - r)/2} = {(p + r) + i・(q - s)}*{(q + s) - i・(p - r)}/4 = {(p + r)*(q + s) - (q - s)*-(p - r)}/4 + i・{(p + r)*-(p - r) + (q - s)*(q + s)}/4 = (p*q + p*s + q*r + r*s + p*q - q*r - p*s + r*s)/4 + i・(- p^2 + r^2 + q^2 - s^2)/4 = (p*q + r*s)/2 + i・{(r^2 + q^2) - (p^2 + s^2)}/4 このように実数部が打ち消し合って簡単な式になる。表記を p,q,r,s から T(m) に戻すと、 Z(m) = (Re(T(m))*Im(T(m)) + Re(T(N-m))*Im(T(N-m)))/2 + i・{(Re(T(N-m))^2 + Im(T(m))^2) - (Re(T(m))^2 + Im(T(N-m))^2)}/4 となる。 　m=0 の場合について考えると、 Z(0) = Re(T(0))*Im(T(0)) である。また、N が偶数のとき、m=N/2 では、 Z(N/2) = Re(T(N/2))*Im(T(N/2)) である。 　この式を使うと、x(n),y(n) という 2 組の実数列の離散フーリエ変換を同時 に行うだけでなく、それぞれの離散フーリエ変換 X(m),Y(m) を算出せずに直接 Z(m) = X(m)*Y(m) を求めることになるので、FFT を用いた多倍長乗算の計算量 を減らすことができる。 ━─────────────────────────────────── ●他の多倍長演算のアプローチ 　詳しいことは次の機会に改めて書きたいが、多倍長除算について少しだけ触れ ておこう。 　何万桁というオーダーになったとき、多倍長の数値データ同士の除算をまとも にやろうとするのは間違いである。多倍長の数値データ同士の除算は乗算よりも 遥かに大きな計算量を必要とするので、なるべく除算をしなくて済むようにする べきだ。 　しかし、どうしても除算が必要なときは、「Newton 法の反復で除数の逆数を 求めてから被除数に掛ける」という手段を使う[1]。 　例えば、定数 A の逆数を求めるには、適切な方法で初期値 x(0) を決定し、 次の反復を行うとよい。 x(n+1) = 2*x(n) - A*x(n)^2 この反復は 1/A に 2 次の収束を示す（1 回の反復で有効桁数がほぼ 2 倍にな る）。 　多倍精度実数の平方根を求める場合にも Newton 法が有効だ。 　当然、Newton 法の過程で必要になる多倍長の整数同士の乗算には FFT を使う。 ━─────────────────────────────────── ●おわりに 　「X68k で円周率を計算しました」なんて言うと、理論をよく理解もせずにた だ速いだけのパソコンを使っている……というよりパソコンに使われているよう な人に「そんな遅いパソコンで円周率を計算して何が楽しいのか」なんて言われ て頭にくることがあるので、もしも円周率の計算に関して書くとしたら、プログ ラムを示すことよりも、可能な限り自分で理解した上で解説文を書くことに重点 を置きたいと思う。 　それにしても AGM（算術幾何平均）の理論[7]は難しいのぉ。私にゃいまだに よくわからん。 ━─────────────────────────────────── ●参考文献 [1] D.E.Knuth, 中川圭介訳, KNUTH The Art of Computer Programming 第4分冊 準数値算法／算術演算, サイエンス社, ISBN4-7819-0426-2 [2] 五十嵐善英, 整数の乗算法, 情報処理 1983 Apr. Vol.24 No.4 pp.358-362 [3] 小池慎一, Ｃによる科学技術計算, ＣＱ出版社, ISBN4-7898-3032-2 [4] 奥村晴彦, Ｃ言語による最新アルゴリズム辞典, 技術評論社, ISBN4-87408-414-1 [5] H.P.スウ, 佐藤平八訳, 工学基礎演習シリーズ1 フーリエ解析, 森北出版, ISBN4-627-93010-0 [6] 丹明彦, そこにπがあるから, Oh!X 1988.8 pp.51-60 [7] JONATHAN M. BORWEIN/PETER B. BORWEIN, Pi and the AGM, A Wiley-Interscience Publication, ISBN0-471-83138-7 （ＥＯＦ）