ここはミラーサイトです。本物は http://homepage2.nifty.com/m_kamada/math/fftmul.htm
counterSince June 16, 2000STUDIO KAMADAJapanese to English by @nifty
戻る | サイトマップ | ホーム
離散フーリエ変換を用いた多倍長乗算の話2007-09-22(Sat) 02:16

はじめに

　CPU のレジスタに入り切らないような巨大あるいは高精度な数値データの演算を行うとき、四則演算をやるにもそのためのルーチンを用意しなければならない。(*1)
　ここでは、「CPU のレジスタではちょっとはみ出す」という程度の数値データのことは考えない。何万桁、何百万桁という普通では目にしないような巨大な数値データを扱う場合に掛け算をどうすればよいかについて考える。(*2)
　はっきり言ってそんな巨大な数が扱えたところで日常生活にこれといった恩恵があるわけではないが、敢えて常識を覆すことによって掛け算の本質が見えてくるかも知れない……などと論文調にこじつけてみたけれど、実は単なる好奇心だったりする。ただ、常識が覆されるのは事実だ。私は実際に巨大な数の計算を必要としているのではなくて、常識を覆えすことに快感を覚えるのである。(*3)
　先に断っておくと、以下の解説を全部理解するにはそれなりの数学の知識が必要だと思われる。数式も出てくるので、数学が苦手な人は注意されたし。細かい定理や法則までは解説しないので、わからない点は各自で調べていただきたい。私だって、こんなことを学校で習ったわけではない。好奇心からその手の本を読み漁って覚えた程度なので、専門に勉強している人からのご指摘はお手柔らかに願いたい。(*4)


(*1) そのような多倍長の数値データの四則演算を行うライブラリとしては、電脳倶楽部でも Vol.26 と Vol.28 に LONG.C というプログラムが掲載されている。が、 LONG.C では今回解説するような高度なアルゴリズムは使われていない。

(*2) こんなの使うの、円周率の計算くらいだよね。

(*3) 勿論、アルゴリズムの検証のために実際に計算してみることは重要だ。

(*4) FFT で掛け算をやろうという時点で既に FFT の使い方を間違っているという話もある。

▲

目次

　　　１． 筆算の方法を用いた多倍長乗算
　　　２． 分割統治法を用いた多倍長乗算
　　　３． 離散フーリエ変換
　　　４． FFT を用いた多倍長乗算
　　　５． 2 組の実数列の同時 FFT
　　　６． 2 組の実数列の同時 FFT を用いた多倍長乗算
　　　７． 他の多倍長演算のアプローチ
　　　８． おわりに
　　　９． 参考文献
　　１０． 関連リンク

▲

１． 筆算の方法を用いた多倍長乗算

　筆算の方法で N 桁の整数同士の乗算を行うには、だいたい N2 に比例する時間が必要である。
　B 進法（B は 2 以上の任意の整数）で表現したとき N 桁になる整数 X および Y を、それぞれ次のように表現する。

　　　　 X =
N-1
Σ x(n)・Bn
n=0

Y =
N-1
Σ y(n)・Bn
n=0

ここで、 x(n),y(n);n=0,1,…,N-1 は B 進法で 1 桁の数とする。
　X と Y の積を筆算の方法で求めると、次のようになる。

　　　　 X*Y =
N-1 ｛ N-1 ｝
Σ Σ x(k)・Bk *y(n)・Bn
n=0 k=0

=
N-1 N-1
Σ Σ x(k)*y(n)・Bk+n
n=0 k=0

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

▲

２． 分割統治法を用いた多倍長乗算

　分割統治法を用いると、N 桁の整数同士の乗算を Nlog2(3)≒N1.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・B2 + （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・B2 + （x1*y0 + x0*y1）・B + x0*y0

=
x1*y1 ・B2 + ｛ 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 = 2k と仮定したとき、B 進法で N 桁の整数同士の乗算を分割統治法で行うと、B 進法で 1 桁の整数同士の乗算が 3k 回必要となる。 k = log2(N) より、この方法での B 進法で N 桁の整数同士の乗算のコストは、 3k = 3log2(N) = Nlog2(3) ≒ N1.585 程度で済むことになる。実際には N が 2k で表現できるとは限らないので、 2k で表現できないときは N1.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 とする。すなわち、

　　　　 X(m) =
N-1
Σ ω-m*n*x(n)　　　　　m=0,1,…,N-1
n=0

　　　　 （一般的に離散フーリエ変換には係数 1/N が付くが、ここでは省略する）

　ここで、ω は 1 の N 乗根の 1 つであり、次のように定義される。

