こんにちは、20B の @tatyam です。この記事は、traP夏のブログリレー 33日目の記事です。

高速フーリエ変換 といえば、競技プログラミングで数列の畳み込みなどに使われるテクニックですね。
今日は高速フーリエ変換と高速フーリエ変換を利用した畳み込みを習得していきましょう。

ポエム

一般的なフーリエ変換の 定義 / 説明 と異なる部分があるかもしれませんが、競プロであればこれで十分だと思います。
AtCoder では AC Library さえ使えれば FFT は理解しなくて良いことになっていますね。
まあでも、中身を理解するといろいろと応用ができるようになるので、理解する価値はかなり高いでしょう。

離散フーリエ変換 (DFT)

離散フーリエ変換 とは、高速フーリエ変換の元になっているものです。離散フーリエ変換の特殊な場合を高速化したものが高速フーリエ変換というわけですね。

$$\begin{array}{rl}& \text{DFT}\\ & \downarrow \text{高速化}\\ & \text{FFT}\end{array}$$

まずは離散フーリエ変換から理解しましょう。

定義

$2\pi$ と書くのはちょっと長いので、$\tau =2\pi$ とします。
${\zeta }_N$ を $1$ の原始 $N$ 乗根のうちの $1$ つ、${\zeta }_N=e^{i\tau /N}=\cos (\frac{\tau }{N})+i\sin (\frac{\tau }{N})$ とします。
例えば、${\zeta }_2=-1,{\zeta }_3=\frac{\sqrt{3}-i}{2},{\zeta }_4=i,{\zeta }_{2N}^2={\zeta }_N,{\zeta }_N^N=1$ です。

複素平面上の単位円を $N$ 等分する点

数列 $A=(A_0,A_1,\dots ,A_{N-1})$ から $N-1$ 次多項式 $f_A(x)=A_0+A_{1}x+\cdots +A_{N-1}{x}^{N-1}$ を作ります。
畳み込みと多項式の積は対応しているので、この 2 つは同一視しがちです。

$$\begin{array}{rl}{B}_0& =f_A({\zeta }_N^0)={\zeta }_N^0{A}_0+{\zeta }_N^0{A}_1+{\zeta }_N^0{A}_2+\cdots +{\zeta }_N^0{A}_{N-1}\\ B_1& =f_A({\zeta }_N^1)={\zeta }_N^0{A}_0+{\zeta }_N^1{A}_1+{\zeta }_N^2{A}_2+\cdots +{\zeta }_N^{N-1}{A}_{N-1}\\ B_2& =f_A({\zeta }_N^2)={\zeta }_N^0{A}_0+{\zeta }_N^2{A}_1+{\zeta }_N^4{A}_2+\cdots +{\zeta }_N^{2(N-1)}{A}_{N-1}\\ & ⋮\end{array}$$

