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

高速フーリエ変換
(fast Fourier transform)
1
• 高速フーリエ変換とは？
– 簡単に言うとフーリエ変換を効率よく計算する
方法
– アルゴリズムの設計技法は分割統治法に基
づいている
• 今回の目的は？
– 多項式の積を求める問題を取り上げ、高速
フーリエ変換のアルゴリズムを用いた解法を
導いて、計算時間を驚くほど短縮できることを
示す
2
多項式の積の計算の例１
• 以下のような2次多項式のｐ（ｘ）とｑ（ｘ）をかける
p( x )  1  x  x 2
q( x )  2  x  x 2
p( x )q( x )  (1  x  x 2 )(2  x  x 2 )
 2  x  x 2  2x  x 2  x 3  2x 2  x 3  x 4
 2  x  2x 2  x 4
この場合２次式のため簡単に出来たが、一般的
にN-1次多項式の時は、このような方法では大変
時間 N2回の乗算 がかかる.どうすればよいのか？
3
多項式の積の計算を高速化するに
は？
• N-1次の多項式はN個の異なる点における
関数値で完全に決まるという事実を利用
する（N-1次式の項の数は定数も含むから
N個）
⇒２つのN-1次多項式同士の積は２N-1個の点
でその多項式の値が分かっていれば、その積
を完全に決定できる
4
２つのN-1次多項式の積を求める
基本的な手順
• 評価：入力された多項式の値を２N-1個の
異なる点で計算する
• 乗算：各点でえられた値を掛ける
• 補間：上で求めた２N-1個の値を補間する
– 補間とは、多項式のN個の異なる点における
値が与えられた時、その多項式の係数を求め
ること
5
２つのN-1次多項式
P(x)
Q(x)
評
価
乗算
N2の計算時間
フーリエ変換
分割統治法を用いて
効率よくしたのが高速フーリエ変換
p(a1 ),, pa 2 N1 
q(a1 ),, qa 2 N1 
乗算
2N-1次の多項式
R(x)=P(x)Q(x)
補
間
逆フーリエ変換
r ( x )  p ( x )q ( x )
r (a1 ),, r (a 2 N1 )
6
• 例２）
p( x )  1  x  x 2
q( x )  2  x  x 2
r x   p( x )q( x )を求めたい
x  2,1,0,1,2においてpx , qx を計算すると
p(2), p(1), p(0), p(1),p(2)  3,1,1,3,7
q(2), q(1), q(0), q(1),q(2)  8,4,2,2,4
となる．すると
r ( x )は
r(2), r(1), r(0), r(1),r(2)  24,4,2,6,28
7
ラグランジュの公式を用いると
x 1 x  0 x 1 x  2
r ( x )  24
 2 1  2  0  2 1  2  2
x  2 x  0 x 1 x  2
4
1 2 1  0 11 1  2
x  2 x 1 x 1 x  2
2
0  2 0 1 0 1 0  2
x  2 x 1 x 1 x  2
6
1 2 11 1 0 1 2
x  2 x 1 x  0 x 1
 28
2  2 2 1 2  0 2 1
   2  x  2x 2  x 4
となる．
8
ラグランジュの公式
N  1次多項式 P (x)のN個の点 x1 , x 2 ,, x Nと
そこでの値y1 , y 2 ,, y Nが与えられた時、
P( x1 )  y1 , P( x 2 )  y 2 , , P( x N )  y N
を満たす唯一の N  1次多項式を求めたい．（補間する）
この時の古典的な解はラグランジュの補間公式を用いること
x  xi
P( x )   y j 
1 j N 1i  N x j  x i
（ラグランジュの補間
公式）
i j
9
2
N
• 例２の場合、評価と補間それぞれでの計
算時間がかかってしまう
• しかし、多項式の値を計算する２N-1個の
点はどこの点を選んでもよいことに注目す
ると、多項式の積の計算を簡単にする点
集合を選べば、効率のよく計算時間を短縮
できるのでは？
10
１の複素累乗根
• 多項式の積の計算をする点として、１の複
素累乗根を用いると計算を高速化できる
２つの複素数の積は
(a  bi )(c  di )  (ac  bd )  (ad  bc )i
となる．次の例のよう
に乗算の結果、実数部や虚数部が０になることもある．
(1  i)(1  i)  2i
(1  i)8  16
8
上の式を16  2 で割ると
8
i 
 1


 1
2
 2
