高速フーリエ変換 (fast Fourier transform)

高速フーリエ変換(fast Fourier transform)

高速フーリエ変換とは？ • 簡単に言うとフーリエ変換を効率よく計算する方法 • アルゴリズムの設計技法は分割統治法に基づいている • 今回の目的は？ • 多項式の積を求める問題を取り上げ、高速フーリエ変換のアルゴリズムを用いた解法を導いて、計算時間を驚くほど短縮できることを示す

多項式の積の計算の例１ • 以下のような2次多項式のｐ（ｘ）とｑ（ｘ）をかける　この場合２次式のため簡単に出来たが、一般的にN-1次多項式の時は、このような方法では大変時間　　　　　　がかかる.どうすればよいのか？

多項式の積の計算を高速化するには？ • N-1次の多項式はN個の異なる点における関数値で完全に決まるという事実を利用する（N-1次式の項の数は定数も含むからN個） ⇒２つのN-1次多項式同士の積は２N-1個の点でその多項式の値が分かっていれば、その積を完全に決定できる

２つのN-1次多項式の積を求める　　　　　　　　　　　　　基本的な手順２つのN-1次多項式の積を求める　　　　　　　　　　　　　基本的な手順 • 評価：入力された多項式の値を２N-1個の　　　　異なる点で計算する • 乗算：各点でえられた値を掛ける • 補間：上で求めた２N-1個の値を補間する • 補間とは、多項式のN個の異なる点における値が与えられた時、その多項式の係数を求めること

２つのN-1次多項式 P(x)Q(x) 2N-1次の多項式 R(x)=P(x)Q(x) 乗算評価補間フーリエ変換分割統治法を用いて効率よくしたのが高速フーリエ変換逆フーリエ変換乗算

例２）

ラグランジュの公式

例２の場合、評価と補間それぞれで　の計算時間がかかってしまう • しかし、多項式の値を計算する２N-1個の点はどこの点を選んでもよいことに注目すると、多項式の積の計算を簡単にする点集合を選べば、効率のよく計算時間を短縮できるのでは？

１の複素累乗根 • 多項式の積の計算をする点として、１の複素累乗根を用いると計算を高速化できる

N-1次多項式の値を　　　　　　　　１のN乗根において計算するN-1次多項式の値を　　　　　　　　１のN乗根において計算する • 1のN乗根における値をそのまま計算したのではあまり意味がない • 分割統治法を用いる • 低次の項から順に交互にとって２つに分ける • N-１次多項式の値をN個の点で計算するために、まずN/2個の係数の２つの多項式に分割する。そしてそれぞれの多項式の値をN/2個の点で計算して全体のN点における値を求めればよい

なぜ分割統治法がいいのか？ • 根を２乗すると別の根が得られる • 特にNが偶数の時は、　１のN乗根を２乗すると、１のN/2乗根が得られる例）実際にN=8で考えてみる

１の累乗根の関数値の補間 • １の累乗根を用いて高速で多項式の値が得られた ⇒その値を高速で補間すれば、高速で多項式の積が求まったことになる．補間する点集合が１の累乗根ならば、関数値を計算する手順で補間多項式も求まる！フーリエ変換とその逆変換を使っているから！！なぜなら

性質４１．２１のN乗根における値が与えられたN-1次多項式の係数を、約NlgN回の乗算で求めることが出来る．性質４１．２１のN乗根における値が与えられたN-1次多項式の係数を、約NlgN回の乗算で求めることが出来る．　これは補間の計算は多項式の値の計算の逆であるから性質４１．１より乗算は約NlgN回となる．　多項式の値の計算と補間の計算の乗算の回数は NlgN回であるが、あくまで１のN乗根に対してだけ　あてはまる

eval(p, outN, 0); eval(q, outN, 0); for (i = 0; i <= outN; i++) r[i] = p[i]*q[i]; eval(r, outN, 0); for (i = 1; i <= N; i++) { t = r[i]; r[i] = r[outN+1-i]; r[outN+1-i] = t; } for (i = 0; i <= outN; i++) r[i] = r[i]/(outN+1); ・大域変数outNに２N-1が与えられている．・p,q,rは添字の範囲が０から２N-1の複素数の配列と仮定されている．・多項式p,qはN-1次であり、N次より上位の係数は０に初期化されている．

手続きeval • １番目の引き数 • 多項式の係数を、１の累乗根における多項式の値に書き換える • ２番目の引き数 • 多項式の次数（係数の数、１の根の数より１少ない）手続きevalの目的 N+1個の係数が連続して入っている配列を入力として、同じ配列にN+1個の関数値を入れて返すこと

しかし、再帰呼出しでは連続していない偶数番目と奇数番目の係数の配列を入れかえなければならない．⇒ここで用いられるのが、完全シャッフルしかし、再帰呼出しでは連続していない偶数番目と奇数番目の係数の配列を入れかえなければならない．⇒ここで用いられるのが、完全シャッフル • 完全シャッフルをすることにより配列の前半に連続して偶数番目の係数を入れ、後半に奇数番目の係数を入れることが出来る．

完全シャッフル P[0] P[1] P[2] P[3] P[4] P[5] P[6] P[7] P[0] P[2] P[4] P[6] P[1] P[3] P[5] P[7]

高速フーリエ変換のプログラム • 実際計算する時の1のN乗根の値は、 • 配列wに1の(outN+1)乗根が入っている　上のことも考慮すると次のような高速フーリエ変換のプログラムとなる．

eval(struct complex p[], int N, int k) { int i, j; if (N == 1) { p0 = p[k]; p1 = p[k+1]; p[k] = p0+p1; p[k+1] = p0-p1; } else { for (i = 0; i <= N/2; i++) {j = k+2*i; t[i] = p[j]; t[i+1+N/2] = p[j+1]; } for (i = 0; i <= N; i++) p[k+i] = t[i]; eval(p, N/2, k); eval(q, N/2, (k+1+N)/2); j = (outN+1)/(N+1); for (i = 0; i <= N/2; i++) { p0 = w[i*j]*p[k+(N/2)+1+i]; t[i] = p[k+i]+p0; t[i+(N/2)+1] = p[k+i]-p0; } for (i = 0; i <= N; i++) p[k+i] = t[i]; } }

性質４１．３2つのN-1次多項式の積は2NlgN+O(N)回の複素数の計算によって求められる．性質４１．３2つのN-1次多項式の積は2NlgN+O(N)回の複素数の計算によって求められる．　性質４１．１と性質４１．２より、N-1次多項式を1のN乗根の値への変換と補間の乗算の回数はそれぞれNlgN回である．また1のN乗根に変換した値を掛け合わせるのにO(N)の乗算が必要だから、積の回数は2NlgN+O(N)回となる．

高速フーリエ変換 (fast Fourier transform)