$$\left[\begin{array}{c}{B}_0\\ B_1\\ B_2\\ ⋮\\ B_{N-1}\end{array}\right]=\left[\begin{array}{ccccc}{\zeta }_N^0& {\zeta }_N^0& {\zeta }_N^0& \cdots & {\zeta }_N^0\\ {\zeta }_N^0& {\zeta }_N^1& {\zeta }_N^2& \cdots & {\zeta }_N^{N-1}\\ {\zeta }_N^0& {\zeta }_N^2& {\zeta }_N^4& \cdots & {\zeta }_N^{2(N-1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\zeta }_N^0& {\zeta }_N^{N-1}& {\zeta }_N^{2(N-1)}& \cdots & {\zeta }_N^{(N-1{)}^2}\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_1\\ A_2\\ ⋮\\ A_{N-1}\end{array}\right]$$

離散フーリエ変換の行列表現

vector<complex<double>> DFT(vector<complex<double>> A){
    const int N = A.size();
    vector<complex<double>> B(N);
    for(int i = 0; i < N; i++) for(int j = 0; j < N; j++){
        B[i] += A[j] * polar(1, 2 * M_PI * i * j / N);
    }
    return B;
}

C++ で離散フーリエ変換

逆離散フーリエ変換

逆に、$B=(f_A({\zeta }_N^0),f_A({\zeta }_N^1),\dots ,f_A({\zeta }_N^{N-1}))$ から $A=(A_0,A_1,\dots ,A_{N-1})$ を復元する操作を逆離散フーリエ変換と言います。
逆離散フーリエ変換はどのような形になっているでしょうか？

ここで、

$$\begin{array}{rl}{B}_0+B_1+\cdots +B_{N-1}=& ({\zeta }_N^0+{\zeta }_N^0+{\zeta }_N^0+\cdots +{\zeta }_N^0)A_0\\ & +({\zeta }_N^0+{\zeta }_N^1+{\zeta }_N^2+\cdots +{\zeta }_N^{N-1})A_1\\ & +({\zeta }_N^0+{\zeta }_N^2+{\zeta }_N^4+\cdots +{\zeta }_N^{2(N-1)})A_2\\ & +\cdots \\ & +({\zeta }_N^0+{\zeta }_N^{N-1}+{\zeta }_N^{2(N-1)}+\cdots +{\zeta }_N^{(N-1{)}^2})A_{N-1}\\ =& N{A}_0\end{array}$$

という式が成り立っています。（ $A_0$ 以外の係数は複素平面上の単位円をグルグル回っているので、全部足すと $0$ ）

よって、$A_0=\frac{1}{N}(B_0+B_1+\cdots +B_{N-1})$ と復元できます。
また、

$$\begin{array}{rl}{\zeta }_N^0{B}_0+{\zeta }_N^{-1}{B}_1+\cdots +{\zeta }_N^{-N+1}{B}_{N-1}=& ({\zeta }_N^0+{\zeta }_N^{-1}+{\zeta }_N^{-2}+\cdots +{\zeta }_N^{-N+1})A_0\\ & +({\zeta }_N^0+{\zeta }_N^0+{\zeta }_N^0+\cdots +{\zeta }_N^0)A_1\\ & +({\zeta }_N^0+{\zeta }_N^1+{\zeta }_N^2+\cdots +{\zeta }_N^{N-1})A_2\\ & +\cdots \\ & +({\zeta }_N^0+{\zeta }_N^{N-2}+{\zeta }_N^{2(N-2)}+\cdots +{\zeta }_N^{(N-1)(N-2)})A_{N-1}\\ =& N{A}_1\end{array}$$

なので、$A_1=\frac{1}{N}({\zeta }_N^0{B}_0+{\zeta }_N^{-1}{B}_1+\cdots +{\zeta }_N^{-N+1}{B}_{N-1})$ と復元できます。

同様の計算により、逆離散フーリエ変換は、離散フーリエ変換の ${\zeta }_N$ を ${\zeta }_N^{-1}$ に置き換えて、最後に $N$ で割ったものであることが分かります。

$$\left[\begin{array}{c}{A}_0\\ A_1\\ A_2\\ ⋮\\ A_{N-1}\end{array}\right]=\frac{1}{N}\left[\begin{array}{ccccc}{\zeta }_N^0& {\zeta }_N^0& {\zeta }_N^0& \cdots & {\zeta }_N^0\\ {\zeta }_N^0& {\zeta }_N^{-1}& {\zeta }_N^{-2}& \cdots & {\zeta }_N^{-N+1}\\ {\zeta }_N^0& {\zeta }_N^{-2}& {\zeta }_N^{-4}& \cdots & {\zeta }_N^{2(-N+1)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\zeta }_N^0& {\zeta }_N^{-N+1}& {\zeta }_N^{2(-N+1)}& \cdots & {\zeta }_N^{-(N-1{)}^2}\end{array}\right]\left[\begin{array}{c}{B}_0\\ B_1\\ B_2\\ ⋮\\ B_{N-1}\end{array}\right]$$

逆離散フーリエ変換の行列表現

vector<complex<double>> IDFT(vector<complex<double>> B){
    const int N = B.size();
    vector<complex<double>> A(N);
    for(int i = 0; i < N; i++) for(int j = 0; j < N; j++){
        A[i] += B[j] * polar(1, -2 * M_PI * i * j / N);
    }
    for(int i = 0; i < N; i++) A[i] /= N;
    return A;
}

C++ で逆離散フーリエ変換

離散フーリエ変換を利用した畳み込み

数列 $A=(A_0,A_1,\dots ,A_{|A|-1}),B=(B_0,B_1,\dots ,B_{|B|-1})$ から、

$$C_k=\sum _{i+j=k}{A}_i{B}_j$$

で定義される数列 $C=(C_0,C_1,\dots ,C_{|A|+|B|-2})$ を作ることを畳み込みと言います。例えば、$A=(3,1,4)$ と $B=(1,5,9)$ を畳み込むと、$C=(3,15+1,27+5+4,9+20,36)=(3,16,36,29,36)$ となります。

畳み込みは定義式をそのまま計算することで $\Theta (|A||B|)$ で計算することができますが、離散フーリエ変換を経由しても計算することができます。

畳み込みと多項式の積

畳み込みの操作は、多項式の積と対応しています。

列 $X$ に対して $f_X(x)$ を $|X|-1$ 次多項式 $f_X(x)=X_0+X_{1}x+\cdots +X_{|X|-1}{x}^{|X|-1}$ とします。

$A$ と $B$ を畳み込んで $C$ ができるとき、$f_A(x)\times f_B(x)=f_C(x)$ が成り立ちます。

$i+j=k$ の条件を $x$ の指数にすることで、$x^i$ と $x^j$ を掛けたときに $x^{i+j}$ と足し算ができて、同類項が足されて畳み込みになるというわけですね。

多点評価と多項式補間

グラフ上で、$1$ 点を通る $0$ 次多項式、$2$ 点を通る $1$ 次多項式、$3$ 点を通る $2$ 次多項式、… は一意に定まりますね。
$f_C(x)$ は $|A|+|B|-2$ 次多項式なので、$|A|+|B|-1$ 個の点での$f_C(x)$ の値が分かれば $f_C(x)$ が分かるということになります。

$N=|A|+|B|-1$ と置きます。
$N$ 個の点 ${\zeta }_N^0,{\zeta }_N^1,\dots ,{\zeta }_N^{N-1}$ を選び、$f_C({\zeta }_N^0),f_C({\zeta }_N^1),\dots ,f_C({\zeta }_N^{N-1})$ を求めて $C$ を復元することを考えます。これが畳み込みを高速にする秘訣です。

$f_A(x)\times f_B(x)=f_C(x)$ であることを思い出すと、$f_A({\zeta }_N^0),f_A({\zeta }_N^1),\dots ,f_A({\zeta }_N^{N-1})$ と $f_B({\zeta }_N^0),f_B({\zeta }_N^1),\dots ,f_B({\zeta }_N^{N-1})$ が分かれば、$f_C({\zeta }_N^0),f_C({\zeta }_N^1),\dots ,f_C({\zeta }_N^{N-1})$ は求まります。

$f_A({\zeta }_N^0),f_A({\zeta }_N^1),\dots ,f_A({\zeta }_N^{N-1})$ を求めるには、離散フーリエ変換です。長さ $N$ になるまで $A$ の末尾に $0$ を追加してから離散フーリエ変換をすれば $f_A({\zeta }_N^0),f_A({\zeta }_N^1),\dots ,f_A({\zeta }_N^{N-1})$ が出てきます。
同様にして $f_B({\zeta }_N^0),f_B({\zeta }_N^1),\dots ,f_B({\zeta }_N^{N-1})$ を求めて、掛け算をして $f_C({\zeta }_N^0),f_C({\zeta }_N^1),\dots ,f_C({\zeta }_N^{N-1})$ を求めます。
最後に、$f_C({\zeta }_N^0),f_C({\zeta }_N^1),\dots ,f_C({\zeta }_N^{N-1})$ を逆離散フーリエ変換して、$C$ が求まりました。

計算量は、離散フーリエ変換が支配的なので、ここが高速になれば高速な畳み込みが実現できそうです。

高速フーリエ変換 (FFT)

$N$ が $2$ 冪なら高速になるらしいですよ…

実験してみよう

試しに $N=4$ のときで考えてみましょう。

$$\begin{array}{rl}{B}_0& =A_0+A_1+A_2+A_3\\ B_1& =A_0+i{A}_1-A_2-i{A}_3\\ B_2& =A_0-A_1+A_2-A_3\\ B_3& =A_0-i{A}_1-A_2+i{A}_3\end{array}$$

どこかまとめて計算できる場所はあるでしょうか？

$B_0$ と $B_2$ 、$B_1$ と $B_3$ が似ていますね…

$$\begin{array}{rl}{B}_0& =A_0+A_1+A_2+A_3\\ B_2& =A_0-A_1+A_2-A_3\\ B_1& =A_0+i{A}_1-A_2-i{A}_3\\ B_3& =A_0-i{A}_1-A_2+i{A}_3\end{array}$$

$$\begin{array}{rl}{B}_0& =(A_0+A_2)+(A_1+A_3)\\ B_2& =(A_0+A_2)-(A_1+A_3)\\ B_1& =(A_0-A_2)+i(A_1-A_3)\\ B_3& =(A_0-A_2)-i(A_1-A_3)\end{array}$$

$A_0+A_2,A_0-A_2,A_1+A_3,A_1-A_3$ の $4$ つにまとまりました。

$N=8$

$N=8$ でもやってみましょう。表記が面倒なので行列を使います。

$$\left[\begin{array}{c}{B}_0\\ B_4\\ B_1\\ B_5\\ B_2\\ B_6\\ B_3\\ B_7\end{array}\right]=\left[\begin{array}{cccccccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^1& {\zeta }_8^2& {\zeta }_8^3& {\zeta }_8^4& {\zeta }_8^5& {\zeta }_8^6& {\zeta }_8^7\\ {\zeta }_8^0& {\zeta }_8^5& {\zeta }_8^2& {\zeta }_8^7& {\zeta }_8^4& {\zeta }_8^1& {\zeta }_8^6& {\zeta }_8^3\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6& {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2& {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\\ {\zeta }_8^0& {\zeta }_8^3& {\zeta }_8^6& {\zeta }_8^1& {\zeta }_8^4& {\zeta }_8^7& {\zeta }_8^2& {\zeta }_8^5\\ {\zeta }_8^0& {\zeta }_8^7& {\zeta }_8^6& {\zeta }_8^5& {\zeta }_8^4& {\zeta }_8^3& {\zeta }_8^2& {\zeta }_8^1\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\\ A_7\end{array}\right]$$

似ている部分を揃えて…

$$\left[\begin{array}{c}{B}_0\\ B_4\\ B_1\\ B_5\\ B_2\\ B_6\\ B_3\\ B_7\end{array}\right]=\left[\begin{array}{cccccccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& -{\zeta }_8^0& {\zeta }_8^0& -{\zeta }_8^0& {\zeta }_8^0& -{\zeta }_8^0& {\zeta }_8^0& -{\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^1& {\zeta }_8^2& {\zeta }_8^3& {\zeta }_8^4& {\zeta }_8^5& {\zeta }_8^6& {\zeta }_8^7\\ {\zeta }_8^0& -{\zeta }_8^1& {\zeta }_8^2& -{\zeta }_8^3& {\zeta }_8^4& -{\zeta }_8^5& {\zeta }_8^6& -{\zeta }_8^7\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6& {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& -{\zeta }_8^2& {\zeta }_8^4& -{\zeta }_8^6& {\zeta }_8^0& -{\zeta }_8^2& {\zeta }_8^4& -{\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^3& {\zeta }_8^6& {\zeta }_8^1& {\zeta }_8^4& {\zeta }_8^7& {\zeta }_8^2& {\zeta }_8^5\\ {\zeta }_8^0& -{\zeta }_8^3& {\zeta }_8^6& -{\zeta }_8^1& {\zeta }_8^4& -{\zeta }_8^7& {\zeta }_8^2& -{\zeta }_8^5\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_1\\ A_2\\ A_3\\ A_4\\ A_5\\ A_6\\ A_7\end{array}\right]$$

マイナスを取り出して…

$$\left[\begin{array}{c}{B}_0\\ B_4\\ B_1\\ B_5\\ B_2\\ B_6\\ B_3\\ B_7\end{array}\right]=\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_2\\ A_4\\ A_6\end{array}\right]+\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ -{\zeta }_8^0& -{\zeta }_8^0& -{\zeta }_8^0& -{\zeta }_8^0\\ {\zeta }_8^1& {\zeta }_8^3& {\zeta }_8^5& {\zeta }_8^7\\ -{\zeta }_8^1& -{\zeta }_8^3& -{\zeta }_8^5& -{\zeta }_8^7\\ {\zeta }_8^2& {\zeta }_8^6& {\zeta }_8^2& {\zeta }_8^6\\ -{\zeta }_8^2& -{\zeta }_8^6& -{\zeta }_8^2& -{\zeta }_8^6\\ {\zeta }_8^3& {\zeta }_8^1& {\zeta }_8^7& {\zeta }_8^5\\ -{\zeta }_8^3& -{\zeta }_8^1& -{\zeta }_8^7& -{\zeta }_8^5\end{array}\right]\left[\begin{array}{c}{A}_1\\ A_3\\ A_5\\ A_7\end{array}\right]$$

同じ部分が出来ています！

$$\left[\begin{array}{c}{B}_0\\ B_4\\ B_1\\ B_5\\ B_2\\ B_6\\ B_3\\ B_7\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_2\\ A_4\\ A_6\end{array}\right]+\left[\begin{array}{cccc}{\zeta }_8^0& 0& 0& 0\\ -{\zeta }_8^0& 0& 0& 0\\ 0& {\zeta }_8^1& 0& 0\\ 0& -{\zeta }_8^1& 0& 0\\ 0& 0& {\zeta }_8^2& 0\\ 0& 0& -{\zeta }_8^2& 0\\ 0& 0& 0& {\zeta }_8^3\\ 0& 0& 0& -{\zeta }_8^3\end{array}\right]\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\end{array}\right]\left[\begin{array}{c}{A}_1\\ A_3\\ A_5\\ A_7\end{array}\right]$$

同じ部分をまとめた結果、

$\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\end{array}\right]\left[\begin{array}{c}{A}_0\\ A_2\\ A_4\\ A_6\end{array}\right]$ と $\left[\begin{array}{cccc}{\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0& {\zeta }_8^0\\ {\zeta }_8^0& {\zeta }_8^2& {\zeta }_8^4& {\zeta }_8^6\\ {\zeta }_8^0& {\zeta }_8^4& {\zeta }_8^0& {\zeta }_8^4\\ {\zeta }_8^0& {\zeta }_8^6& {\zeta }_8^4& {\zeta }_8^2\end{array}\right]\left[\begin{array}{c}{A}_1\\ A_3\\ A_5\\ A_7\end{array}\right]$

を計算すれば良いということが分かりました。
これはなんと、離散フーリエ変換の形をしています。(すごい)
長さが半分になっていくので、$N$ が $2$ 冪であれば、これを繰り返し適用できて、高速になりそうです。

一般化

$N$ を $2$ 冪とします。$A=(A_0,A_1,\dots ,A_{N-1})$ を高速フーリエ変換して $B=(B_0,B_1,\dots ,B_{N-1})$ を作ります。

1. $A$ を奇数番目と偶数番目に分けて、$A_e=(A_0,A_2,\dots ,A_{N-2}),A_o=(A_1,A_3,\dots ,A_{N-1})$ を作ります。
2. $A_e,A_o$ をそれぞれ高速フーリエ変換して、$B_e,B_o$ とします。 (長さ $1$ の場合そのまま)
3. $(B_o{)}_i$ (0-indexed) を ${\zeta }_N^i$ 倍します。
4. $B_i=(B_e{)}_i+(B_o{)}_i,B_{N/2+i}=(B_e{)}_i-(B_o{)}_i$ の式で $B$ を作ります。

長さ $N$ の処理を $1$ 回、長さ $N/2$ の処理を $2$ 回、… 、長さ $1$ の処理を $N$ 回行うので、計算量は $O(N\log N)$ です。

vector<complex<double>> FFT(vector<complex<double>> A){
    const int N = A.size();
    if(N == 1) return A;
    vector<complex<double>> even(N / 2), odd(N / 2);
    for(int i = 0; i < N / 2; i++){
        even[i] = A[i * 2];
        odd[i] = A[i * 2 + 1];
    }
    even = FFT(even);
    odd = FFT(odd);
    for(int i = 0; i < N / 2; i++){
        odd[i] *= polar(1.0, 2 * M_PI * i / N);
        A[i] = even[i] + odd[i];
        A[N / 2 + i] = even[i] - odd[i];
    }
    return A;
}

C++ で高速フーリエ変換

FFT には定数倍高速化の技術がたくさんあり、ここから何倍か速くできるのですが、それは他の記事等に任せることにします。

逆高速フーリエ変換

高速フーリエ変換の ${\zeta }_N$ を ${\zeta }_N^{-1}$ に置き換えて、最後に $N$ で割ったものが逆高速フーリエ変換です。

高速フーリエ変換を利用した畳み込み

基本的には離散フーリエ変換を利用した畳み込みと同じで、$N$ が $2$ 冪になるように、$N=|A|+|B|-1$ とするのではなく、$|A|+|B|-1$ を $2$ 冪に切り上げたものを $N$ とすれば良いです。

おわり

🎉 これであなたも FFT の理解者！ 🎉

おまけ

巡回畳み込み

離散フーリエ変換で畳み込みをするとき、$|A|+|B|-1$ の長さを確保しますが、もしこれが足りなかったら、すなわち $N<|A|+|B|-1$ のときどうなるでしょうか？

$f_A({\zeta }_N^0),f_A({\zeta }_N^1),\dots ,f_A({\zeta }_N^{N-1})$ と $f_B({\zeta }_N^0),f_B({\zeta }_N^1),\dots ,f_B({\zeta }_N^{N-1})$ には特に影響はないので、$f_C({\zeta }_N^0),f_C({\zeta }_N^1),\dots ,f_C({\zeta }_N^{N-1})$ までは正しく計算できますね。しかし、$f_C$ は $N-1$ 次多項式として扱われ、長さ $N$ の $C$ が逆離散フーリエ変換によって出てきます。このときの $C$ はどのようなものになっているでしょうか？

例えば、$f_A(x)\times f_B(x)=x^N$ のとき、$f_C({\zeta }_N^0)=f_C({\zeta }_N^1)=\cdots =f_C({\zeta }_N^{N-1})=1$ で、これは $f_A(x)\times f_B(x)=1$ の場合と区別ができません。そのため、$C=(1,0,\dots ,0)$ となります。
また、$f_A(x)\times f_B(x)=x^{N+1}$ のとき、$f_C({\zeta }_N^0)={\zeta }_N^0,f_C({\zeta }_N^1)={\zeta }_N^1,\dots ,f_C({\zeta }_N^{N-1})={\zeta }_N^{N-1}$ で、これは $f_A(x)\times f_B(x)=x$ の場合と区別ができません。そのため、$C=(0,1,0,\dots ,0)$ となります。

このように、$N$ 次以上の項は、$x^N=1$ のルールにより $x^0$ の方から再び現れます。
このような場合の畳み込み

$$C_k=\sum _{i+j=k(\mathrm{m}\mathrm{o}\mathrm{d}N)}{A}_i{B}_j$$

を巡回畳み込みと言います。離散フーリエ変換は畳み込みをするための変換というより、巡回畳み込みをするための変換(十分な長さを確保することで通常の畳み込みになる)なんですね。

明日の担当は @liquid1224 さんです！楽しみ〜