FFT(Fast Fourier Transform),高速フーリエ変換についての記事です。
非常に理解が難しいアルゴリズムとして有名なので、どうにかして理解したいという方に向けての記事です。
説明はなるべく丁寧に行っていますので、わからないところがあったらコメント欄へどうぞ。
数Ⅱまでの知識を前提にして、それ以降は個別に説明しているつもりですが、抜けているところがあれば言ってください。
この記事を書くにあたり、AtCoderの高速フーリエ変換を参考にさせていただきました。AtCoder社の方々、ありがとうございます。
高速フーリエ変換とは
離散フーリエ変換という処理を高速に行うアルゴリズムです。これを利用して、多項式乗算が高速に行えます。ここでは、多項式乗算を行う方法論をメインに、高速フーリエ変換の解説をします。また、ふたつの多項式の次数は便宜上同じとします。（違う場合、高次部分の係数を $0$ にするなどして合わせればよいです。）
求めるべきもの
f(x)=\sum_{i=0}^{N} A_ix^i,g(x)=\sum_{i=0}^{N} B_ix^i
を乗算することを考えます。
この積は、
f(x)g(x)=\sum_{i=0}^{N} \sum_{j=0}^{N} A_iB_jx^{i+j}
であり、$A_n$ から $A_{2n-1}$, $B_n$ から $B_{2n-1}$ を便宜上 $0$ として、
C_k=\sum_{i=0}^{k} A_iB_{k-i}
とすると、
f(x)g(x)=\sum_{i=0}^{2N} C_ix^i
となります。
この $C_i$ を、$0\leq i\leq 2N$ を満たすすべての $i$ について求めるのが今回の目標です。
それぞれ愚直に求めると、$f,g$ の全項を組み合わせて参照することになるので、 $O(N^2)$ です。これをどうにかして高速化します。
多項式補間
愚直な乗算は難しそうなので、$C_i$ の値を、多項式補間を用いて算出することを考えます。
多項式補間とは、多項式の変数に実際にいくつかの値を代入し、多項式を計算した値から、多項式の係数を決定する手法です。
たとえば、$f(x)=ax+b$ という $1$ 次関数があるとします。
$a$ と $b$ の値は分かりませんが、$f(3)=5,f(7)=-3$ がわかっているものとします。
実際に $3,7$ を代入してみると、
- $3a+b=5$
- $7a+b=-3$
と、連立方程式が立ち、$a,b$ の値が求められます。
一般的には、$n$ 次多項式について、$n+1$ 個の異なる変数値においての値が定まっていれば、多項式を決定できることが知られています。
ガウスの消去法などを使って愚直に連立方程式を解くと $O(n^3)$ かかりますが、これはなんとかしてあとで改善します。
以上より、高速な多項式乗算の全体的な流れは次のようになります。
- $f,g$ の次数の和より大きい適当な値 $n$ をとり、$n$ 個の点 $x_0,x_1\cdots x_{n-1}$ をうまく選ぶ。
- $f(x_0),f(x_1)\cdots f(x_{n-1}),g(x_0),g(x_1)\cdots g(x_{n-1})$ を計算する。
- $f(x_0)g(x_0),f(x_1)g(x_1)\cdots f(x_{n-1})g(x_{n-1})$ を計算する。
- 3.で求めた値をもとに、なんとかして $C_i (0\leq i\leq 2N)$ を復元する。
高速フーリエ変換を用いて、2. と4. を高速に行うのが実際のアルゴリズムです。
評価する点の選び方
さて、多項式補間に使う $n$ と点 $x_i (0\leq i<n)$ はどう選べばいいのでしょうか。
結論から言うと、$n$ はなんらかの2べきの数とし、$x_i$ には、 $1$ の $n$ 乗根を使います。
多くの解説だと（たぶん）数Ⅲの知識が前提にされているのでこのあたりは軽く飛ばされますが、今回はその説明も詳細にやります。
複素平面
$\sqrt{-1}$ を虚数単位 $i$ として定め、$a+bi$ として表される数が複素数と呼ばれるのは、ご存じだと思います。
$x,y$ が実数のとき、$x+yi$ を単なる実数の組としてみなし、平面の直交座標 $(x,y)$ の点に対応付けるのです。
この複素数を、平面上の点、または位置ベクトルに対応させるとき、この平面を複素平面と呼びます。
複素平面上では、複素数に対する代数的演算が、平面上での幾何的な操作と対応付けられます。
たとえば、加法は平行移動、実数倍は原点中心の拡大・縮小に対応するなどします。
複素平面と積
加法や実数倍と同じように、複素数同士の積も、複素平面上で幾何的な操作として扱えます。
複素数 $z=x+yi$ を、$r(\cos\theta +i\sin\theta)$ として表すことを考えます。この表示を $z$ の極形式と呼びます。
このとき、$r=\sqrt{x^2+y^2}$ が成り立ちます。つまり、$r$ は $z$ の絶対値となります。
また、$\theta$ を $z$ の偏角と呼び、 $\arg z$ と表されます。
$r,\theta$ から $z$ は一意に定まります。具体的には、$z$ は複素平面上で、$(r,0)$ を $\theta$ 回転して得られます。
ふたつの複素数 $z$,$w$ の極形式をそれぞれ $r(\cos\alpha +i\sin\alpha)$, $s(\cos\beta +i\sin\beta)$ とすると、三角関数の加法定理より、
$(\cos\alpha +i\sin\alpha)(\cos\beta +i\sin\beta)=\cos(\alpha +\beta)+i\sin(\alpha +\beta)$
が成り立つので、
$zw=rs(\cos(\alpha +\beta)+i\sin(\alpha +\beta))$ となり、
$zw$ は、点 $z$ を、原点を中心に $\beta$ 回転、$s$ 倍に拡大して得られる点だとわかります。
このように、複素数の積は、複素平面上で、点の回転と拡大を組み合わせて表すことができます。
また、互いに共役な複素数は実軸に対して対称な位置にある、などの性質もあります。
1 の n 乗根
ここまでの話を踏まえると、$1$ の $n$ 乗根は次のように考えられます。
$n\theta=2m\pi (m\in\mathbb{Z},0\leq\theta<2\pi)$ をみたすような $\theta$ は $n$ 個存在し、
$z=\cos\theta+i\sin\theta$ を満たす $z$ は、$1$ の $n$ 乗根である。
ここで、 $\theta=\frac{2\pi}{n}$ であるときの $z$ を $\zeta_n$ とすると、$1$ の $n$ 乗根全体は $\zeta_n^i (0\leq i<n)$ で表されます。
1 の n 乗根の性質
$\zeta_n$ は次のような性質を持ちます。
\zeta_n^i=\zeta_n^j\Leftrightarrow i\equiv j\pmod n\\
\sum_{i=0}^{n-1}\zeta_n^{i(j-k)}=\left\{
\begin{array}{ll}
n & {\rm if}\ j\equiv k\pmod n\\
0 & \rm otherwise
\end{array}
\right.
上は自明でしょう。
下の式は証明します。
$j=k$ のとき、偏角の合計は $2m\pi(m\in\mathbb{N})$ で、絶対値の和は $n$ なので、これが $n$ と等しいことは明らかです。
そうでない場合、この値は、初項 $1$, 公比 $\zeta_n^{j-k}$ の等比数列の最初から $n$ 項の和になります。
これを等比数列の和の公式、「初項 $a,$ 公比 $r$ の等比数列最初の $n$ 項の和は $a(1-r^n)/(1-r)$」に代入すると、$r$ が $1$ の $n$ 乗根であることより、$0$ になるので、証明できました。
また、$\zeta_n$ を $\zeta_n^{-1}$ で置き換えてもこれらの性質は変わりません。これは後に重要になってきます。
離散フーリエ変換
多項式 $f(x)$ に対し、その離散フーリエ変換(DFT,Discrete Fourier Transform), $\widehat{f}(x)$を次のように定義します。
\widehat{f}(t)=\sum_{i=0}^{n-1}f(\zeta_n^i)t^i
つまり、$\widehat{f}(t)$ は、多項式 $f(x)$ の $x$ に $1$ の $n$ 乗根を代入し、得られた値を係数にもつ多項式です。
Wikipediaなどを見に行くと、離散フーリエ変換は複素関数から複素関数への写像である、というふうな記述があると思いますが、これでずいぶん分かりやすくなったと思います。
離散フーリエ変換の性質
f(x)=\sum_{j=0}^{n-1}c_jx^j
とすると、その離散フーリエ変換、$\widehat{f}$ は、
\begin{eqnarray}
\widehat{f}(t)
&=&\sum_{i=0}^{n-1}f(\zeta_n^i)t^i\\
&=&\sum_{i=0}^{n-1}(\sum_{j=0}^{n-1}c_j(\zeta_n^i)^j)t^i\\
&=&\sum_{j=0}^{n-1}c_j\sum_{i=0}^{n-1}(\zeta_n^jt)^i
\end{eqnarray}
になります。
ここで、$t$ に $\zeta_n^{-k}$ を代入してみます。
\begin{eqnarray}
\widehat{f}(\zeta_n^{-k})
&=&\sum_{j=0}^{n-1}c_j\sum_{i=0}^{n-1}(\zeta_n^j\zeta_n^{-k})^i\\
&=&\sum_{j=0}^{n-1}c_j\sum_{i=0}^{n-1}\zeta_n^{i(j-k)}
\end{eqnarray}
となります。
ここで、前に挙げた $\zeta_n$ の性質を使います。
\sum_{i=0}^{n-1}\zeta_n^{i(j-k)}=\left\{
\begin{array}{ll}
n & {\rm if}\ j\equiv k\pmod n\\
0 & \rm otherwise
\end{array}
\right.
と合わせて考えると、結局、
\widehat{f}(\zeta_n^{-k})=nc_k
です。
よって、離散フーリエ変換をした関数に、離散フーリエ変換と同じ処理を、 $\zeta_n$ の部分を $\zeta_n^{-1}$ に置き換えて適用し、係数を $\frac{1}{n}$ 倍するともとの多項式に戻ります。これを、離散フーリエ逆変換と呼びます。
また、離散フーリエ変換後の係数は、$1$ の $n$ 乗根を多項式に代入した値なので、$f$ と $g$ の積の離散フーリエ変換は、単純に $f,g$ それぞれの離散フーリエ変換の、同次の係数を掛け合わせたものになります。
これで、離散フーリエ変換の性質を利用して方針が固められました。
- $f,g$ の次数の和より大きい適当な値 $n$ をとり、$f,g$ それぞれについて、$1$ の $n$ 乗根を使って離散フーリエ変換を計算する。
- 変換後の関数について同次項の係数を掛け合わせ、新しい多項式を構成する。これは、$f$ と $g$ の積の離散フーリエ変換と一致する。
- これに対して離散フーリエ逆変換を行い、$f$ と $g$ の積が復元できる。
こうして、離散フーリエ変換を高速に計算することで多項式乗算が行えることが示せました。FFTの本質はここからです。
高速フーリエ変換
ここまでで多項式乗算の準備が整いました。あとは離散フーリエ変換を高速に計算するだけです。
$n$ を $2$ の冪乗として、$n-1$ 次多項式 $f(x)=\sum_{i=0}^{n-1}c_ix^i$ に対し、$f_0,f_1$ を次のように定義します。
f_0(x)=\sum_{i=0}^{n/2-1}c_{2i}x^i\\
f_1(x)=\sum_{i=0}^{n/2-1}c_{2i+1}x^i
このとき、$f(x)=f_0(x^2)+xf_1(x^2)$ が成り立ちます。
よって、$f_0(\zeta_n^0),f_0(\zeta_n^2),\cdots f_0(\zeta_n^{2(n-1)})$ と$f_1(\zeta_n^0),f_1(\zeta_n^2),\cdots f_1(\zeta_n^{2(n-1)})$ が求められれば、$f(\zeta_n^0),f(\zeta_n^1),\cdots f(\zeta_n^{n-1})$ も求めることができます。
$\zeta_n^2=\zeta_{n/2}$ なので、（これは明らかです）
$f_0(\zeta_{n/2}^0),f_0(\zeta_{n/2}^1),\cdots f_0(\zeta_{n/2}^{n-1})$ と $f_1(\zeta_{n/2}^0),f_1(\zeta_{n/2}^1),\cdots f_1(\zeta_{n/2}^{n-1})$ が求められればいいことになります。
そして、$\zeta_{n/2}$ は $n/2$ 乗すると $1$ になることも考えると、上に挙げたものはそれぞれ前半と後半が同じです。よって、求めなければいけないのは結局、
$f_0(\zeta_{n/2}^0),f_0(\zeta_{n/2}^1),\cdots f_0(\zeta_{n/2}^{n/2-1})$ と $f_1(\zeta_{n/2}^0),f_1(\zeta_{n/2}^1),\cdots f_1(\zeta_{n/2}^{n/2-1})$ になります。
これは、元の問題のちょうど半分のサイズの全く同じ問題です。
よって、元の問題の半分のサイズの同じ問題を2つ解くと、$O(n)$ で元の問題の答えも求まります。
こうして問題を再帰的に解いていけば、$O(n\log n)$ で離散フーリエ変換が行えます！！！！！（マージソートなどと同じような仕組みの分割統治法なので、わからなければ調べてみましょう。）
これが高速フーリエ変換というアルゴリズムです。
これでここまでの長い解説もやっと終わり。$O(n\log n)$ での多項式乗算が完成しました！
ここまでをまとめます。
- 次数の和が $n$ 以下の多項式 $f,g$ の積を求めたい。
- $f,g$ 両方について $1$ の $n$ 乗根の値を代入し、出てきた値を求める。これは離散フーリエ変換で求められる。
- 上で求めた値を掛け合わせると、$f,g$ の積の離散フーリエ変換が求められる。
- 離散フーリエ逆変換を $f,g$ の積の離散フーリエ変換に行い、求めたい多項式を復元できる。
こうして、多項式乗算を $O(N\log N)$ で行うことができました。お疲れさまでした。
サンプルコードは以下です。（ここでは、高速化のために自作した複素数クラスを用いています。その点では互換性がないので、基本的にはstd::complexを使っていただければ結構です。とはいえ、必要最低限の機能を持った複素数クラスを自作すると実行時間が半分程度まで削れるので、それもお勧めします。）
class FastFourierTransform {
private:
static void dft(vector<mycomplex>& func, int inverse) {
int sz = func.size();
if (sz == 1)return;
vector<mycomplex> veca, vecb;
rep(i, sz / 2) {
veca.push_back(func[2 * i]);
vecb.push_back(func[2 * i + 1]);
}
dft(veca, inverse); dft(vecb, inverse);
mycomplex now = 1, zeta = polar(1.0, inverse * 2.0 * acos(-1) / sz);
rep(i, sz) {
func[i] = veca[i % (sz / 2)] + now * vecb[i % (sz / 2)];
now *= zeta;
}
}
public:
template<typename T>
static vector<double> multiply(vector<T> f, vector<T> g) {
vector<mycomplex> nf, ng;
int sz = 1;
while (sz < f.size() + g.size())sz *= 2;
nf.resize(sz); ng.resize(sz);
rep(i, f.size()) {
nf[i] = f[i];
ng[i] = g[i];
}
dft(nf, 1);
dft(ng, 1);
rep(i, sz)nf[i] *= ng[i];
dft(nf, -1);
vector<double> res;
rep(i, sz)res.push_back(nf[i].real() / sz);
return res;
}
};
おわりに
長い記事ですが、読んでいただけたでしょうか。
自分もこれを理解するまでには長い時間を要しましたが、理解できた時にはかなりすっきりとしました。
みなさんにも、これが理解の助けとなる時が来てほしい、という思いで記事を書きました。
今はわからなくても、これが少しでも理解に貢献できたなら幸いです。
そして、この記事の題名にもあるように、「完全に理解する」という意気で、さまざまなアルゴリズムに貪欲に挑戦していってほしいと思います。
ありがとうございました。