# 快速傅立葉變換FFT
本篇僅談及FFT在多項式乘法及大數乘法的用途，不含信號轉換相關部分，有興趣可自行找其他資料閱讀
## 多項式乘法
假設多項式
$A(x)=2x^2-x+1,\ B(x)=2x-3$
則計算$A\times B$時我們會這樣做：
$\begin{align}
&\ \ \ \ \ \ \ \ \ \ \ 2x^2-x+1\\
\times)&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x+3\\
& --------\\
& \ \ \ \ \ \ \ \ \ \ 6x^2-3x+3\\
& \ 2x^3-x^2+\ \ x\\
& --------\\
& 2x^3+5x^2-2x+3
\end{align}$
若$A$和$B$的次數中較大者為$N$，則此計算方法的時間複雜度為$O(N^2)$
**好慢**
所以我們要加速
## 點值表示法
在表達一個多項式時，除了用習慣的$A(x)=2x^2-x+1$稱為"係數表示法"以外，我們發現如果對於一個$k-1$次多項式
$A(x)=a_0+a_1x+\dots+a_{k-2}x^{k-2}+a_{k-1}x^{k-1}$
選定$k$個點$x_0\sim x_{k-1}$便可算出$a_0\sim a_{k-1}$，亦即決定唯一的多項式$A(x)$
所以我們也可以將多項式轉為表達成數對
$A=\{(x_0,\ y_0),\ (x_1,\ y_1),\dots,\ (x_{k-1},\ y_{k-1}),\ (x_{k}, y_{k})\}$
這種表示法稱為"點表示法"或"點值表示法"
## 點值表示法的乘法
因為$C(x)=A(x)B(x)$，故$C(t)=A(t)B(t)$，即要用點值表示法將兩個多項式相乘時，只要固定$x_i$，把$y_i$對應相乘即可
一樣用上面的例子
$A(x)=1-x+2x^2,\ B(x)=-3+2x$
因為較高次為$A$的$2$次，因此改為點值表示法時需$3$個點，為方便計算取$0、1、2$
故
$A(x)=\{(0,\ 1),\ (1,\ 2),\ (2,\ 7)\}$、$B(x)=\{(0,\ -3),\ (1,\ -1),\ (2,\ 1)\}$
設$A(x)\times B(x)=C(x)$，則$C(x)=\{(0,\ -3),\ (1,\ -2),\ (2,\ 7)\}$
但我們發現$C$應該至少要有$3$次項，也就是至少需$4$個點來表示
更進一步來說，兩個多項式相乘，他們的次數差不多會加倍
所以真正計算兩個次數為$n$的多項式乘法時，應該各取$2n$個點（取太多是沒關係的，但不能少取）
回到上面的例子，應取$2\times 2=4$個點才能計算：
$A(x)=\{(0,\ 1),\ (1,\ 2),\ (2,\ 7),\ (3,\ 16)\}$、$B(x)=\{(0,\ -3),\ (1,\ -1),\ (2,\ 1),\ (3,\ 3)\}$
$C(x)=\{(0,\ -3),\ (1,\ -2),\ (2,\ 7),\ (3,\ 48)\}$
到這裡我們統整一下利用點值表示法做多項式乘法的過程
1. 將係數表示轉為點值表示，對於兩個係數為$n$的多項式要取$2n$個點，因為每個點都須跑過整個多項式去計算出答案，故時間複雜度$O(n^2)$
2. 點乘：將對應的點相乘，時間複雜度$O(n)$
3. 內插：用某種方法把得到點值表示的答案多項式轉回係數表示
這樣的話總時間複雜度$O(n^2)$，沒有達到加速的效果
而步驟2是沒辦法加速的，我們能改進的只有步驟1和3
這時就可以用到快速傅立葉變換FFT了，他的特點在於取那個$x_i$的值，並非整數甚至實數，而是複數。因此這時我們要先了解一些複數的東西 (當然已經會的可以跳過)
## 複數及複平面
### 先備知識
定義$i=\sqrt{-1}$
而複數就是表達為$z=a+bi$的數，其中$a、b\in\mathbb{R}$
### 複平面
如果將座標平面的橫軸令為實數軸、縱軸為虛數軸，那麼這個平面稱為複平面，且所有複數$z=a+bi$均可在這上面找到對應的點$Z(a,\ b)$
而原點$O$與$Z$形成向量$\overrightarrow{OZ}$，$z$的模定義為$|z|=|\overrightarrow{OZ}|=\sqrt{a^2+b^2}$、與橫軸形成的有向角稱為輻角，表為$\arg(z)$
因為有同界角，故$\arg(z)$有無限個，其中只有一個主輻角$\text{Arg}(z)$，滿足$0\le \text{Arg}(z)\le 2\pi$（有些地方會寫$-\pi\le \text{Arg}(z)\le \pi$），可知
$\arg(z)=\text{Arg}(z)+2k\pi,\ k\in \mathbb{Z}$
根據三角函數，若複數$z=a+bi$模為$|z|=r$、主輻角為$\theta$，則$Z(r\cos\theta,\ r\sin\theta)$，$z=r(\cos\theta+i\sin\theta)$
若$z$在單位圓上，即$r=1$，則$z=\cos\theta+i\sin\theta$
以下所提及的複數都在單位圓上
### 複數的次方
:::info
棣美弗定理：$(\cos\theta+i\sin\theta)^n=\cos (n\theta)+i\sin(n\theta)$
:::spoiler 證明
(一)
利用數學歸納法，先證明$n=0$及正整數，再推至負整數
1.若$n=0$，左式$=$右式$=1$
2.若$n=1$，左式$=$右式$=\cos\theta+i\sin\theta$
3.若$n=k$時成立，$(\cos\theta+i\sin\theta)^k=\cos(k\theta)+i\sin(k\theta)$
則當$n=k+1$時，
$\begin{align}
(\cos\theta+i\sin\theta)^{k+1}&=(\cos\theta+i\sin\theta)^k\cdot(\cos\theta+i\sin\theta)\\
&=(\cos(n\theta)+i\sin(n\theta))(\cos\theta+i\sin\theta)\\
&=\cos(k\theta)\cos\theta-\sin(k\theta)\sin\theta+i(\sin\theta\cos(k\theta)+\sin(k\theta)+\cos\theta)\\
&\overset{A}{=}\cos(k\theta+\theta)+i\sin(k\theta+\theta)\\
&=\cos((k+1)\theta)+i\sin((k+1)\theta)
\end{align}$
其中$A$處使用和角公式
故依數學歸納法，知：對於$\forall n\in\mathbb{Z}$且$n\ge 0$成立
又，由恆等式
$(\cos(n\theta)+i\sin(n\theta))(\cos(-n\theta)+i\sin(-n\theta))=\cos0+i\sin0=1$
即
$(\cos\theta+i\sin\theta)^n(\cos(-n\theta)+i\sin(-n\theta))=1$
可知，
$(\cos(-n\theta)+i\sin(-n\theta))=(\cos\theta+i\sin\theta)^{-n}$
故對於$n\in\mathbb{Z}$皆成立
(二)
根據歐拉公式：
$\cos \theta+i\sin \theta=e^{i\theta}$
故
$\begin{align}(\cos\theta+i\sin\theta)^n&=e^{(i\theta)n}\\
&=e^{i(n\theta)}\\
&=\cos(n\theta)+i\sin(n\theta)
\end{align}$
:::
利用棣美弗定理可以很簡單地計算方程式$w^n=z$如下：
設$z=\cos\theta+i\sin\theta$，$w=\cos\varphi+i\sin\varphi$
因為$w^n=z$，根據棣美弗定理：
$w^n=(\cos\varphi+i\sin\varphi)^n=\cos (n\varphi)+i\sin(n\varphi)=\cos\theta+i\sin\theta$
即
$n\varphi=\theta+2k\pi,\ \varphi=\displaystyle\frac{\theta+2k\pi}{n},\ k\in\mathbb{Z}$
因此有公式：
:::info
$w=\cos\displaystyle\frac{\theta+2k\pi}{n}+i\sin\displaystyle\frac{\theta+2k\pi}{n},\ k\in\mathbb{Z}$
:::
在$[0,\ 2\pi]$的範圍內，取$k=0、1、\dots、n-1$，則$w$有$n$個解
:::warning
:::spoiler 例如：求1的三次複數根
解：即求$w^3=1$的解
因為$1$就是實數軸上的$1$，設$z=\cos\theta+i\sin\theta$，$\theta=0$
故$w=\cos\displaystyle\frac{2k\pi}{n}+i\sin\displaystyle\frac{2k\pi}{n},\ k\in\mathbb{Z}$
在$[0,\ 2\pi]$的範圍內，$w$有$3$個解：
$w_0=\cos 0+i\sin 0=1$
$wa_1=\cos\displaystyle\frac{2\pi}{3}+i\sin\displaystyle\frac{2\pi}{3}=\frac{-1+\sqrt3 i}{2}$
$w_2=\cos\displaystyle\frac{4\pi}{3}+i\sin\displaystyle\frac{4\pi}{3}=\frac{-1-\sqrt3 i}{2}$
:::
畫成圖的話就是將單位圓平分成$n$塊。這也可以解釋為什麼只有$n$個不重複的解
3次根的圖：
## 取值
現在回到FFT
上面說到對於兩個次數為$m$的多項式$A(x),\ B(x)$要相乘需要取$2m$個點
以下令$n=2m$
那麼我們要取的這$n$點就是$w^n=1$的$n$個解。記為$w^0_n,\ w^1_n,\ \dots\ ,w^{n-1}_n$
以下以$A$舉例，$B$也是一樣的做法
$A(x)=a_0+a_1x+a_2x^2+\dots+a_{n-1}x^{n-1}=\sum\limits_{i=0}^{n-1}a_ix^i$
所以我們需要計算出$n$個值$y_0\sim y_{n-1}$，其中$y_k=A(w^k_n)=\sum\limits_{i=0}^{n-1}a_iw^{ki}_n$
如果將各$a_i$表達為$n$維向量$\vec a=<a_0,\ a_1,\ \dots\ ,\ a_{n-1}>$、計算出的$n-1$個$y$值也寫成向量$y=<y_0,\ y_1,\ \dots\ ,\ y_{n-1}>$，那麼稱$y$為$a$的離散傅立葉變換DFT(Discrete Fourier Transform)
我們要找一種能快速計算出$y$的方法，這種方法就稱為快速傅立葉變換FFT(Fast Fourier Transform)
FFT利用了分治法，時間複雜度為$O(n\log n)$
## 分治
首先，FFT使用的是分治法。他的做法是將一個次數為$2T$的多項式分成兩個次數為$T$的小多項式，故我們必須確保$\deg(A(x))$為$2$的次方。
又，更上面說到我們在做乘法時應該將次數加倍
所以要正式套用FFT前，我們應該要將$A(x)$的次數先擴展到比他大、且離他最近的$2$的次方，再$\times 2$。而多的項的係數皆為$0$。
例如次數為$5$，那就先變成$8$，再加倍成$16$
現在我們假設$\deg(A(x))=n$為$2$的次冪
FFT對$A(x)$定義了兩個新的多項式$A^{[0]}(x),\ A^{[1]}(x)$，其中
$A^{[0]}(x)=a_0+a_2x+a_4x^2+\dots+a_{n-2}x^{\frac{n}{2}-1}$
$A^{[1]}(x)=a_1+a_3x+a_5x^2+\dots+a_{n-1}x^{\frac{n}{2}-1}$
這時把$k^2$代入$A^{[0]}(x)$及$A^{[1]}(x)$，得到
$A^{[0]}(k^2)=a_0+a_2k^2+a_4k^4+\dots+a_{n-2}k^{n-2}$
$A^{[1]}(k^2)=a_1+a_3k^2+a_5k^4+\dots+a_{n-1}k^{n-2}$
得到$A(k)=A^{[0]}(k^2)+k\cdot A^{[1]}(k^2)$
這時，我們要計算$A(w^0_n),\ A(w^1_n),\ \dots\ ,\ A(w^{n-1}_n)$就簡化成計算$A^{[0]}(w^0_n),\ A^{[0]}(w^1_n),\ \dots\ ,\ A^{[0]}(w^{n-1}_n)$及
$A^{[1]}(w^0_n),\ A^{[1]}(w^1_n),\ \dots\ ,\ A^{[1]}(w^{n-1}_n)$。
雖然好像變多了，但就算使用傳統的$O(n^2)$計算每個多項式也可以將總複雜度由$O(n^2)$減為$O((\frac n2)^2+(\frac n2)^2)=O(\frac{n^2}{2})$
舉例$n=4$就是
$y_0=A(w^0_4)=A^{[0]}((w^0_4)^2)+w^0_4\cdot A^{[1]}((w^0_4)^2)$
$y_1=A(w^1_4)=A^{[0]}((w^1_4)^2)+w^1_4\cdot A^{[1]}((w^1_4)^2)$
$y_2=A(w^2_4)=A^{[0]}((w^2_4)^2)+w^2_4\cdot A^{[1]}((w^2_4)^2)$
$y_3=A(w^3_4)=A^{[0]}((w^3_4)^2)+w^3_4\cdot A^{[1]}((w^3_4)^2)$
值得注意的是$(w^k_n)^2=w^{2k}_n$
:::success
:::spoiler 為什麼
由$w$的公式解：$\displaystyle w^k_n=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}$
由棣美弗定理，
$\displaystyle(w^k_n)^a=\cos\frac{2ka\pi}{n}+i\sin\frac{2ka\pi}{n}=w^{ka}_n$
:::
這時我們引入**相消定理**
:::info
**相消定理**：
若$d\in\mathbb{N}$、$k\ge 0$、$n\ge 0$，則$w^{dk}_{dn}=w^k_n$
:::spoiler 證明
$\displaystyle w^{dk}_{dn}=\cos\frac{2dk\pi}{dn}+i\sin\frac{2dk\pi}{dn}=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}=w^k_n$
:::
## 整理
$1.$ $(w^k_n)^2=w^{2k}_n$
根據**相消定理**，可以得到推論：
$2.$ $\forall$偶數$v$，$\displaystyle w^{\frac{v}{2}}_v=w^1_2$(上下同除以$\frac{v}{2})=-1(w^2=1$的第二個解)
由$2.$延伸可知
$3.$ $w^{k+\frac n2}_n=w^k_n\cdot w^{\frac2n}_n=-w^k_n$
由$3.$延伸可知
$4.$ $(w^{k+\frac n2}_n)^2=(-w^k_n)^2=(w^k_n)^2$
再者，回想上面$w^n=1$的公式解，可知
$5.$ $w^k_n=w^{k\mod n}_n$
最終得到結論：
對於$p<\displaystyle\frac{n}{2}$，寫$k=p$，得
$A(p)=A(w^k_n)=A^{[0]}(w^{2k}_n)+w^k_n\cdot A^{[1]}(w^{2k}_n)=A^{[0]}(w^{k\ mod\ \frac n2}_\frac n2)+w^k_n\cdot A^{[1]}(w^{k\ mod\ \frac n2}_\frac n2)$
對於$p\ge\displaystyle\frac{n}{2}$，寫$k=p-\displaystyle\frac{n}{2}$，得
$A(p)=A(w^{k+\frac n2}_n)=A^{[0]}(w^{2k}_n)-w^k_n\cdot A^{[1]}(w^{2k}_n)=A^{[0]}(w^{k\ mod\ \frac n2}_\frac n2)-w^k_n\cdot A^{[1]}(w^{k\ mod\ \frac n2}_\frac n2)$
回到剛剛$n=4$的例子就是
$y_0=A(w^0_4)=A^{[0]}(w^0_4)+w^0_4\cdot A^{[1]}(w^0_4)=A^{[0]}(w^0_2)+w^0_4\cdot A^{[1]}(w^0_2)$
$y_1=A(w^1_4)=A^{[0]}(w^2_4)+w^1_4\cdot A^{[1]}(w^2_4)=A^{[0]}(w^1_2)+w^1_4\cdot A^{[1]}(w^1_2)$
$y_2=A(w^2_4)=A^{[0]}(w^4_4)+w^2_4\cdot A^{[1]}(w^4_4)=A^{[0]}(w^0_2)-w^0_4\cdot A^{[1]}(w^0_2)$
$y_3=A(w^3_4)=A^{[0]}(w^6_4)+w^3_4\cdot A^{[1]}(w^6_4)=A^{[0]}(w^1_2)-w^1_4\cdot A^{[1]}(w^1_2)$
理解之後就可以寫出虛擬碼：
```python=
def FFT( A[] ) :
n := A.length()
if n == 1:
return A
A0 := { A[0], A[2], ... , A[n - 2] }
A1 := { A[1], A[3], ... , A[n - 1] }
y0 := FFT( A0 )
y1 := FFT( A1 )
y := {}
ang := 2 * pi / n
w := 1
wn := cos(ang) + i * sin(ang)
for i in [0, n / 2 - 1] do:
y[i] := y0[i] + w * y1[i]
y[i + n / 2] := y0[i] - w * y1[i]
w := w * wn
return y
```
## 求答案
到這裡之後可以簡單地$O(n)$求出$C(x)$的點值表示法，只差把$C(x)$的係數算出來了，這個步驟稱為快速傅立葉逆變換IFFT(Inverse Fast Fourier Transform)
首先，因為計算DFT為一線性計算$(y_k=\sum\limits_{i=0}^{n-1}a_iw^{ki}_n)$，故可以將原式(係數)與結果(值)以一矩陣$V$表示。
$\begin{bmatrix}
y_0\\y_1\\y_2\\\vdots\\y_{n-1}
\end{bmatrix}=
\begin{bmatrix}
w^0_n&w^0_n&w^0_n&\cdots&w^0_n\\
w^0_n&w^1_n&w^2_n&\cdots&w^{n-1}_n\\
w^0_n&w^2_n&w^4_n&\cdots&w^{2(n-1)}_n\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
w^0_n&w^{n-1}_n&w^{2(n-1)}_n&\cdots&w^{(n-1)(n-1)}_n
\end{bmatrix}
\begin{bmatrix}
a_0\\a_1\\a_2\\\vdots\\a_{n-1}
\end{bmatrix}$
那要用$y$矩陣返回計算$a$矩陣只需找到$V$的反矩陣乘在左邊即可
(找到$V^{-1}$使$V^{-1}V=I_n$，$I_n$為$n$階單位矩陣)
因為$V$的特殊性，發現若$V(i,\ j)=w^{ij}_n$則$V^{-1}(i,\ j)=\displaystyle\frac 1nw^{-ij}_n$
:::info
:::spoiler 證明
要證明：
藉由矩陣的乘法規則，得知
$\begin{align}
\sum\limits_{k=0}^{n-1}V^{-1}(x,\ k)\cdot V(k,\ y)&=\sum\limits_{k=0}^{n-1}\displaystyle\frac 1nw^{-xk}_n\cdot w^{ky}_n\\
&=\displaystyle\frac1n\sum\limits_{k=0}^{n-1}w^{k(y-x)}_n\\
&=\begin{cases}
\displaystyle\frac1n\sum\limits_{k=0}^{n-1}w^0_n=\frac1n\cdot n=1&,\ \text{if}\ x=y\\\\
\displaystyle\frac1n\sum\limits_{k=0}^{n-1}w^{k(y-x)}_n=\frac1n(\sum\limits_{k=0}^{\frac n2-1}w^k_n+\sum\limits_{k=\frac n2}^{n-1}w^k_n)\\\\\hspace{2.95cm}=\displaystyle\frac1n(\sum\limits_{k=0}^{\frac n2-1}w^k_n+\sum\limits_{k=0}^{\frac n2-1}w^{k+\frac n2}_n)&,\ \text{if}\ x\not=y\\\\\hspace{2.95cm}=\displaystyle\frac1n(\sum\limits_{k=0}^{\frac n2-1}w^k_n-\sum\limits_{k=0}^{\frac n2-1}w^k_n)=0
\end{cases}\\
&=I(x,\ y)
\end{align}$
得證
:::
因此
$\begin{bmatrix}
a_0\\a_1\\a_2\\\vdots\\a_{n-1}
\end{bmatrix}
=\begin{bmatrix}
\frac1nw^{-0}_n&\frac1nw^{-0}_n&\frac1nw^{-0}_n&\cdots&\frac1nw^{-0}_n\\
\frac1nw^{-0}_n&\frac1nw^{-1}_n&\frac1nw^{-2}_n&\cdots&\frac1nw^{-(n-1)}_n\\
\frac1nw^{-0}_n&\frac1nw^{-2}_n&\frac1nw^{-4}_n&\cdots&\frac1nw^{-2(n-1)}_n\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
\frac1nw^{-0}_n&\frac1nw^{-(n-1)}_n&\frac1nw^{-2(n-1)}_n&\cdots&\frac1nw^{-(n-1)(n-1)}_n
\end{bmatrix}
\begin{bmatrix}
y_0\\y_1\\y_2\\\vdots\\y_{n-1}
\end{bmatrix}$
即
$a_k=\displaystyle\frac1n\sum\limits_{i=0}^{n-1}y_iw^{-ki}_n$
比對計算$y$的式子：
$y_k=\displaystyle\sum\limits_{i=0}^{n-1}a_iw^{ki}_n$
可以發現到，我們只需把原本FFT中的$\displaystyle\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}$改成$\displaystyle\cos\frac{2\pi}{n}-i\sin\frac{2\pi}{n}$
之後每項$\times \displaystyle\frac1n$就好了
因為差別不大，我們可以將FFT與IFFT寫在同一個function，多傳入一個變數$inv$，若是FFT就傳0，IFFT則傳入1，以下是新的虛擬碼：
```python=
def FFT( A[], inv ) :
n := A.length()
if n == 1:
return A
A0 := { A[0], A[2], ... , A[n - 2] }
A1 := { A[1], A[3], ... , A[n - 1] }
y0 := FFT( A0, inv )
y1 := FFT( A1, inv )
y := {}
ang := 2 * pi / n
w := 1
wn := cos(ang) + i * sin(ang) * (-1 if inv else 1)
for i in [0, n / 2 - 1] do:
y[i] := y0[i] + w * y1[i]
y[i + n / 2] := y0[i] - w * y1[i]
if inv:
y[i] := y[i] / 2
y[i + n / 2] := v[i + n / 2] / 2
w := w * wn
return y
```
跟上面一開始的FFT只有兩個不同：
第一個是$\mathrm{line\ 12}$的`i * sin(ang)`如果是IFFT就多乘一個$-1$，這個在上面解釋過了
第二個是$\mathrm{line\ 16\sim18}$，如果是IFFT則/2。這是為了達到把每項$\times\frac1n$，比較直觀的做法是在呼叫FFT function之後手動遍歷陣列去乘$\frac1n$。
這裡的做法則是是因為在遞迴中，陣列的每個位置會被碰到$\log_2 n$次，每次都/2的話共會除以$2^{\log_2 n}=n$，就剛好達到目的了
## 蝴蝶操作
在FFT中，將$y_0[i]+w\times y_1[i]$存入$y[i]$、$y_0[i]-w\times y_1[i]$存入$y[i+\frac n2]$這個動作稱為蝴蝶操作(butterfly operation)
他畫成圖長這樣：
## 時間複雜度
到這裡可以來分析一下為什麼最上面說時間複雜度是$O(n\log n)$了($\log$以$2$為底)
其實很簡單：FFT會遞迴$\log n$次，而每次中不論是分離$A_0$和$A_1$、進行操作都是$O(n)$，因此總時間複雜度$O(n\log n)$
## 迭代型優化
用遞迴寫FFT還是太慢了。觀察FFT遞迴形成的樹(以$n=8$舉例)：
```graphviz
digraph{
node[shape=ellipse]"(a0, a1, a2, a3, a4, a5, a6, a7)"->{"a0, a2, a4, a6", "a1, a3, a5, a7"}
"a0, a2, a4, a6"->{"a0, a4", "a2, a6"}
"a1, a3, a5, a7"->{"a1, a5", "a3, a7"}
"a0, a4"->{"a0", "a4"}
"a2, a6"->{"a2", "a6"}
"a1, a5"->{"a1", "a5"}
"a3, a7"->{"a3", "a7"}
}
```
發現如果先將原本$a_0\sim a_{n-1}$轉成樹根左到右出現的順序(例如$n=8$就是上面圖中的$a_0,\ a_4,\ a_2,\ \dots ,\ a_7$)，就可以藉由**從下往上，取出該層各元素進行蝴蝶操作後放回去**，從而用迭代實現FFT
這裡基本上只有一個特別的，就是要如何把$0\sim n-1$轉成樹根的次序呢
可以觀察轉換前後的數的二進位
$0(000)\rightarrow 0(000)$
$1(001)\rightarrow 4(100)$
$2(010)\rightarrow 2(010)$
$3(011)\rightarrow 6(110)$
$4(100)\rightarrow 1(001)$
$5(101)\rightarrow 5(101)$
$6(110)\rightarrow 3(011)$
$7(111)\rightarrow 7(111)$
會發現轉換一個數就是把他的二進位反轉。
其實也不難理解，把二進位從右邊開始為第一位的話，第一次分治就是把第一位為$0$的分到左邊，$1$分到右邊；第二次是用第二位分...但正常$0\sim n-1$會是以最左邊為第一位然後用一樣的方式分，所以最後的結果就會剛好轉置
所以，又可以寫出虛擬碼：
```python=
def rev( x, n ) :
k := 0
for t in [0, n - 1] do:
if x <= 0:
break
k := k * 2
if x % 2 == 1:
k := k + 1
# 上面這三行用位元運算就是：k := k << 1 | (x & 1)
x := x / 2
return k
def rev_array( A[] ) :
n := sizeof(A)
r := log_2(n)
p := 0
for i in [0, n - 1] do:
p := rev(i, r)
if p > i: # 確保不會重複交換
swap(A[i], A[p])
return
def FFT( A[], inv ) :
n := sizeof(A)
lgn := log_2(n)
rev_array(A)
for s in [1, lgn] do: # 第幾層
m := 2 ^ s
ang := 2 * pi / m * (-1 if inv else 1)
wn := cos(ang) + i * sin(ang)
for k in [0, n - 1] by m do:
w := 1
for j in [0, m / 2 - 1] do:
u := A[k + j]
t := w * A[k + j + m / 2]
A[k + j] := u + t
A[k + j + m / 2] := u - t
if inv:
for i in A do:
i := i / n
return
```
用迭代的時間複雜度仍為$O(n\log n)$，但常數小了不少，且空間複雜度也縮小許多(不須開額外陣列)
## 大數乘法
至於要怎麼用FFT實現大數乘法呢？
其實很簡單，例如$236\times 345$想成$(2x^2+3x+6)\times(3x^2+4x+5)$，答案就是處理進位後的各係數
## 題目
[CSES Signal Processing](https://cses.fi/problemset/task/2113/)
C++程式碼：
:::spoiler `Recursive version`
```cpp=
#include <bits/stdc++.h>
using namespace std;
#define vec vector
#define cp complex<double>
void FFT(vec<cp> &a, bool inv) {
int n = a.size();
if(n == 1) return;
vec<cp> a0(n >> 1), a1(n >> 1);
for(int i = 0, p = 0; i < n; i += 2, p++) {
a0[p] = a[i];
a1[p] = a[i + 1];
}
FFT(a0, inv), FFT(a1, inv);
double ang = 2 * M_PI / n;
cp w(1, 0), wn(cos(ang), sin(ang) * (inv? -1 : 1));
for(int i = 0; i < n / 2; i++, w *= wn) {
cp t = w * a1[i];
a[i] = a0[i] + t;
a[i + n / 2] = a0[i] - t;
if(inv) {
a[i] /= 2;
a[i + n / 2] /= 2;
}
}
return;
}
int main() {
int n, m;
cin >> n >> m;
int maxn = 1;
while(maxn < (n + m)) maxn <<= 1;
maxn <<= 1;
vec<cp> a(maxn), b(maxn);
for(int i = 0; i < n; i++) cin >> a[i];
for(int i = m - 1; ~i; i--) cin >> b[i];
FFT(a, 0), FFT(b, 0);
for(int i = 0; i < maxn; i++) a[i] *= b[i];
FFT(a, 1);
for(int i = 0; i < (m + n - 1); i++) cout << llround(a[i].real()) << ' ';
}
```
:::
結果：**#6 TLE**
:::spoiler `Iterative version`
```cpp=
#include <bits/stdc++.h>
using namespace std;
#define vec vector
#define cp complex<double>
int rev(int x, const int &n) {
int k, t;
for(k = 0, t = 0; x || t < n; t++, x >>= 1) {
k = (k << 1 | (x & 1));
}
return k;
}
void rev_array(vec<cp> &a) {
int n = a.size(), r = log2(n);
for(int i = 0; i < n; i++) {
int x = rev(i, r);
if(x > i) swap(a[i], a[x]);
}
}
void FFT(vec<cp> &a, bool inv) {
int n = a.size();
rev_array(a);
for(int m = 2; m <= n; m <<= 1) {
double ang = 2 * M_PI / m * (inv? -1 : 1);
cp wn(cos(ang), sin(ang));
for(int k = 0; k < n; k += m) {
cp w = 1;
int s = (m >> 1);
for(int j = 0; j < s; j++, w *= wn) {
cp t = w * a[k + j + s], u = a[k + j];
a[k + j] = u + t;
a[k + j + s] = u - t;
}
}
}
if(inv) {
for(auto &i : a) i /= n;
}
}
int main() {
int n, m;
cin >> n >> m;
int maxn = 1;
while(maxn < (n + m)) maxn <<= 1;
maxn <<= 1;
vec<cp> a(maxn), b(maxn);
for(int i = 0; i < n; i++) cin >> a[i];
for(int i = m - 1; ~i; i--) cin >> b[i];
FFT(a, 0), FFT(b, 0);
for(int i = 0; i < maxn; i++) a[i] *= b[i];
FFT(a, 1);
for(int i = 0; i < (m + n - 1); i++) cout << llround(a[i].real()) << ' ';
}
```
:::
結果：**AC，max time：0.64s**
也放個大數乘法的程式碼：
:::spoiler `Big number multiply`
```cpp=
#include <bits/stdc++.h>
using namespace std;
#define vec vector
#define cp complex<double>
int rev(int x, const int &n) {
int k, t;
for(k = 0, t = 0; x || t < n; t++, x >>= 1) {
k = (k << 1 | (x & 1));
}
return k;
}
void rev_array(vec<cp> &a) {
int n = a.size(), r = log2(n);
for(int i = 0; i < n; i++) {
int x = rev(i, r);
if(x > i) swap(a[i], a[x]);
}
}
void FFT(vec<cp> &a, bool inv) {
int n = a.size();
rev_array(a);
for(int m = 2; m <= n; m <<= 1) {
double ang = 2 * M_PI / m * (inv? -1 : 1);
cp wn(cos(ang), sin(ang));
for(int k = 0; k < n; k += m) {
cp w = 1;
int s = (m >> 1);
for(int j = 0; j < s; j++, w *= wn) {
cp t = w * a[k + j + s], u = a[k + j];
a[k + j] = u + t;
a[k + j + s] = u - t;
}
}
}
if(inv) {
for(auto &i : a) i /= n;
}
}
vec<int> mul(string &a, string &b) {
if(a == "0" || b == "0") {
vec<int> ans = {0};
return ans;
}
int maxn = max(a.size(), b.size()), n = 1;
while(n < maxn) n <<= 1;
n <<= 1;
vec<cp> A(n), B(n);
//降冪轉成升冪
for(int i = 0, p = a.size() - 1; i < a.size(); i++, p--) A[i] = a[p] - '0';
for(int i = 0, p = b.size() - 1; i < b.size(); i++, p--) B[i] = b[p] - '0';
FFT(A, 0);
FFT(B, 0);
for(int i = 0; i < n; i++) A[i] *= B[i];
FFT(A, 1);
vec<int> ans(n);
for(int i = 0; i < n; i++) ans[i] = round(A[i].real());
int carry = 0;
for(int i = 0; i < n; i++) {
ans[i] += carry;
carry = ans[i] / 10;
ans[i] %= 10;
}
int k = n - 1;
while(k >= 0 && !ans[k]) k--;
return ans.resize(k + 1), ans;
}
void print(vec<int> v) {
for(int i = v.size() - 1; ~i; i--) cout << v[i];
cout << '\n';
}
int main() {
string a, b;
cin >> a >> b;
print(mul(a, b));
}
```
:::
## 心得
這是我第一次寫這麼長篇而且挺難的東西，解鎖許多成就，例如寫虛擬碼。
其實我也是邊寫邊學的，想說這樣可以比較方便理解，而的確有效。
不過畢竟是第一次寫，加上是數學容易使描述不當，或我以為很明顯的東西對其他人卻很難理解，所以我只能盡量把自己在學時有遇到的瓶頸都寫出來，希望還算可以。
不過寫完才發現應該是$\omega$的我全部寫$w$
有錯誤還煩請告知
## 資料來源
* [CSDN 大數乘法（快速傅立葉變換）上](https://blog.csdn.net/u013351484/article/details/48739415)
* [CSDN 大數乘法（快速傅立葉變換）下](https://blog.csdn.net/u013351484/article/details/48809943)
* [知乎 快速傅立葉變換（FFT）求解多項式乘法](https://zhuanlan.zhihu.com/p/76622485)
* [維基百科 棣美弗公式](https://zh.wikipedia.org/zh-tw/%E6%A3%A3%E8%8E%AB%E5%BC%97%E5%85%AC%E5%BC%8F)
* [OI Wiki 快速傅立葉變換](https://oi-wiki.org/math/poly/fft/)
* [知乎 快速傅立葉變換的迭代法程式碼實現](https://zhuanlan.zhihu.com/p/422925896?utm_id=0)
*
---
> [name= 編者：柏霖 Boling]
Sign in with Wallet
Connect another wallet