捲積
別捲了
(離散)捲積是啥
卷積定義:
$${\displaystyle y[n]=f[n]*g[n]=\sum _{m=0}^{M-1}f[n-m]g[m]}$$
我超,多項式乘法
所以今天講多項式乘法
多項式乘法
直接做就是 \(O(n\times m)\)
好慢喔
點值表示法
我超,\(O(n)\)
多項式
常見的多項式表示法都是係數表示法
例如說 \(3x^2+2x+3\)
但其實有另一種表示法,叫做點值表示法
\(N+1\)個點唯一決定一個\(N\)次多項式
所以我可以用(-1,4), (0,3),(1,8)來表示上面那個方程式
簡單乘法
有甚麼好處?
我可以把\(f(x)\)與\(g(x)\)的點值相乘,再把新的座標插值成一個\(h(x)\),而\(h(x)=f(x)\times g(x)\)
但是
直接算點值是 \(O(n^2)\)
插值回去也是\(O(n^2)\)
但是我們可以找一些好處理的點
單位根
我超,複平面
單位根
單位根表示為 \({\omega}^n_k\)
\({\omega}^n_k = e^{i \times 2 \pi \times \frac{n}{k}}\)
其實就是在單位圓上,從\(1+0i\)開始逆時針繞\(\frac{n}{k}\)圈
complex
C++有個資料型態叫complex
可以存複數的
然後可以自己處理複數運算
但是有點慢
所以要速度的話可以用pair然後自己寫複數運算
喔對了後面提到的所有跟\(e^{i\pi}\)有關的東西
在程式裡面都要改成三角函數的形式
DFT -> FFT
我超,分治
DFT
所以我們現在想把把單位根帶進多項式裡求點值
好其實這就是 DFT
$${\displaystyle {\hat {x}}[k]=\sum _{n=0}^{N-1}e^{-i{\frac {2\pi }{N}}nk}x[n]\qquad k=0,1,\ldots ,N-1.}$$
稍微換個位置
$${\displaystyle {\hat {x}}[k]=\sum _{n=0}^{N-1}e^{{(-i{\frac {2\pi }{N}}k)}^n}x[n]\qquad k=0,1,\ldots ,N-1.}$$
\(e^{-i \frac{2\pi}{N}k}\)就是單位根
加個負號只是換個方向轉而已
還是 \(O(n^2)\)
怎麼加速
平方
\(f(x)=a_0+a_1 x + a_2 x^2 + a_3 x^3... a_{n-1} x^{n-1}\)
=\((a_0 + a_2 x^2 + a_4 x^4...) + (a_1 x+ a_3 x^3 + a_5 x^5...)\)
=\((a_0 + a_2 x^2 + a_4 x^4...) + x(a_1 + a_3 x^2 + a_5 x^4...)\)
定義\(f_{even}(x)=a_0+a_2 x + a_4 x^2...\)
\(f_{odd}(x)=a_1+a_3 x + a_5 x^2...\)
所以 \(f(x) = f_{even}(x^2)+x \times f_{odd}(x^2)\)
砍半
這樣就只要代單位根的平方進去
阿有差嗎
差了一半


原本所有單位根加起來轉一圈,現在轉動的量多一倍,變成轉兩圈,並且第一圈與第二圈的點完全重疊
所以最後的運算就只要算一半單位根的點值
\(sequence(e^{i2 \pi \frac{k}{10}},k,0,9)\)
\(sequence(e^{(i2 \pi \frac{k}{10})^2},k,0,9)\)
代回去
第一圈的單位根(\(\omega_n^k, k < \frac{n}{2}\))
\(f(\omega_n^k) = f_{even}(\omega^k_{n/2})+\omega^k_nf_{odd}(\omega^k_{n/2})\)
第二圈的單位根(\(\omega_n^{k+n/2}, k < \frac{n}{2}\))
\(f(\omega_n^{k+n/2}) = f_{even}(\omega^k_{n/2})-\omega^k_nf_{odd}(\omega^k_{n/2})\)
\(\omega_n^{k+n/2}=-\omega_n^{k}\)
因為少轉半圈
所以現在只要求\(f_{even}(\omega^k_{n/2})\)跟\(f_{odd}(\omega^k_{n/2})\)
哇喔,剩一半了
分治繼續砍下去,\(O(n log(n))\),複雜度算法跟 merge sort 一樣
IDFT -> FFT
我超,分治
正交性
$$\sum_{k=0}^{n-1}\omega^{m k}={\binom{n,\ \ \mathrm{if}\ n\ |\ m}{\mathrm{0, otherwise}}}$$
一堆均勻散佈在單位圓上的點,向量加起來全部抵銷
神秘推導
$$\sum_{k=0}^{n-1}y_{k}\,\omega^{-j k}=\sum_{k=0}^{n-1}\left(\sum_{\ell=0}^{n-1}c_{\ell}\,\omega^{\ell k}\right)\omega^{-j k}=\sum_{\ell=0}^{n-1}c_{\ell}\sum_{k=0}^{n-1}\omega^{(\ell-j)k}=n\cdot c_{j}$$
DFT,\(c_j\)為係數,\(y_k\)為點值
$$y_k=\sum_{j=0}^{n-1}c_j \omega^{jk}$$
$$\sum_{k=0}^{n-1}y_{k}\omega^{-jk}=n\cdot c_{j}$$
$$\frac{1}{n}\sum_{k=0}^{n-1}y_{k}\omega^{-jk}=c_{j}$$
神秘推導
$$c_{j}=\frac{1}{n}\sum_{k=0}^{n-1}y_{k}\omega^{-jk}$$
哇喔,跟DFT長得好像喔
所以其實跟DFT是用一樣的方法
$${\displaystyle x\left[n\right]={1 \over N}\sum _{k=0}^{N-1}e^{i{\frac {2\pi }{N}}nk}{\hat {x}}[k]\qquad n=0,1,\ldots ,N-1.}$$
遞迴
偷上上屆會爛掉的實作
using ld = long double;
using CD = complex<ld>;
CD OMEGA (int i, int n) {
ld angle = i;
angle = angle / n * (2 * pi);
return CD(cos(angle), sin(angle));
}
void FFT (vector <CD>& a, bool inv = 0) {
int n = a.size();
if (n == 1) return;
vector <CD> even(n / 2), odd(n / 2);
for (int i = 0; i < n / 2; i++) even[i] = a[i * 2], odd[i] = a[i * 2 + 1];
FFT(even, inv); FFT(odd, inv);
for (int i = 0; i < n / 2; i++) {
CD z = OMEGA(i * (inv?-1:1), n);
a[i] = even[i] + odd[i] * z;
a[i + n / 2] = even[i] - odd[i] * z;
if (inv) a[i] /= 2, a[i + n / 2] /= 2;
}
return;}題目
三次轉二次
我超,速度爆打NTT
複數乘法
DFT把係數轉成複數點值後,要先把兩個函式的點值互乘,再用IDFT
但是注意到\((a+bi)^2=a^2-b^2+2abi\)
如果我把\(f(x)\)的係數放在實數,\(g(x)\)放在虛數
然後做DFT後自己平方,\(2ab\)就會出現在虛數位
最後做一次IDFT再多除以二後,把虛數位抓出來也是答案
這樣就只要做一次DFT與IDFT
迭代式
我超,PUPA
FFT
By ck11300768鄭博軒
FFT
- 9