をえる．
11
一般に、その数を累乗して１となる複素数は
たくさんあり、
１の累乗根とよばれる．
実際、自然数Nに対して、
Z N  1となる複素数ZがN個ある．
これらの中の１つは「
１の原始N乗根」とよばれ、
それをWNと表わすと、
全ての根は k  0,1,2, , N  1に対して、
WN をk乗してえられる．
例）１の８乗根
W80 , W81 , W82 , W83 , W84 , W85 , W86 , W87
W80  1で、W81が原始８乗根である　.
Nが偶数の時は WN 2は  1であるから、
W84  1となる．
N
12
N-1次多項式の値を
１のN乗根において計算する
• 1のN乗根における値をそのまま計算した
のではあまり意味がない
• 分割統治法を用いる
– 低次の項から順に交互にとって２つに分ける
– N-１次多項式の値をN個の点で計算するため
に、まずN/2個の係数の２つの多項式に分割
する。そしてそれぞれの多項式の値をN/2個
の点で計算して全体のN点における値を求め
ればよい
13
• なぜ分割統治法がいいのか？
– 根を２乗すると別の根が得られる
– 特にNが偶数の時は、
１のN乗根を２乗すると、１のN/2乗根が得ら
れる
例） ( W 1 ) 2  W 1
8
4
( W41 ) 2  W21
実際にN=8で考えてみる
14
例３）N  8の時
p( x )  p 0  p1x  p 2 x 2  p 3 x 3  p 4 x 4  p 5 x 5  p 6 x 6  p 7 x 7
 (p 0  p 2 x 2  p 4 x 4  p 6 x 6 )  x (p1x  p 3 x 2  p 5 x 4  p 7 x 6 )
 p e ( x 2 )  xpo ( x 2 )
上のように２つに分割する．ここで、７次多
項式 p( x )を１の８乗根で計算する
W8 : w 80 , w 18 , w 82 , w 83 , w 84 , w 85 , w 86 , w 87
まず、w 84  1だから
W8 : w 80 , w 18 , w 82 , w 83 , w 80 , w 18 , w 82 , w 83
となる．次に各項を２
乗すると
W82 : w 04 , w14 , w 24 , w 34 , w 04 , w 14 , w 24 , w 34
になる．
15
p( x )  p e ( x 2 )  xpo ( x 2 )
に注目すると以下のような次式が導かれる．
p( w 80 )  p e ( w 04 )  w 80 p o ( w 04 ),
p( w 18 )  p e ( w 14 )  w 18 p o ( w 14 ),
p( w 82 )  p e ( w 24 )  w 82 p o ( w 24 ),
p( w 83 )  p e ( w 34 )  w 83 p o ( w 34 ),
p( w 84 )  p e ( w 04 )  w 80 p o ( w 04 ),
p( w 85 )  p e ( w 14 )  w 18 p o ( w 14 ),
p( w 86 )  p e ( w 24 )  w 82 p o ( w 24 ),
p( w 87 )  p e ( w 34 )  w 83 p o ( w 34 ),
16
実際に、p( w18 )  p e ( w14 )  w 18 p o ( w14 )を計算してみる．
p( w18 )  p e ( w14 )  w18 p o ( w14 )
 (p 0  p 2 w14  p 4 w 24  p 6 w 34 )  w18 (p1  p 3 w14  p 5 w 24  p 7 w 34 )
 (p 0  p 4 w 24 )  w14 (p 2  p 6 w 24 )  w18{( p1  p 5 w 24 )  w14 (p 3  p 7 w 24 )}
 (p 0  p 4 w12 )  w14 (p 2  p 6 w12 )  w18{( p1  p 5 w12 )  w14 (p 3  p 7 w12 )}
ここで、
w12  1だから
 (p 0  p 4 )  w14 (p 2  p 6 )  w18{( p1  p 5 )  w14 (p 3  p 7 )}
またw  ( w )
1
w  (w )
1
1
4
1
8
1
2
1
4
2
2
 (1)
 (i)
1
2
1
2
 i,
 i だから
 (p 0  p 4 )  i(p 2  p 6 )  i{( p1  p 5 )  i(p 3  p 7 )}
17
一般的に、１のN乗根におけるp( x )の計算をする
p e x とp o x とを1のN / 2乗根において再帰的に計算する．
このような乗算を
N回繰り返す．
（ただし、再帰的な計
算でN / 2, N / 4,となるので、
Nは2のベキ乗であると仮定する．）
再帰は w 2xのようにN  2の時に打ち切る．
このときxはw12  1, ( w12 ) 2  1 のどちらかだから
p 0  p1xはp 0  p1かp 0  p1を得る．
あとは1つずつ上にあげて計算し、全体のp( x )を求める．
18
性質４１．１
約N lg N回の乗算によって、N  1次多項式の値を１のN乗根において計算できる
証明）第６章漸化式４
を参照
乗算の回数を M ( N)とする．分割統治法の再帰式は
M ( N )  2 M ( N / 2)  N
これを解くと、
N  2 n とおく
M (2 n )  2M (2 n 1 )  2 n
M (2 n ) M (2 n 1 )

1
n
n 1
2
2
M(2 n  2 )

11
n 2
2
M ( 2 n 3 )

