フーリエ変換は大学の情報工学の授業でも習うと思うのですが，数式は追えるが何をやっているのか理解しづらいものの代表な気がします(そうじゃないですか？)．FFTは競プロにおいて数列の畳み込みなどの応用問題として出題されることがあります．今回競プロの勉強と院試勉強を兼ねてFFT関連の知識をまとめてみようと思います．

フーリエ変換はややこしいけど実はみんな大好き分割統治法の一種に過ぎないんだってことを伝えることがこの記事の目的です．

フーリエ変換が活躍する場面

数列の畳み込み
画像や音声の周波数成分分解やフィルター

公式

標本点 $x=0, 1, ..., N-1$ での関数 $f(x)$ に対する離散フーリエ変換は

F(t) = \sum_{x=0}^{N-1} f(x)\omega_N^{tx} (t = 0, 1, ..., N-1)

で，フーリエ逆変換は

f(x) = \frac{1}{N} \sum_{t=0}^{N-1}F(t)\omega_N^{-tx} (x = 0, 1, ..., N-1)

です．但し，

N

は

N=2^n

で表される数とし，

\omega_N

は1の原始

N

乗根(つまり，N乗して初めて1になる数)です．復元した

f(x)

は元の関数と同じでは無いのですが，標本点において値は一致します．フーリエ変換，フーリエ逆変換は

\mathcal{F}, \mathcal{F}^{-1}

と表記されることがあります．Texでは\mathcal{F}で書けます．

表記方法は色々あり，フーリエ逆変換の $\frac{1}{N}$ の代わりに両変換に $\frac{1}{\sqrt{N}}$ を掛けて表記する場合もあります．数学的にはそちらの方が都合がいいようですが，コンピュータで扱う際は誤差の観点から $\frac{1}{N}$ とする事が多いです． $\omega_N$ の右上にマイナスが付いていたりいなかったりしますが，フーリエ変換とフーリエ逆変換で符号が逆という事が本質で，符号は逆でもいいです( $\omega_N^{-1}$ も原始 $N$ 乗根なので)．

数列の畳み込みという観点で見ると数列 $a_0, a_1, ..., a_{N-1}$ は $f(0), f(1), ..., f(N-1)$ に対応させることでうまくいきます．畳み込みは

(f * g)(x) = \sum_{n=0}^{N-1} f(n)g(x-n)

で定義され，

\mathcal{F}_N(f*g) = \mathcal{F}_N(f) \mathcal{F}_N(g)

となります．但し， $f, g$ は周期 $N$ の周期関数とします．

何故うまくいくのか

具体的な方法を説明する前に，FFTがどうして上手く動くのかの説明をします．FFTは1の原始 $N$ 乗根が以下に示すようないい性質を持っていることを利用しています．

$0<i<N$ で $\omega_N^0, \omega_N^i, ..., \omega_N^{i(N-1)}$ を全てを足すと0となる
原始 $2^n$ 乗根の2乗は原始 $2^{n-1}$ 乗根となり，上の性質を満たす．
フーリエ変換のシグマの偶数番目と奇数番目を分離可能

1つめは，

\sum_{x=0}^{N-1} \omega_N^{ix} = \frac{\omega_N^{Ni} - 1}{\omega_N^i - 1} = 0 (0<i<N)

から分かりますね．

i = 0

の時は上式の分母が0となるので成立しませんが，任意の数の0乗は1なので足すとNになります．

2つめは $\omega_{N}^2 = \omega_{\frac{N}{2}}$ となります．ただのN乗根の場合、例えば1もN乗根ですが、上の性質を満たさないのでダメです.また、Nが奇数の場合もこれは成り立たないので $N=2^n$ が必要ということです.

3つめは $N>1$ で

F(t) = \sum_{x=0}^{\frac{N}{2}-1} f(2x)\omega_{N}^{2xt} + f(2x+1)\omega_{N}^{(2x+1)t}

= \sum_{x=0}^{\frac{N}{2}-1} f(2x)\omega_{\frac{N}{2}}^{xt} + \omega_{N}^t \sum_{x=0}^{\frac{N}{2}-1} f(2x+1)\omega_{\frac{N}{2}}^{xt}

$$ = (偶数番目のフーリエ変換) + \omega_{N}^t (奇数番目のフーリエ変換) $$

と表すことができ，奇数番目と偶数番目を分解することが出来ます．分解した後の式は $\frac{N}{2}$ 要素のFFTなので $0 \leq t < \frac{N}{2}$ ですが， $\omega_{\frac{N}{2}}^{\frac{N}{2}} = 1$ より周期性があり， $t' = t + \frac{N}{2}$ を代入すると右側の項のみ負になることが分かります．

これは分割統治(が分からない人はマージソートを想像して)で $O(N log N)$ で解けるのは想像できると思います．計算量で引っかかる人は， $N$ 要素の数列を $\frac{N}{2}$ 要素の数列2つに分割しているのがポイントです．

T(N) = 2 T(\frac{N}{2}) + \Theta (N)

となり，マージソートの計算量と同じだから $T(N) = \Theta (N \log N)$ となりますね．(漸化式の評価面倒なのでちょっとサボった)

厳密には $N=2^n$ である必要は無いようで，小さな素数の積で表せるなら実行可能なようですが， $N=2^n$ とする方が実装は簡単だと思います．

FFT

