FSYo講數學+FFT，Orz

---

前置

傅立葉變換

（這裡傅立葉變換不理解不影響FFT的學習）

先看 3B1B的傅立葉變換

泰勒展開，是用一個多項式去擬合一個函式 \(x_0\) 處的值，在較小的範圍內能夠比較接近。所以需要做到每次求導都和原函式相同，於是有

\[g(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)(x-x_0)^2}{2}+\frac{f'''(x_0)(x-x_0)^3}{6}+\cdots+\frac{f^{(n)}(x_0)(x-x_0)^n}{n!}+\cdots \]

分母是後面的 \((x-x_0)^n\) 求導所產生的常數帶來的影響。這樣求每階導都與原函式相同。

我們泰勒展開 \(e^{ix},isinx,cosx\)，其中 \(i\) 是複數，分別得到：

\[\begin{align*} &e^{ix}=1+ix-\frac{x^2}{2!}-\frac{ix^3}{3!}+\frac{x^4}{4!}+\frac{ix^5}{5!}\cdots\\ &i\sin x=\frac{ix}{1!}-\frac{ix^3}{3!}+\frac{ix^5}{5!}-\frac{x^7}{7!}\cdots\\ &\cos x=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}\cdots \end{align*} \]

可以發現 \(e^{ix}=\cos x+i\sin x\)

，這就是尤拉公式。

將 \(\pi\) 帶入 \(x\) 就可以得到 \(e^{i\pi}+1=0\)

\(e^{ix}=\cos x+i\sin x\) 在複平面上就可以形象地理解成一個圓。隨著 \(x\) 增大，\(e^{ix}\) 就在繞著圓旋轉。

於是我們就能更清晰地理解傅立葉變換公式

\[\int_{t1}^{t2}g(t)e^{-2\pi ift}\mathrm dt \]

如果寫成離散形式並用 \(O(n^2)\) 計算就是離散傅立葉變換（DFT），使用神奇的單位根優化做到 \(O(n\log n)\) 計算就是我們的快速傅立葉變換（FFT）。

---

複數運算

一個複數 \(z=a+bi\)

，其中 \(a\) 是他的實部，\(b\) 是他的虛部，\(i\) 稱為虛部單位。

一個複數 \(z=a+bi\) 的共軛複數是 \(\bar{z}=a-bi\) 。

複數有四則運算

\[\begin{align*} &(a\pm bi)\pm (c\pm di)=(a\pm c)\pm (b\pm d)i\\ &(a+b_i)(c+d_i)=(ac-bd)+(ad+bc)i\\ &\frac{a+b_i}{c+d_i}=\frac{(a+b_i)(c-d_i)}{(c+d_i)(c-d_i)}=\frac{(ac+bd)+(bc-ad)i}{c^2+d^2} \end{align*} \]

---

多項式卷積

\[h(k)=\sum_{i+j=k}f(i)g(j) \]

直接做顯然是 \(O(n^2)\).

---

拉格朗日插值

（這裡拉格朗日插值不理解不影響FFT的學習）

根據點值插出多項式的演算法.

構造出 \(n\) 個函式，使得第 \(i\) 個函式只在 \(x=x_i\) 時有值為 \(y_i\)，再將這 \(n\) 個函式求和。

由於次數為 \(n-1\) 且過這 \(n\) 個點的多項式函式只有一個，所以這個求和的函式就是那個正確的。

時間複雜度 \(O(n^2)\)，使用高斯消元是 \(O(n^3)\).

\[f(x)=\sum_{i=1}^n y_i\prod_{j\ne i}^n \frac{x-x_j}{x_i-x_j} \]

---

單位復根

\(n\) 次單位復根是滿足 \(w^n=1\) 的複數 \(w\)，顯然這樣的 \(w\) 有 \(n\) 個，均勻分佈在複平面的單位圓上。

我們在尤拉公式 \(e^{ix}=\cos x+i\sin x\) 中，把 \(2\pi\) 帶入 \(x\)，可以得到 \(e^{2i\pi}=1\)，那麼 \(w=e^{\frac{2\pi i}{n}}\)，將這個 \(w\) 記做主次單位根 \(w_n^1\)，那麼其他的單位根就是 \(w_n^k=e^{\frac{2\pi ik}{n}}(0\leq k<n)\)，顯然 \(w_n^0=1\)。

---

FFT

快速傅立葉變換（FFT）用於 \(O(n\log n)\) 時間內求兩個多項式之積。

我們有兩個多項式，（項數不同用 \(0\) 補足）

\[A(x)=a_0+a_1x+a_2x^2+a_3x^3+\cdots a_nx^n\\ B(x)=b_0+b_1x+b_2x^2+b_3x^3+\cdots b_nx^n \]

直接做是無法 \(O(n\log n)\) 做的，於是我們找到兩個多項式的點值表示式

\[A(x)=(x_0,y_0),(x_1,y_1),\cdots,(x_n,y_n)\\ B(x)=(x_0,y'_0),(x_1,y'_1),\cdots,(x_n,y'_n) \]

那麼他們相乘得到的多項式是

\[C(x)=(x_0,y_0y'_0),(x_1,y_1y'_1),\cdots,(x_n，y_ny'_n) \]

也就是能夠在 \(O(n)\) 時間內做乘法。

所以我們只需要能夠將係數表示式快速轉成點值表示式，以及將點值表示式快速轉成係數表示式的辦法。

---

\(O(n\log n)\) 求點值

我們需要將大於 \(2n\) 個點帶入多項式計算點值，這裡我們將單位根帶入多項式求點值，並強制選 \(2^t\) 個單位根，也就是保證 \(n\) 是 \(2\) 的冪，那麼這些單位根是 \(w_{2^t}^k=e^{\frac{2\pi ik}{2^t}}(0\leq k<2^t)\)，看起來很麻煩，所以還是把 \(2^t\) 寫成 \(n\)。

\[w_n^k=e^{\frac{2\pi ik}{n}}(0\leq k<n). \]

顯然有 \(w_n^j\times w_n^k=w_n^{(j+k)\%n}\)。

那麼多項式 \(A(x)\) 的點值就是：

\[\begin{align*} &A(w_n^0)=a_0+a_1+a_2+\cdots+a_{n-1}\\ &A(w_n^1)=a_0+a_1w_n^1+a_2w_n^2+\cdots+a_{n-1}w_n^{n-1}\\ &A(w_n^2)=a_0+a_1w_n^2+a_2w_n^4+\cdots+a_{n-1}w_n^{2(n-1)}\\ &A(w_n^3)=a_0+a_1w_n^3+a_2w_n^6+\cdots+a_{n-1}w_n^{3(n-1)}\\ &\vdots \end{align*} \]

考慮將每個點值按奇偶分開，比如

\[A(w_n^i)=(a_0+a_2w_n^{2i}+a_4w_n^{4i}+\cdots+a_{n-2}w_n^{(n-2)i})+(a_1w_n^i+a_3w_n^{3i}+\cdots+a_{n-1}w_{n-1}^{(n-1)i}) \]

發現右邊可以提一個 \(w_n^i\)，於是

\[A(w_n^i)=(a_0+a_2w_n^{2i}+a_4w_n^{4i}+\cdots+a_{n-2}w_n^{(n-2)i})+w_n^i(a_1+a_3w_n^{2i}+\cdots+a_{n-1}w_{n-1}^{(n-2)i}) \]

然後可以把左邊視為係數為 \(a_0,a_2,a_4,\cdots,a_{n-2}\)，帶入 \(x\) 的為 \(w_n^{2i}\) 的一個項數是 \(\frac n2\) 的多項式，記為 \(FL(w_n^{2i})\)，把右邊視為係數為 \(a_1,a_3,a_5,\cdots,a_{n-1}\)，帶入 \(x\) 的為 \(w_n^{2i}\) 的一個項數是 \(\frac n2\) 的多項式，記為 \(FR(w_n^{2i})\)，然後 \(w_n^{2i}=w_{n/2}^{i}\) ，所以相當於帶入 \(w_{n/2}^{i}\)。（這就是強制 \(n=2^t\) 的原因）

所以我們要算 \(A(w_n^i)\) 的點值，只需要算 \(FL(w_{n/2}^i)\) 和 \(FR(w_{n/2}^i)\) 的點值，然後

\[A(w_n^i)=FL(w_{n/2}^i)+w_n^iFR(w_{n/2}^i) \]

我們以 \(n=8\) 為例把每個點值的算式寫出來

\[\begin{align*} &A(w_n^0)=FL(w_{n/2}^0)+w_n^0FR(w_{n/2}^0)\\ &A(w_n^1)=FL(w_{n/2}^1)+w_n^1FR(w_{n/2}^1)\\ &A(w_n^2)=FL(w_{n/2}^2)+w_n^2FR(w_{n/2}^2)\\ &A(w_n^3)=FL(w_{n/2}^3)+w_n^3FR(w_{n/2}^3)\\ &A(w_n^4)=FL(w_{n/2}^0)+w_n^0FR(w_{n/2}^0)\\ &A(w_n^5)=FL(w_{n/2}^1)+w_n^1FR(w_{n/2}^1)\\ &A(w_n^6)=FL(w_{n/2}^2)+w_n^2FR(w_{n/2}^2)\\ &A(w_n^7)=FL(w_{n/2}^3)+w_n^3FR(w_{n/2}^3)\\ \end{align*} \]

注意到並不需要每個單位根的點值都這樣算一遍，因為前面一半的 \(FL，FR\) 和後面一半是一樣的，所以我們每次的計算規模是減半的，把所有單位根的點值一起計算，每次計算時 \(O(n)\)，所以總的時間複雜度是 \(T(n)=2T(\dfrac n2)+O(n)=O(n\log n)\)，我們就成功地在 \(O(n\log n)\) 的時間內把係數表示式轉成點值表示式了。以上稱為快速傅立葉變換。

---

\(O(n\log n)\) 求係數

·設原本的多項式是

\[A(w_n^x)=\sum_{i=0}^{n-1}a_iw_n^{xi} \]

我們對他進行傅立葉變換，得到 \(n\) 個點值，設 \(A(w_n^i)\) 為 \(b_i\)，我們把這些 \(b_i\) 當做係數再做一次 FFT，得到

\[C(w_n^x)=\sum_{i=0}^{n-1}b_iw_n^{xi} \]

把 \(b_i\) 換成多項式的式子

\[C(w_n^x)=\sum_{i=0}^{n-1}w_n^{xi}\sum_{j=0}^{n-1}a_jw_n^{ij} \]

套路交換求和號

\[C(w_n^x)=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}w_n^{i(j+x)} \]

利用單位根的性質，發現在複平面的單位元上，只有 \(j+x=n\) 時，後面的 \(\sum\) 才有值為 \(n\)，否則所有單位根之和總是 \(0\)，（特殊情況是 \(x=0\)，此時只有 \(j=0\) 才有值為 \(n\)），所以我們發現

\[\begin{align*} &C(w_n^0)=na_0\\ &C(w_n^1)=na_{n-1}\\ &C(w_n^2)=na_{n-2}\\ &\vdots\\ &C(w_n^{n-1})=na_1 \end{align*} \]

所以我們對用快速傅立葉變換所求出的點值當成係數再做一遍快速傅立葉變換，再得到的點值就可以表示出原本的係數了。以上稱為快速傅立葉逆變換。

---

優化

如果遞迴地FFT，常數會很大，我們把 FFT 的遞迴寫出來

發現每個位置在 FFT 最後一層的位置是他二進位制反轉後的位置。

我們可以 \(O(n)\) 預處理出每個位置最終的位置 rev[i] 然後從下往上合併上去，會比遞迴快很多。

預處理

int bit=0;
while ((1<<bit)<n)bit++;
fo(i,0,n-1){
   rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
   if (i<rev[i])swap(a[i],a[rev[i]]);//不加這條if會交換兩次（就是沒交換）
}

反正很神妙就對了。

還有不會的蝴蝶變換

---