11.8聽課記錄 && FFT

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) \]

這里有一個消去引理，就是 \(w_n^{k+\frac n2}=-w_n^k\)，因為 \(w_n^{k+\frac n2}=e^{\frac{2\pi i(k+\frac n2)}{n}}=e^{\frac{2\pi ik}{n}+\pi i}\)，而\(e^{\pi i}=-1\)，於是就得證了。

我們以 \(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++;
for(int i=0;i<n;i++){
	rev[i]=(rev[i>>1]>>1)|((i&1)<<(bit-1));
	if(i<rev[i])swap(a[i],a[rev[i]]);//不加這條if會交換兩次（就是沒交換）
}

反正很神妙就對了。

---