Understanding Fast Fourier Transform from scratch — to solve Polynomial Multiplication.

Motivation

TFHEの実装眺めていたらFFTが使われていた.

整数多項式で高速フーリエ変換(FFT)

2多項式の積

$A(x)=a_0+a_{1}x+a_2{x}^2\cdots a_n{x}^n$

$B(x)=b_0+b_{1}x+b_2{x}^2\cdots b_n{x}^n$

のとき

$C(x)=A(x)\ast B(x)=c_0+c_{1}x+\cdots +c_{2n}{x}^{2n}$

ただし $c_i=a_0{b}_i+a_1{b}_{i-1}\cdots a_i{b}_0$.

Example

E.g., multiply (x^2 + 2x + 3)(2x^2 + 5) = 

			              2x^2 + 0x + 5
			              1x^2 + 2x + 3
                        -------------
                        6      0   15
                 4      0      10
           2     0      5
         ----------------------------
          2x^4 + 4x^3 + 11x^2 +10x+15

これを $A$ と$B$ の畳み込み(Convolution) といい、 $A\otimes B$ と表す. 素朴な計算量は $O(n^2)$

フーリエ変換と畳込み理論

関数 $f$, $g$ の畳み込みはフーリエ変換 $\zeta$ とするとpoint-wise multiplicationに書き換えられる(後述).

$$\zeta (f\otimes g)=\zeta (f).\zeta (g)$$

多項式(polynomials)の表現

多項式は係数表現とpoint-value representation があります.

point-value form

$n-1$ 次多項式は最低 $n$ 個の異なる点として表すことができる.

$A(x)=a_0+a_{1}x+a_2{x}^2\cdots x^n$ は以下のように表される.

• n個のペア $\{(x_0,y_0),(x_1,y_1),\cdots (x_n,y_n)\}$
• s.t ${}^{\forall }i\ne j,x_i\ne x_j$
• s.t. $y_k=A(x_k)$

point-value form では $C(x)=A(x)B(x)$ は $C(x_k)=A(x_k).B(x_k)$ for any point $x_k$

$A,B$ が $n$次元以下(degree-bound ‘n’)の時 $C$ は $2n$ 次となるので$A,B$ の点を $2n$ 個に増やす.

Example

$A(x)=x^3-2x+1$, $B(x)=x^3+x^2+1$

$x_k=(-3,-2,-1,0,1,2,3)$ の7個の係数が必要

$A:(-3,-17),(-2,-3),(-1,1),(0,1),(1,0),(2,5),(3,22)$
$B:(-3,-20),(-2,-3),(-1,2),(0,1),(1,3),(2,13),(3,37)$
$C:(-3,340),(-2,9),(-1,2),(0,1),(1,0),(2,65),(3,814)$

これが point-wise multiplication, 計算量は $\Theta (n)$

$\Theta (n)$ とは

計算量の上界 $O(n)$ と下界 $\Omega (n)$ が一致する時, e.g. $2n^2+5n+1000=\Theta (n^2)$

逆に点集合から多項式を求めることを 補完(iterpolation) という.

係数表現とpoint-value representationの高速な変換が必要. FFTの出番.

アルゴリズムの概観

$m=2n-1,\mathrm{degree}\mathrm{of}\mathrm{C}\le \mathrm{m}$

1. $x_0\cdots x_{m-1}$ から $m$ 点を賢く選ぶ
2. $A(x_0),\cdots ,A(x_{m-1})$ を計算
3. $B$ も同様に
4. $A(x)\ast B(x)$ で $C(x_0)\cdots C(x_{m-1})$ を計算
5. Cから係数を補完

• 手順2, 3は $O(n^2)$ かかりそうだがFFTで高速化可能
• 賢く: root of unity から $m$ 点選ばないとだめ

迂回: 1の冪根(Roots of unity)

複数回かけて $1$ なる数のこと. 1以外の Roots of unity は複素数

複素数

$$r{e}^{i\theta }=(r\cos \theta ,r\sin \theta )$$ $$\omega =\cos (\frac{2\pi }{m})+i\sin (\frac{2\pi }{m})$$

$({\omega }_n^0,{\omega }_n^1\cdots ,\times )\simeq ({\mathbb{Z}}_n,+)modn$

FFTは以下の性質を使う

• ${\omega }_n^n=1$
• ${\omega }_n^n+k={\omega }_n^k$
• ${\omega }_n^{n/2}=-1$
• ${\omega }_n^{n/2}+k=-{\omega }_n^k$

離散フーリエ変換(DFT)

目的は $m$ 点($m=2^k$)の $A(x),\mathrm{degree}<\mathrm{m}$ を $O(m\log m)$ で評価すること

クーリー–テューキー型FFTアルゴリズム

Examples

$m=8$,
$A(x)=a_0+\cdots a_7{x}^7$, vec: $A=[a_0,\cdots a_7]$

奇数と偶数で分ける.

$A_{\mathrm{even}}=[a_0,a_2,a_4,a_6]$,
$A_{\mathrm{odd}}=[a_1,a_3,a_5,a_7]$

多項式表現では

$A_{\mathrm{even}}=a_0+a_{2}x+a_4{x}^2+a_6{x}^3$
$A_{\mathrm{odd}}=a_1+a_{3}x+a_5{x}^2+a_7{x}^3$

次元は $<m/2$, でA(x), A(-x) は次のようにかける.

$A(x)=A_{\mathrm{even}}(x^2)+x{A}_{\mathrm{odd}}(x^2)$
$A(-x)=A_{\mathrm{even}}(x^2)-x{A}_{\mathrm{odd}}(x^2)$

1つの $m$ 次多項式を計算する代わりに $m/2$ 次多項式を2つ計算することになる.(分割統治法)

再帰的にできるので $O(m\log m)$

Example

$A(x)=2+2x+3x^2+4x^3$

のとき
$A(x)=A_{\mathrm{even}}(x^2)+x\cdot A_{\mathrm{odd}}(x^2)$
を用いて
$A_{\mathrm{even}}(x)=3+3x$
$A_{\mathrm{odd}}(x)=2+4x$
より

• $A({\omega }_4^0)=A(1)=A_{\mathrm{even}}(1)+1\cdot A_{\mathrm{odd}}(1)=3+3+2+4=12$
• $A({\omega }_4^1)=A(i)=A_{\mathrm{even}}(i)+i\cdot A_{\mathrm{odd}}(i)=3-3+2i+4i=6i$
• $A({\omega }_4^2)=A(-1)=A_{\mathrm{even}}(-1)-1\cdot A_{\mathrm{odd}}(-1)=3+3-2-4=0$
• $A({\omega }_4^3)=A(-i)=A_{\mathrm{even}}(-i)-i\cdot A_{\mathrm{odd}}(-i)=3-3-2i+4i=2i$