捲積

別捲了

(離散)捲積是啥

卷積定義:

$${\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;}

題目

luoguP3803(裸題)

luoguP1919(多個進位)

CSES2111(想想)

TIOJ1064(神秘優化題)

三次轉二次

我超,速度爆打NTT

複數乘法

DFT把係數轉成複數點值後,要先把兩個函式的點值互乘,再用IDFT

但是注意到\((a+bi)^2=a^2-b^2+2abi\)

如果我把\(f(x)\)的係數放在實數,\(g(x)\)放在虛數

然後做DFT後自己平方,\(2ab\)就會出現在虛數位

最後做一次IDFT再多除以二後,把虛數位抓出來也是答案

這樣就只要做一次DFT與IDFT

迭代式

我超,PUPA

FFT

By ck11300768鄭博軒