FFTは原始根を複素数の $\omega_N = e^{\frac{2\pi}{N}}$ とする方法です．

逆変換は $\omega_N^N = 1$ を用いると，

f(x) = \frac{1}{N} \sum_{t=0}^{N-1}F(t)\omega_N^{(N-t)x}= \frac{1}{N} \sum_{t=0}^{N-1}F(N-t)\omega_N^{tx}

となるので，

t=1, 2, ..., N-1

を反転させるとFFTのアルゴリズムで動きます．これは下で説明するNTTでも同様です．

NTT

FFTは実数を用いて計算を行うため，誤差が生じます．NTTはmod上でフーリエ変換を行う手法です．

整数論でおなじみのフェルマーの小定理より素数 $p$ のmod上でpで割り切れない任意の整数 $g$ は

g^{p-1} \equiv 1 \bmod p

となります．また，原始根の存在定理より原始根 $r$ に対して

r^{p-1} \equiv 1 \bmod p

となります(原始根は $r^0, r^1, ..., r^{p-2}$ が全て異なるようなr)．3以上の素数は $p = m 2^n + 1$ という形で表すことができるため，

(r^m)^{2^n} \equiv 1 \bmod p

となります．原始根の定義と合わせると，1の原始 $2^n$ 乗根は $r^m$ ということになります．つまり素数pと原始根が分かっていれば，それを $\frac{p-1}{N}$ 乗したものが $N$ 乗根となるわけです．賢いですね．

例えば998244353の場合は $2^{23} \times 119 + 1$ と表すことができ，原始根は3です．要素数 $8\times 10^6$ 程度の畳み込みが可能です．1e9+7の場合は $2^1 \times 500000003+1$ となり，2乗根までしか存在しません．

C++によるソース

重要な部分だけ書きます．実際に動かす時はTypeをFFTの場合はcomplexに，NTTの場合はmodintに変える必要があります．

vector<Type> fft(vector<Type>data, int n, Type w) { // 数列, 要素数, 原始N乗根
	if (n == 1) return data;
	vector<Type> even(n/2), odd(n/2);
	for (int i = 0; i < n/2; i++) {
		even[i] = data[2*i];
		odd[i] = data[2*i+1];
	}
	vector<Type> e = ftt(even, n/2, w*w), o = ftt(odd, n/2, w*w);
	Type wmul = 1;
	for (int i = 0; i < n; i++) {
		data[i] = e[i % (n/2)] + wmul * o[i % (n/2)];
		wmul *= w;
	}
	return data;
}

応用1 任意mod畳み込み

複数のmodに対してNTTを行い，中国剰余定理で復元させれば良いです．Garnerのアルゴリズムという名前が付いているらしいです．1回の畳み込みで2個くらいmod取れば良さそうです．

あるいはFFTで整数を15bitずつに分解すると3回のFFTで誤差無く実行できるようです．

応用2 形式冪級数

多項式の整数性を利用して，組み合わせを多項式で表現する手法です．例えば1回で2, 3, 5歩進む時，n回でm歩進んでいる場合の数を $(x^2+x^3+x^5)^n$ の $x^m$ の係数とみなせます．

多項式の掛け算をFFTなどを用いた畳み込みで行うことも出来ます．詳しくはmaspyさんの解説が詳しいです．

応用3 GCD畳み込み

フーリエ変換全く関係ないです．

c_k = \sum_{gcd(i, j) = k} a_i b_i

この畳み込みの計算を高速に解く方法を考えます． $f(x)_k$ を数列xの中で添字がkの倍数の和を取る関数を定義すると，

f(c)_i = f(a)_i \ast f(b)_i

となりますね． $\ast$ は要素ごとの積です． $f^{-1}$ は上から引いていけばよいので求めることができ，

c = f^{-1}(f(a) * f(b))

で求めることが出来ます．ん〜釈然としない...より詳細の解説はnoshiさんのブログから．一般化すると，緩い条件で簡単に数え上げることが出来るなら，それで計算して後から逆変換を掛けられないか考えるってところかな...

応用4 mod畳み込み

c_k = \sum_{ij \equiv k (\bmod p)} a_i b_j

を求めよ．出典AGC047 C

まずは添字がmodの個数の数列に変換し，これの畳み込みの形を考える．普通の畳み込みは添字部分は足し算だけど，今回はmod上の掛け算となっている．

掛け算は扱いづらいのでどうにか足し算に出来ないか考えてみる．原始根のi乗と考えると，○乗部分は足し算になるので畳み込める．頭いいですね．本番では思いつけなかった．

バタフライ演算

みんな大ッキライバタフライ演算です．とは言え分割統治が頭にあると，再帰により下の桁から順番に分岐することが分かるので，逆算する時は上の桁から分岐させればいいことが分かります．なので2進数でbitを反対から取ったbit-reverseの順に最初に並べれば良さそうな感じがします．bitが立っている方(下図で赤色)の結果には $\omega^?$ を掛けます．

図はこのサイトから引用しました．言われたら分かるけど，絶対自分では思い付けねぇ...

FFTを人間が解くのは無意味だし，コンピュータが解くなら分割統治で良いと思うんだけど，バタフライ演算の方が定数倍が高速とかあるのかな？分からないです．

終わりに

疲れました．間違いがあったら教えて下さい．

modの問題で原始根に帰着させると見通しが良くなることがあるので知っておくといいと思いました。