高速フーリエ変換 (FFT)

再帰型の FFT はやっぱり遅いので、非再帰型の FFT について調べた。『アルゴリズムイントロダクション 第 3 巻』の「多項式FFT」を参考にしている。多項式の次数が n といったら、この記事の中では $a_0+a_1x + \cdots + a_{n-1} x^{n-1}$ で表せる多項式のことだと思ってください。

多項式FFT

多項式を表現するとき、係数を持つ方法と点での値を持つ方法の 2 通りが考えられる。 点表現というのは $f(x_0),f(x_1),\ldots,f(x_{n-1})$ という $n$ 個の具体値で多項式を表現する方法である。 $(fg)(x_i)=f(x_i)g(x_i)$ なので点表現は乗算が $O(n)$ でできるという利点がある。

係数表現で乗算をするとき、係数表現から点表現への高速な変換・逆変換があれば、点表現に移す→乗算する→係数表現に戻すというプロセスによって高速に多項式乗算ができる。 実際に、多項式 $f$ の次数を $n$、$w=e^{2\pi/ni}$ としたとき、$f(w^0),f(w^1),f(w^2),\ldots,f(w^{n-1})$ が高速に求められる。 これを高速フーリエ変換 (FFT) と呼ぶ。

再帰FFT

多項式の次数は 2 の冪であるものとする。$w_n=e^{2\pi/ni}$ とする。$w_n^{i}$ での値を求めることを考えよう。

$$f(w_n^i)=a_0+a_1w_n^i+a_2w_n^{2i}+\cdots+a_{n-1}w_n^{(n-1)i}$$

偶数次と奇数次で分離する。

$$f(w_n^i)=(a_0+a_2w_n^{2i}+\cdots+a_{n-2}w_n^{(n-2)i}) +w_n^i(a_1+a_3w_n^{2i}+\cdots+a_{n-1}w_n^{(n-2)i})$$

$w_n^{2i}=w_{n/2}^{i}$であることを考えると次のように変形できる。

$$f(w_n^i)=(a_0+a_2w_{n/2}^{i}+\cdots+a_{n-2}w_{n/2}^{(n-2)/2i}) +w_n^i(a_1+a_3w_{n/2}^{i}+\cdots+a_{n-1}w_{n/2}^{(n-2)/2i})$$

これを見ると、多項式

$$ f_0(x)=a_0+a_2x+\cdots+a_{n-2}x^{(n-2)/2i}$$ $$ f_1(x)=a_1+a_3x+\cdots+a_{n-1}x^{(n-2)/2i}$$

について、単位円を $n/2$ 等分する点での値さえ分かれば計算できることが分かる。具体的には以下のようになる。

$$ f(w_n^i)=f_0(w_{n/2}^i)+w_n^{i}f_1(w_{n/2}^i) $$

$f_0,f_1$ は $n/2$ 次の多項式である。次数 $n$ の多項式FFT するのに掛かる時間を $T(n)$ とすると、$T(n)=2T(n/2)+O(n)=O(n \log n)$ となる。

再帰型を考える前に、もう少しだけこの式を変形する。$w_n^{i+n/2}=-w_n^{i}$ を使う。

$$ f(w_n^i)=f_0(w_{n/2}^i)+w_n^{i}f_1(w_{n/2}^i) $$ $$ f(w_n^{i+n/2})=f_0(w_{n/2}^i)-w_n^{i}f_1(w_{n/2}^i) $$

再帰FFT

再帰の様子を以下の図に示した。

f:id:pekempey:20161024162527p:plain

再帰で書くために、まず次のような置換をする。

$$(0,1,2,3,4,5,6,7)\to(0,4,2,6,1,5,3,7)$$ 何が起きてるのか分かりづらいが、二進表現で見てみると $$(000,001,010,011,100,101,110,111)\to(000,100,010,110,001,101,011,111)$$ となり、ビットを逆転させた位置に移動させているだけである。下位ビットの0/1で振り分けていることを考えれば何故こうなるのかは分かるだろう。

さて 4 段目から 3 段目への変換を考えてみよう。

f:id:pekempey:20161024163356p:plain

これは次のような計算になる。+は加算器、×は乗算器を表す。

f:id:pekempey:20161024165311p:plain

3段目から2段目への変換は下図のようになる。

f:id:pekempey:20161024170119p:plain

最後 2 段目から 1 段目への変換をくっつければ完成である。

f:id:pekempey:20161024171154p:plain

これはよく見る図。実際は一番最初にビット逆転をする。

おまけ

Number Theoretic Transform という mod 上の FFT みたいなのもある。サイズ n の NTT をするとき mod の原始 n 乗根が必要になるので、原始 220 乗根くらいがある mod でないと使えない。ただし畳み込みをするだけなら任意の mod でも可能。

おまけ(2)

1,w,w2,w3,... という列を得るには 1 に w を何回も掛けていけば求まるには求まるんだけど、誤差が乗るので注意。最終的な値が int に収まるなら多分大丈夫だけど、codechef の問題で誤差死したことがあったので一応。

$ cat a.cpp
#include <iostream>
#include <complex>
#include <cmath>

using namespace std;

const double pi = acos(-1);

int main() {
  int n;
  scanf("%d", &n);
  complex<double> w = polar(1.0, 2 * pi / n);
  complex<double> ws = 1;
  for (int i = 0; i < n; i++) {
    ws *= w;
  }
  printf("%.16f %.16f\n", ws.real(), ws.imag());
}
$ echo "100000" | ./a.out
0.9999999999964125 0.0000000000000005
$

おまけ(3)

再帰の節の数式がごちゃごちゃしてると思うので、もう少し見やすい(?)ものを書いておく。

let n = 8
let f(x)  = a[1]x^0 + a[1]x^1 + a[2]x^2 + a[3]x^3 + a[4]x^4 + a[5]x^5 + a[6]x^6 + a[7]x^7
let f0(x) = a[0]x^0             a[2]x^1           + a[4]x^2           + a[6]x^3
let f1(x)             a[1]x^0           + a[3]x             + a[5]x^2           + a[7]x^3

have
f(x) = f0(x^2) + x*f1(x^2)

let w = exp(i * 2pi/n)

have
f(w^0) = f0(w^0) + w^0*f1(w^0)
f(w^1) = f0(w^2) + w^1*f1(w^2)
f(w^2) = f0(w^4) + w^2*f1(w^4)
f(w^3) = f0(w^6) + w^3*f1(w^6)
f(w^4) = f0(w^0) - w^0*f1(w^0)
f(w^5) = f0(w^2) - w^1*f1(w^2)
f(w^6) = f0(w^4) - w^2*f1(w^4)
f(w^7) = f0(w^6) - w^3*f1(w^6)

感想

再帰は理解してたけど非再帰は空で書ける気がしなかったので理解した。spaghetti source のライブラリの bit reverse のところがどうなってるのかが良く分かってなくて、いまでも良く分かってない。(2018-03-09: 二進カウンタを逆から回しているだけのようです)

変更履歴

  • 2016/10/24 21:43:誤差に弱いコードだったみたいなので修正
  • 2018/03/30 05:24:参考にしても仕方がない実装を載せていたので消しました。