　　　　 ω = e2π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 をすべて求めるのに乗算を N2 回必要とするが、これを 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 とする。すなわち、

　　　　 X(m) =
N-1
Σ ω-m*n*x(n)　　　　　m=0,1,…,N-1
n=0

Y(m) =
N-1
Σ ω-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

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

　　　　 z(n) =
N-1
Σ 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 の Bn の位を数をそれぞれ x(n),y(n);n=0,1,…,M-1 としたとき、

　　　　 X =
M-1
Σ x(n)・Bn
n=0

Y =
M-1
Σ y(n)・Bn
n=0

であり、X と Y の積は、

　　　　 X*Y =
2M-1
Σ z(n)・Bn
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・B2 までの数だが、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・（- p2 + r2 + q2 - s2）/4

= （p*q + r*s）/2
+ i・｛（r2 + q2） - （p2 + s2）｝/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 を用いた多倍長乗算の計算量を減らすことができる。



補足
　2組の実数列の同時FFTを用いる方法は式変形が面白いので書きましたが、この方法を用いるよりも実数列のFFT（長さが半分の複素数列に置き換える）を2回使う方法のほうが速く計算できます。


▲

７． 他の多倍長演算のアプローチ

　詳しいことは次の機会に改めて書きたいが、多倍長除算について少しだけ触れておこう。
　何万桁というオーダーになったとき、多倍長の数値データ同士の除算をまともにやろうとするのは間違いである。多倍長の数値データ同士の除算は乗算よりも遥かに大きな計算量を必要とするので、なるべく除算をしなくて済むようにするべきだ。
　しかし、どうしても除算が必要なときは、「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] 小池慎一, Cによる科学技術計算, CQ出版社, ISBN4-7898-3032-2

[4] 奥村晴彦, C言語による最新アルゴリズム辞典, 技術評論社, 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

▲

１０． 関連リンク

・Ooura's Mathematical Software Packages
　大浦拓哉氏のページです。高速で信頼性も高いことで定評がある「汎用 FFT パッケージ」、FFT の定義から高速化のテクニックまでを網羅した「FFT (高速フーリエ・コサイン・サイン変換) の概略と設計法」、個人のパソコンで計算できる円周率計算プログラムとしてはこれ以上速いものはないのではないかと思われる「FFT と AGM による円周率計算プログラム」などを公開されています。

・SN library
　津留和生氏のページです。「アセンブラを使わないＣ＋＋による多倍長計算ライブラリ」を公開されています。超越関数まで網羅した充実した多倍長計算ライブラリです。

▲

更新履歴
(1999-12-10) 月刊電脳倶楽部128号（1999年1月号）の読み物横町のコーナーに収録
(2000-08-22) HTML化、加筆修正
(2000-10-10) HTMLを更新
(2001-04-19) HTMLを更新、関連リンクを追加
(2002-11-10) URLを変更
(2007-09-22) 関連リンクを更新
▲
戻る | サイトマップ | ホーム
Copyright (C) 1999-2007 Makoto Kamada All Rights Reserved.Return from mirror
[PR] ｜ インプラント｜ 債務整理　相談｜ 債務整理｜ インプラント｜ 転職サイト｜ SEO｜ アクセス解析｜ ハウスメーカー｜ レンタルオフィス｜ SEO対策｜ 消費者金融｜ 不動産担保ローン｜ 時計｜ 車 買取｜ ハワイ挙式｜ アスクル｜ 転職｜ 生命保険｜ テンプレート｜ 沖縄旅行｜ 動画｜ 免許合宿｜ 二輪｜ 引越し｜ 消費者金融｜ 税理士｜ ゴルフ会員権｜ 留学｜ レーシック｜ マッサージ｜ FX｜ 投資信託｜ くりっく365｜ アフィリエイト｜ 育毛剤｜ FX｜ ホームページ制作｜ デイトレード｜ FX｜ ホノルルマラソン｜ ベスト ハワイ ホテル レーツ｜ バリ島｜ ハワイウエディング｜ Hawaii hotels｜ Hawaii Activities｜ bhhr｜
【運営会社「パラダイムシフト」サービス】 ハワイ現地オプショナルツアー（リラックマ.ハワイホテルガイド） - ビジネスクラス航空券 - 格安航空券(1) - 格安航空券(2) - 海外ホテル - 韓国旅行
無料ホームページ作成 - レンタルサーバー - 携帯ホームページ - ブログ - ホテル 予約 - タイムシェア - ヴィラ - ハワイ　コンドミニアム - バリ島　ホテル - ハワイ　不動産 - プーケット　ホテル
[PR] 1番高くバイクを売れる業者を探しませんか？厳選サイトをご紹介