【ePrint解説】Liberating TFHE 3/4 にて
いよいよ次回でこの Liberating TFHE の連載もラストです
その前にFFTやNTTに関する記事を独立で１本出したいと思います
と最後に書いたのですが，あれから４ヶ月も経ってしまった・・・
FFTやNTTに触れることになったので，これを機に今までやったことのあるフーリエ変換を復習しつつ，QDFTまで統一的に扱っていきます．
\mathbb{R} 上のFT
まずはオーソドックスなやつからです．
\underline{\mathrm{Def(} \mathbb{R} \mathrm{上の FT)}}
f \in L^1(\mathbb{R}) := \left \{ f \ | \ \displaystyle \int_{\mathbb{R}} |f(\xi)| d \xi < \infty \right \} のとき，\xi \in \mathbb{R} を変数とする f(\xi) について， \displaystyle \hat{f}(\xi) = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} f(x) e^{- \sqrt{-1} x \xi} d x と定めるとき，\hat{f} を \mathbb{R} 上の f のFourier 変換（FT）という．\hat{f} のことを \mathcal{F} f とも書く．
また，g \in L^1(\mathbb{R}) のとき，x \in \mathbb{R} を変数とする g(x) について， \displaystyle \check{g}(x) = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} g(\xi) e^{\sqrt{-1} x \xi} d \xi と定めるとき，\check{g} を \mathbb{R} 上の g の逆Fourier 変換（IFT）という．\check{g} のことを \mathcal{F}^{\ast} g とも書く．
ちょっと見慣れない方での定義になるかもしれませんが，量子離散フーリエ変換を見据えると，この形が自然かと
\underline{\mathrm{Exp}}
\displaystyle f(x) = e^{- x^2} とするとき，
\begin{align*}
\displaystyle
\mathcal{F} f(\xi)
& = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} f(x) e^{- \sqrt{-1} x \xi} d x \\
& = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} e^{- x^2} e^{- \sqrt{-1} x \xi} d x \\
& = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} e^{- x^2 - \sqrt{-1} \xi x} d x \\
& = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} e^{- (x + \frac{1}{2} \sqrt{-1} \xi)^2 - \frac{\xi^2}{4}} d x \\
& = \frac{1}{\sqrt{2 \pi}} e^{- \frac{\xi^2}{4}} \int_{\mathbb{R}} e^{- (x + \frac{1}{2} \sqrt{-1} \xi)^2} d x \\
& = \frac{1}{\sqrt{2 \pi}} e^{- \frac{\xi^2}{4}} \int_{\mathbb{R}} e^{- x^2} d x \\
& = \frac{1}{\sqrt{2 \pi}} e^{- \frac{\xi^2}{4}} \sqrt{\pi} \\
& = \frac{1}{\sqrt{2}} e^{- \frac{\xi^2}{4}}
\end{align*}
ゆえ，\displaystyle \mathcal{F} f(\xi) = \hat{f}(\xi) = \frac{1}{\sqrt{2}} e^{- \frac{\xi^2}{4}} である．また，
\begin{align*}
\displaystyle
\mathcal{F}^{\ast} (\mathcal{F} f)(x)
& = \frac{1}{\sqrt{2 \pi}} \int_{\mathbb{R}} g(\xi) e^{\sqrt{-1} x \xi} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} \int_{\mathbb{R}} e^{- \frac{\xi^2}{4}} e^{\sqrt{-1} x \xi} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} \int_{\mathbb{R}} e^{- \frac{1}{4} (\xi^2 + 4 \sqrt{-1} x \xi)} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} \int_{\mathbb{R}} e^{- \frac{1}{4} (\xi^2 + 2 \sqrt{-1} x)^2 - x^2} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} e^{- x^2} \int_{\mathbb{R}} e^{- \frac{1}{4} (\xi + 2 \sqrt{-1} x)^2} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} e^{- x^2} \int_{\mathbb{R}} e^{- \frac{1}{4} \xi^2} d \xi \\
& = \frac{1}{2 \sqrt{\pi}} e^{- x^2} \sqrt{4 \pi} \\
& = e^{- x^2}
\end{align*}
ゆえ，\mathcal{F}^{\ast} \mathcal{F} f = f となりました．
上の具体例から，一般に次が成り立ちます（証明は略）:
\underline{\mathrm{Prop}}
f, \hat{f} \in L^1(\mathbb{R}) に対して，\mathcal{F}^{\ast} \mathcal{F} f = f である．
\mathbb{C} 上の関数のDFT
いわゆる複素関数の離散フーリエ変換です．上の定義の離散化を考えます．
イメージとしては，(- \infty, \infty) の範囲を [0, 2 \pi) に制限して，等間隔に M 個の点を取ります．
\underline{\mathrm{Def(} \mathbb{C} \mathrm{上の関数の DFT)}}
f を \mathbb{C} 上の関数，M を正整数とする．\xi \in \mathbb{C} を変数とする f(\xi) について， \displaystyle \hat{f}(\xi) = \frac{1}{\sqrt{M}} \sum_{x = 0}^{M - 1} f(x) \exp \left(- \sqrt{-1} \frac{2 \pi \xi x}{M} \right) と定めるとき，\hat{f} を f の \mathbb{C} 上の離散フーリエ変換（DFT）という．
また，g を \mathbb{C} 上の関数，x \in \mathbb{C} を変数とする g(x) について，\displaystyle \check{g}(x) = \frac{1}{\sqrt{M}} \sum_{\xi = 0}^{M - 1} g(\xi) \exp \left(\sqrt{-1} \frac{2 \pi x \xi}{M} \right) と定めるとき，\check{g} を g の \mathbb{C} 上の逆離散フーリエ変換（IDFT）という．
\mathbb{C} 上の多項式のDFT
関数で定義できるなら，多項式でも定義できます（大雑把には関数は無限次数の多項式なので）
\underline{\mathrm{Def(} \mathbb{C} \mathrm{上の多項式の DFT)}}
f_0, f_1, \cdots, f_{M - 1} \in \mathbb{C}，M を正整数とする．\xi \in \mathbb{C} を変数とする f(\xi) = f_0 + f_1 \xi + \cdots + f_{M - 1} \xi^{M - 1} について，\displaystyle \hat{f}(\xi) = \frac{1}{\sqrt{M}} \sum_{x = 0}^{M - 1} f_x \exp \left(- \sqrt{-1} \frac{2 \pi \xi x}{M} \right) と定めるとき，\hat{f} を f の \mathbb{C} 上の離散フーリエ変換（DFT）という．
また，g_0, g_1, \cdots, g_{M - 1} \in \mathbb{C} で x \in \mathbb{C} を変数とする g(x) = g_0 + g_1 x + \cdots + g_{M - 1} x^{M - 1} について，\displaystyle \check{g}(x) = \frac{1}{\sqrt{M}} \sum_{\xi = 0}^{M - 1} g_{\xi} \exp \left(\sqrt{-1} \frac{2 \pi x \xi}{M} \right) と定めるとき，\check{g} を g の \mathbb{C} 上の逆離散フーリエ変換（IDFT）という．
ここで後の議論に備えて，\xi の部分と \exp の部分を次のように変形しておきます．
\underline{\mathrm{Def(} \mathbb{C} \mathrm{上の多項式の DFT)}}
f_0, f_1, \cdots, f_{M - 1} \in \mathbb{C}，M を正整数，\displaystyle \zeta_M = e^{\sqrt{-1} \frac{2 \pi}{M}} を1の原始 M 乗根とする．x \in \mathbb{C} を変数とする f(x) = f_0 + f_1 x + \cdots + f_{M - 1} x^{M - 1} について，\displaystyle \hat{f}(x) = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i (\zeta_M^{i})^x と定めるとき，\hat{f} を f の \mathbb{C} 上の離散フーリエ変換（DFT）という．
また，g_0, g_1, \cdots, g_{M - 1} \in \mathbb{C} で x \in \mathbb{C} を変数とする g(x) = g_0 + g_1 x + \cdots + g_{M - 1} x^{M - 1} について，\displaystyle \check{g}(x) = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} g_{i} (\zeta_M^{- i})^x と定めるとき，\check{g} を g の \mathbb{C} 上の逆離散フーリエ変換（IDFT）という．
多項式を配列とみなすなら，DFTというのは \mathbb{C}^M \rightarrow \mathbb{C}^M の関数で (f_0, f_1, \cdots, f_{M - 1}) を \displaystyle \left( \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i, \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^i, \cdots, \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{(M - 1)i} \right) に移すものとなります．
これは線型写像（さすがに明らか）なので，表現行列を考えることができ，次のような形になっています:
\begin{align*}
\begin{pmatrix}
1 & 1 & 1 & \cdots & 1 \\
1 & \zeta_M^{1 \cdot 1} & \zeta_M^{2 \cdot 1} & \cdots & \zeta_M^{(M - 1) \cdot 1} \\
1 & \zeta_M^{1 \cdot 2} & \zeta_M^{2 \cdot 2} & \cdots & \zeta_M^{(M - 1) \cdot 2} \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \zeta_M^{1 \cdot (M - 1)} & \zeta_M^{2 \cdot (M - 1)} & \cdots & \zeta_M^{(M - 1) \cdot (M - 1)} \\
\end{pmatrix}
\end{align*}
この行列は Vandermonde 行列ですので，（行列式が非零ゆえ）正則行列です．よって逆行列が存在し，その逆行列そのものがIDFTでの表現行列に対応します．
FFT
FFTは上記のDFTの計算を高速化したものです．
FFTにより，一般に m bit同士の掛け算の計算回数が O(m^2) から O(m \log m) へと落とすことができます．
どのように高速化するかを考えるために，DFTでの移り先をよく見る必要があります．
DFTでの移り先を観察する
(f_0, f_1, \cdots, f_{M - 1}) を移した \displaystyle \left( \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i, \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^i, \cdots, \frac{1}{\sqrt{M} } \sum_{i = 0}^{M - 1} f_i \zeta_M^{(M - 1)i} \right) の成分をいくつか調べてみます．まずは第２成分から（\zeta_M^M = 1 に注意）.
\begin{align*}
\displaystyle
& \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{2 \cdot i} \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{2 \cdot 1} + f_2 \zeta_M^{2 \cdot 2} + \cdots + f_{M/2 - 1} \zeta_M^{2 \cdot (M/2 - 1)} \\
& + f_{M/2} \zeta_M^{2 \cdot M/2} + f_{M/2 + 1} \zeta_M^{2 \cdot (M/2 + 1)} + \cdots + f_{M - 1} \zeta_M^{2 \cdot (M - 1)} ) \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{2} + f_2 \zeta_M^{4} + \cdots + f_{M/2 - 1} \zeta_M^{M - 2} \\
& + f_{M/2} + f_{M/2 + 1} \zeta_M^{2} + \cdots + f_{M - 1} \zeta_M^{M - 2} ) \\
& = \frac{1}{\sqrt{M}} ( (f_0 + f_{M/2}) + (f_1 + f_{M/2 + 1}) \zeta_M^{2} + \cdots + (f_{M/2 - 1} + f_{M - 1}) \zeta_M^{M - 2}) \\
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M/2 - 1} (f_i + f_{i + M/2}) \zeta_M^{2 \cdot i} \\
\end{align*}
良い感じに変形できました．続いて第３成分です．
\begin{align*}
\displaystyle
& \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{3 \cdot i} \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{3 \cdot 1} + f_2 \zeta_M^{3 \cdot 2} + \cdots + f_{M/2 - 1} \zeta_M^{3 \cdot (M/2 - 1)} \\
& + f_{M/2} \zeta_M^{3 \cdot M/2} + f_{M/2 + 1} \zeta_M^{3 \cdot (M/2 + 1)} + \cdots + f_{M - 1} \zeta_M^{3 \cdot (M - 1)} ) \\
\end{align*}
となり，一見するとどうにもなりそうもありませんが，\zeta_M^{M / 2} = - 1 に注意して，第２成分のときと同じような形に持っていけないかを考えると，
\begin{align*}
\displaystyle
& \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{3 \cdot i} \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{3 \cdot 1} + f_2 \zeta_M^{3 \cdot 2} + \cdots + f_{M/2 - 1} \zeta_M^{3 \cdot (M/2 - 1)} \\
& + f_{M/2} \zeta_M^{3 \cdot M/2} + f_{M/2 + 1} \zeta_M^{3 \cdot (M/2 + 1)} + \cdots + f_{M - 1} \zeta_M^{3 \cdot (M - 1)} ) \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{3 \cdot 1} + f_2 \zeta_M^{3 \cdot 2} + \cdots + f_{M/2 - 1} \zeta_M^{3 \cdot (M/2 - 1)} \\
& + f_{M/2} \zeta_M^{3 \cdot M/2} + f_{M/2 + 1} \zeta_M^{3 \cdot M/2} \cdot \zeta_M^{3} + \cdots + f_{M - 1} \zeta_M^{3 \cdot M / 2} \cdot \zeta_M^{3 \cdot (M / 2 - 1)} ) \\
& = \frac{1}{\sqrt{M}} ( f_0 + f_1 \zeta_M^{3 \cdot 1} + f_2 \zeta_M^{3 \cdot 2} + \cdots + f_{M/2 - 1} \zeta_M^{3 \cdot (M/2 - 1)} \\
& - f_{M/2} - f_{M/2 + 1} \zeta_M^{3} - \cdots - f_{M - 1} \zeta_M^{3 \cdot (M / 2 - 1)} ) \\
& = \frac{1}{\sqrt{M}} ( (f_0 - f_{M/2}) + (f_1 - f_{M/2 + 1}) \zeta_M^{3} + \cdots + (f_{M/2 - 1} - f_{M - 1}) \zeta_M^{3 \cdot (M / 2 - 1)}) \\
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M/2 - 1} (f_i - f_{i + M/2}) \zeta_M^{3 \cdot i} \\
\end{align*}
となり，良い感じに変形できました．
$\zeta_M^{M / 2} = - 1$
\zeta_M とは1の原始 M 乗根のことでした．つまり，\zeta_M^M = 1 であり，任意の 1 \leq i \leq M - 1 に対して \zeta_M^i \neq 1 を満たします．
ここで，X^M = 1 を変形した X^M - 1 = 0 を考えます．M は偶数ですから，(X^{M / 2} - 1)(X^{M / 2} + 1) = 0 となります．
\zeta_M は M 乗すると1なので X^M - 1 = 0 の解の1つです．よって，(X^{M / 2} - 1)(X^{M / 2} + 1) = 0 の解の1つなので，X^{M / 2} - 1 = 0 か X^{M / 2} + 1 = 0 の解になっています．
X = \zeta_M が X^{M / 2} - 1 = 0 の解であるとすると，\zeta_M^{M / 2} = 1 となり，「任意の 1 \leq i \leq M - 1 に対して \zeta_M^i \neq 1」に矛盾します．
よって，X = \zeta_M は X^{M / 2} + 1 = 0 の解であり，\zeta_M^{M / 2} = -1 を得ます．
第４成分以降も同じように変形でき，0 \leq n \leq M / 2 - 1 に対して，次のような予想が立てられます:
\begin{align*}
\displaystyle
\frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{2 n \cdot i}
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M / 2 - 1} (f_i + f_{i + M / 2}) \zeta_M^{2 n \cdot i} \\
\frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{(2 n + 1) \cdot i}
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M / 2 - 1} (f_i - f_{i + M / 2}) \zeta_M^{(2 n + 1) \cdot i}
\end{align*}
実際，この予想は正しいです．証明は上での具体例において，2 \mapsto 2 n, 3 \mapsto 2 n + 1 と変更してください．
ここで１つ注目する点として，\zeta_M^2 = \zeta_{M / 2} が成り立ちます．
上の性質が成り立つ理由
色々な証明は与えられるのですが，最も簡単なものは定義に戻るやり方かなと．
\zeta_M = \exp^{\sqrt{-1} 2 \pi / M} なので，\zeta_M^2 = \exp^{\sqrt{-1} 2 \cdot 2 \pi / M} = \exp^{\sqrt{-1} 2 \pi / (M / 2)} = \zeta_{M / 2} です．
他の証明の方針としては
- 円の M 等分点を考える
-
X = \zeta_M^2 も X^M = 1 の解であることに注目する
などがあったりします．
すると上での予想（定理）は次のように変形できます．
\begin{align*}
\displaystyle
\frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{2 n \cdot i}
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M / 2 - 1} (f_i + f_{i + M / 2}) \zeta_{M / 2}^{n \cdot i} \\
\frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \zeta_M^{(2 n + 1) \cdot i}
& = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M / 2 - 1} (f_i - f_{i + M / 2}) \zeta_M^{i} \cdot \zeta_{M / 2}^{n \cdot i}
\end{align*}
バタフライ演算
上記の定理の係数を無視した部分を考えると，例えば，\sum_{i = 0}^{M - 1} f_i \zeta_M^{2 n \cdot i} が \sum_{i = 0}^{M / 2 - 1} (f_i + f_{i + M / 2}) \zeta_{M / 2}^{n \cdot i} に移っていることがわかります．
つまり，(f_0, f_1, \cdots, f_{M - 1}) という長さ M の配列の 2 n 番目の成分をDFTで移した先と (f_0 + f_{M / 2}, f_1 + f_{M / 2 + 1}, \cdots, f_{M / 2 - 1} + f_{M - 1}) という長さ M / 2 の配列の n 番目の成分をDFTで移した先が同じ，ということです．
同様に(f_0, f_1, \cdots, f_{M - 1}) という長さ M の配列の 2 n + 1 番目の成分をDFTで移した先と ((f_0 - f_{M / 2}), (f_1 - f_{M / 2 + 1}) \zeta_M, \cdots, (f_{M / 2 - 1} - f_{M - 1}) \zeta_M^{M / 2 - 1}) という長さ M / 2 の配列の n 番目の成分をDFTで移した先が同じです．
この性質を利用したのがバタフライ演算です．
例えば，M = 8 では，
-
(f_0, f_1, f_2, f_3, f_4, f_5, f_6, f_7) をDFTで移した長さ8の配列
と
-
(f_0 + f_4, f_1 + f_5, f_2 + f_6, f_3 + f_7) をDFTで移した長さ4の配列と(f_0 - f_4, (f_1 - f_5) \zeta_8, (f_2 - f_6) \zeta_8^2, (f_3 - f_7) \zeta_8^3) をDFTで移した長さ4の配列を結合した長さ8の配列
が「大体」同じになります．
＊ここでいう「大体」は順序を除いて，という意味です（詳しくは下図を）
イメージは次のようになります:
上図でDFTの添字は長さ8の配列に対するDFTか長さ４の配列に対するDFTかを表しています．
Cooley-Tukeyの発想
バタフライ演算で何をやっているのかを振り返ると，長さ8の配列に対するDFTを1回行う操作を長さ4の配列に対するDFTを2回行う操作に変えられるということを主張しています．
長さ M の配列に対するDFTでは，M \times M サイズの表現行列をその配列に掛けるのですから，それが (M / 2) \times (M / 2) サイズの表現行列を掛けることを２回行なっても，後者の方が計算量が少ないメリットがあります．
「長さ8の配列に対するDFTを1回行う操作を長さ4の配列に対するDFTを2回行う操作に変えられる」ということは「長さ４の配列に対するDFTを1回行う操作を長さ２の配列に対するDFTを2回行う操作に変えられる」ということで，「長さ2の配列に対するDFTを1回行う操作を長さ1の配列に対するDFTを2回行う操作に変えられる」ということになります．
長さ1の配列に対するDFTとは結局何もしていないのと同じですので，上のように配列を組み替える操作を繰り返すことで結果的に長さ8の配列に対するDFTを得ることができます．
これを定式化したものが Cooley-Tukey のアルゴリズムと呼ばれるもの（の一部分）です．
長さ8のときのイメージ図は上のようになります（手書きですみません・・・）
bit-reverse
さて，Cooley-Tukey の概要を掴んだところで１つ考えるべきことがあります．
(f_0, f_1, f_2, f_3, f_4, f_5, f_6, f_7) をDFTによって得られる配列は (\hat{f_0}, \hat{f_1}, \hat{f_2}, \hat{f_3}, \hat{f_4}, \hat{f_5}, \hat{f_6}, \hat{f_7}) のようになってほしいのですが，
最終的に得られる配列が (\hat{f_0}, \hat{f_4}, \hat{f_2}, \hat{f_6}, \hat{f_1}, \hat{f_5}, \hat{f_3}, \hat{f_7}) と順序があべこべになっています．
ここで唐突ですが，その配列の添字を2進展開した配列を考えると，(000, 100, 010, 110, 001, 101, 011, 111) となっています．
そして，各成分に対して bit-reverse という操作を行います．これは x_1 x_2 \cdots x_n という文字列に対して x_n \cdots x_2 x_1 のように，元の文字列を「逆順」にする操作のことです．
すると，(000, 001, 010, 011, 100, 101, 110, 111) となり，10進数に戻すと (0, 1, 2, 3, 4, 5, 6, 7) となっています．
何をしているのかというと，(\hat{f_0}, \hat{f_4}, \hat{f_2}, \hat{f_6}, \hat{f_1}, \hat{f_5}, \hat{f_3}, \hat{f_7}) という配列のインデックスに対して bit-reverse を行いその順序に従って並べた配列は (\hat{f_0}, \hat{f_1}, \hat{f_2}, \hat{f_3}, \hat{f_4}, \hat{f_5}, \hat{f_6}, \hat{f_7}) になる，ということです．
つまり，(\hat{f_0}, \hat{f_4}, \hat{f_2}, \hat{f_6}, \hat{f_1}, \hat{f_5}, \hat{f_3}, \hat{f_7}) を (\hat{f_0}, \hat{f_1}, \hat{f_2}, \hat{f_3}, \hat{f_4}, \hat{f_5}, \hat{f_6}, \hat{f_7}) にするには，並べ替えの操作（置換）が必要だが，それは bit-reverse という操作で１発で行える，ということになります．
この性質が一般に成り立つかを調べてみます．つまり，
-
(f_0, f_1, \cdots, f_{M - 1}) なる配列に対して，再帰的にバタフライ演算を作用させることで配列 (\hat{f}_0, \hat{f}_1, \cdots, \hat{f}_{M - 1}) を得る
-
(\hat{f}_0, \hat{f}_1, \cdots, \hat{f}_{M - 1}) の添字に対して bit-reverse による置換を作用させた配列は (f_0, f_1, \cdots, f_{M - 1}) にDFTを作用させた配列と一致する
ことを考えればOKです．そして示すべきことは
1によって得られる配列のインデックスが (0..00, 0..01, ..., 1..11) という配列の各成分に bit-reverse を作用させて10進数にしたものと一致する
です．
＊bit-reverse を2回行うと元の文字列に戻ることに注意してください
これは次のような対応で示すことができます．
先ほどの長さ8でのイメージ図において，各配列は長さ0以上の \{ 0, 1 \} 上の文字列 x を持っているものとします．このとき，
- 初期配列 (f_0. f_1, \cdots, f_7) は空列を持っている
- 文字列 x を持っている配列が
- 赤い点線で囲まれる配列に移る場合，その文字列の先頭に0を追加する
- 青い点線で囲まれる配列に移る場合，その文字列の先頭に1を追加する
という操作を考えます．
先ほどの図で確認してみると
このようになっていて，確かに最終的に得られる長さ１の各配列が持っている文字列は 000, 100, 010, 110, 001, 101, 011, 111 で，
各成分を bit-reverse すると 000, 001, 010, 011, 100, 101, 110, 111 となっていて確かにOKです．
Cooley-Tukey のアルゴリズム
以上をまとめて，Cooley-Tukey のアルゴリズムを与えます．
まずは疑似コード版です．
＊アルゴリズムのインデックスは気にしないでください・・・
上記は再帰を用いて書きましたが，直接ループで回すこともできます．（FFTの部分のみ）そのコードをpythonで書いてみます．
DFT_size = M
list_num = M / DFT_size
while DFT_size >= 1:
for list_index in range(list_num):
list_start_index = list_index * DFT_size
for i in range(DFT_size / 2):
left_index = list_start_index + i
right_index = list_start_index + i + DFT_size / 2
F[left_index], F[right_index] = F[left_index] + F[right_index], (F[left_index] - F[right_index]) * pow(zeta_M, element_index)
bin_index_list[left_index] = '0' + bin_index_list[left_index]
bin_index_list[right_index] = '1' + bin_index_list[right_index]
DFT_size /= 2
list_num *= 2
環上のDFT
\mathbb{C} 上の多項式のDFTの定義において，\displaystyle \zeta_M^{- x \xi} に注目します．ここで，一般の環 R において「1の原始 M 乗根」なるものを考える（厳密にはこのような「1の原始 M 乗根」が存在するような環を考える）ことにより，環上のDFTへと一般化できます:
\underline{\mathrm{Def(環上のDFT)}}
R を環，M を2べき，f_0, f_1, \cdots, f_{M - 1} \in R とする．\displaystyle \zeta_M \in R を1の原始 M 乗根として，\xi \in R を変数とする \displaystyle f(\xi) = \sum_{i = 0}^{M - 1} f_i \xi^i = f_0 + f_1 \xi + \cdots + f_{M - 1} \xi^{M - 1} \in R とする f(\xi) について， \displaystyle \hat{f}(\xi) = \frac{1}{\sqrt{M}} \sum_{i = 0}^{M - 1} f_i \xi^i と定めるとき，\hat{f} を f の環上のDFTという．
逆変換も同様に定義できます．
NTT
一般の環において，1の原始 M 乗根が存在するとは限らなかったのですが，これは体で考えることで解決されます．特に実用上では有限体を考えれば良く，p を奇素数とするとき，\mathbb{Z} / p \mathbb{Z} で考えることで，1の原始 M 乗根が存在するかどうか問題を解決できます．
この\mathbb{Z} / p \mathbb{Z} 上でのDFTをNTT（数論変換）と言います．
QDFT
QDFT は量子版の離散フーリエ変換です．
\underline{\mathrm{Def(QDFT)}}
m を正整数として M = 2^m とする．\mathcal{H}^{\otimes m} の正規直交基底を1つとり，それを (| 0 \rangle, \cdots | M - 1 \rangle) と表す．| f \rangle \in \mathcal{H}^{\otimes m} に対して，| f \rangle = \sum_{i = 0}^{M - 1} f_i | i \rangle とする．ここで \sum_{i = 0}^{M - 1} || f_i ||^2 = 1 である．この | f \rangle に対して，\displaystyle | \hat{f} \rangle = \frac{1}{\sqrt{M}} \sum_{j = 0}^{M - 1} \sum_{i = 0}^{M - 1} f_i \exp \left( \sqrt{-1} \frac{2 \pi i j}{M} \right) なる量子状態に移すユニタリ変換を量子（離散）フーリエ変換という．
逆行列も同様に定義されます．量子フーリエ変換の話は Shorの素因数分解アルゴリズム でも書いているので，もしよろしければそちらも．
今回の内容はここまでです．ここまでご覧になってくださった方々ありがとうございます！
Discussion