111
n 3
2
n
よって ( N  2 n よりn  lg N)
M ( N )  M (2 n )  n 2 n  N lg N
19
１の累乗根の関数値の補間
• １の累乗根を用いて高速で多項式の値が得られた
⇒その値を高速で補間すれば、高速で多項式の積が
求まったことになる．
補間する点集合が１の累乗根ならば、
関数値を計算する手順で補間多項式も求まる！
なぜなら
フーリエ変換とその逆変換を使っているから！！
20
前と同じN  8の例で、補間は、多項式
r ( x )  r0  r1x  r2 x 2  r3 x 3  r4 x 4  r5 x 5  r6 x 6  r7 x 7
が１の８乗根において次の値をとるように係
数ri (0  i  8)を求めることである．
r ( w 80 )  s 0, r ( w 18 )  s1, r ( w 82 )  s 2, r ( w 83 )  s 3,
r ( w 84 )  s 4, r ( w 85 )  s 5, r ( w 86 )  s 6, r ( w 87 )  s 7 .
各点における値を係数とする多項式を
s( x )  s 0  s1x  s 2 x 2  s 3 x 3  s 4 x 4  s 5 x 5  s 6 x 6  s 7 x 7
とする．多項式s( x )の値を１の累乗根の逆数
W81 : w 80 , w 81 , w 8 2 , w 83 , w 8 4 , w 85 , w 86 , w 87 .
において計算すると、r ( x )の係数
r0 , r1 , r2 , r3 , r4 , r5 , r6 , r7 .
が求まる．（
r ( w 18 )  s1だからs( w 81 )  8r1）
ここで、数列
W81は１の８乗根をW8の逆順にしただけ
W81 : w 80 , w 87 , w 86 , w 85 , w 84 , w 83 , w 82 , w 18 .
つまり値の計算すべき点を並び替えるだけで
関数値の計算の手順がそのまま補間の計算に使える．
21
1のN乗根の t番目の値の逆数におけるs( x )の値を計算すると、
s( w Nt ) 

0 j N
 r(w
0 j N

t j
s
(
w
 j N)
j
N
)(w Nt ) j

j i
t j
r
(
w
)
(
w
 i N N)

j i  t 
r
w
i N
 ri
j i  t 
w
 N  Nrt
0 j N 0i  N

0 j N 0i  N

0i  N 0 j N
となる．
i  tの時、 ri  rt ,
0i  N
w 
0 j N
j it 
N
N
i  t  N
w
1
i  tの時、 w Nji  t   Ni  t 
0
w N 1
0 j N
22
1の8乗根の1番目の値の逆数におけるs( x )の値を計算すると、
s( w 81 ) 
1 j
s
(
w
 j 8)
0 j8
 s 0  s1w 81  s 2 ( w 81 ) 2    s 7 ( w 81 ) 7

j
1 j
r
(
w
)(
w
 8 8)
0 j8
   w  rw w 
r w r w


0
8
1
8
1
8
j i
1 j
r
(
w
)
(
w
i 8 8)
0 j8 0i 8

2
8
 
2 2
8
 w 
 r w
7
8
7 7
8

 r0  r1w 80    r7 ( w 80 ) 7  r0  r1w 18    r7 ( w 18 ) 7 w 81


   r0  r1w 87    r7 ( w 87 ) 7 ( w 81 ) 7

ji 1
r
w
i 8

0 j8 0i 8

r w 
i
0 j8 0i 8
j i 1
8
 8r1
23
• 性質４１．２
１のN乗根における値が与えられたN-1次多項式の係数
を、約NlgN回の乗算で求めることが出来る．
これは補間の計算は多項式の値の計算の逆であるから
性質４１．１より乗算は約NlgN回となる．
多項式の値の計算と補間の計算の乗算の回数は
NlgN回であるが、あくまで１のN乗根に対してだけ
あてはまる
24
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次より上位の係数は０に
初期化されている．
25
手続きeval
• １番目の引き数
– 多項式の係数を、１の累乗根における多項式
の値に書き換える
• ２番目の引き数
– 多項式の次数（係数の数、１の根の数より１少
ない）
手続きevalの目的
N+1個の係数が連続して入っている配列を入力と
して、同じ配列にN+1個の関数値を入れて返すこと
26
• しかし、再帰呼出しでは連続していない偶数番目
と奇数番目の係数の配列を入れかえなければな
らない．
⇒ここで用いられるのが、完全シャッフル
• 完全シャッフルをすることにより
配列の前半に連続して偶数番目の係数を入れ、
後半に奇数番目の係数を入れることが出来る．
27
完全シャッフル
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]
28
高速フーリエ変換のプログラム
• 実際計算する時の1のN乗根の値は、
 2j 
 2j 
w  cos
  i sin

 N 1
 N 1
j
N
• 配列wに1の(outN+1)乗根が入っている
上のことも考慮すると
次のような高速フーリエ変換のプログラム
となる．
29
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];
}
30
}
• 性質４１．３
2つのN-1次多項式の積は2NlgN+O(N)回
の複素数の計算によって求められる．
性質４１．１と性質４１．２より、N-1次多項
式を1のN乗根の値への変換と補間の乗
算の回数はそれぞれNlgN回である．また1
のN乗根に変換した値を掛け合わせるの
にO(N)の乗算が必要だから、積の回数は
2NlgN+O(N)回となる．
31
まとめ
◎2つのN  1次多項式の積を求める時、
普通に掛けると、乗算の回数は N 2回
高速フーリエ変換を使うと、
2 N lg N  ON 回
◎高速フーリエ変換の手順
①：
2つのN  1次多項式を 2 N  1個の異なる点で計算する．
（このとき
1のN乗根を使う）
②：変換した値を掛け合わせる．
③：①と同じ計算手順により補間する
32

Download Report