30:T9f0,md`## Cooley-Tukey Algorithm 因為 $w_n$ 的性質，我們得以在 $O(n\log n)$ 的時間內求出 $n$ 次多項式 $A$ 在 $(w_n)^k$, $k = [0 .. (n-1)]$ 這 $n$ 個點的值。${sidenote`指數其實應該是 $-k$，不過順序顛倒過來會殊途同歸，所以這裡寫正的。`} 首先，我們可以在線性時間求多項式代入某給定點的值 (秦九韶方法)，將多項式整理成一乘一加的形式，迭代 $n$ 次即得到結果。 $$A(t) = a_0 + \underbrace{\Big( \cdots \big( (a_n t + a_{n-1}) t + a_{n-2} \big) t \cdots \Big)t}_{n \text{ levels}} $$ 乍看之下，計算 $n$ 個值需要 $O(n^2)$ 的時間，而快速傅立葉是一類能加速這個過程的演算法。這個又被稱為 **Base-2 FFT** 的演算法的基本思路是分治 (divide and conquer)，為了計算的需要，要假設 $n$ 是 $2$ 的整數次方${sidenote`這是因為每次對半切的關係。也可以改成每次切成大於兩塊，只是後續推導會有點不一樣，時間複雜度也會差一點。`}，如果不是的話可以補一些係數為 0 的項，不影響結果。不妨假設 $n = 2^\gamma$。首先，將 $A$ 係數依次方的奇偶性分成兩群，也就是 $$A(x) = A^{\text even} \left( x^2 \right) + x A^{\text odd}\left( x^2 \right)$$ 等號右邊可以看成是兩個關於 $x^2$ 的 $n/2$ 次多項式。這時的 $n$ 個點變為 $({w_n}^0)^2, {w_n}^2, \cdots, ({w_n}^{n-1})^2$，能表示成 $w_{n/2}$ 的冪次 (性質 3.)。次方縮小一半之外，厲害之處在於需要求 **$x^2$ 的數量也變成一半**，這是演算法能達到 $O(n \log n)$ 的關鍵！借用演算法筆記的一張圖： ![](https://web.ntnu.edu.tw/~algo/FastFourierTransform1.png) 再來是合併，注意到一個小細節是左半側和右半側的對應值差一個負號，所以前半和後半可以一起做 (性質 4.)，也就是： $$A \left( {w_n}^k \right) = \begin{cases} A^{\text even} \left( {w_n}^{2k} \right) + {w_n}^k A^{\text odd} \left( {w_n}^{2k} \right) & , k < n/2 \\ A^{\text even} \left( {w_n}^{2k'} \right) - {w_n}^{k'} A^{\text odd} \left( {w_n}^{2k'} \right) & , k = n/2 + k' \text{ for some } k' \ge 0 \end{cases} $$ 這麼做 $\gamma$ 輪之後只會剩下常數項，直接回傳就好。最後一層係數排列的順序，恰好是次方以二進位表示後，將 bits 顛倒過來寫所成的新數為次方排列的結果 (*bit reversal permutation*)。我們可以預先排列好，做完就會是正確的順序。 `31:T5f5,md`## 反向快速傅立葉變換用矩陣來表達上述過程，FFT 是在「分解」(factorize) 以下式子左側，使其易於計算： $$ \begin{bmatrix} 1 & 1 & 1 & \cdots & 1 \\ 1 & w_n & {w_n}^2 & \cdots & {w_n}^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & {w_n}^{n-1} & {w_n}^{2(n-1)} & \cdots & {w_n}^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} y_0 \\ y_1 \\ \vdots \\ y_{n-1} \end{bmatrix} $$ 把式子記成 $V\mathbf{a} = \mathbf{y}$，要從 $\mathbf{y}$ 解回 $\mathbf{a}$ 要做的事是左乘 $V^{-1}$。如果我們算出的 $V^{-1}$ 夠漂亮，那麼 FFT 的技巧就也能用在 inverse FFT。事實上，可以證明 $$V^{-1} = \overline{V}/n.$$ > **Claim**: $V \overline{V} = n I_n$. > **Proof**: > 考慮 $V \overline{V}$ 的第 $(i, j)$ 個元，即 $$\sum_{k=0}^{n-1} \left( {w_n}^{ik} \cdot {w_n}^{-kj} \right) = \sum_{k=0}^{n-1} {w_n}^{k(i-j)} = \begin{cases} n, &i = j;\\ 0, &\text{otherwise.}\end{cases}$$ > > 其中 $i - j = c \ne 0$ 的時候原式總是可以寫成 $$1 + {w_n}^c \cdots + {w_n}^{(n-1)c} = \frac{({w_n}^n)^c - 1}{w_n - 1} = 0 \;\; (\because w_n \ne 1).$$ 所以 inverse FFT 的時候，只需要對 $w_n$ 的次方部分取負號 (相當於取共軛)，再對解出來的向量 $\mathbf{y}$ 的各元同除以 $n$ 即可。${sidenote`訊號處理時比較注重 normalizing，而在正向反向都同乘或同除以 $\sqrt n$，但實務上常常考量計算簡潔性而不這麼做。`} `32:T78e,md`## 數論轉換在一個有限體下，我們也可以找到一些點來做快速傅立葉變換，這就是數論轉換。這裡為了簡潔起見，先只考慮階數 (order) 為 $P = c \cdot 2^k + 1$ 的有限體，$P$ 為**質數**。首先要找一個數 $r$ 滿足 $$\{r \bmod P, r^2 \bmod P, \cdots, r^{P-1} \bmod P\} = \{1, 2, \cdots, (P-1)\}. $$費馬小定理 (Fermat's Little Theorem) 已經保證在 $r$ 與 $P$ 互質的前提下，$r^{P-1} = 1$，可推論對所有 $m \in \{ 1\ .. (P-2) \}$，$r^m \ne 1 \mod P$，也就是說 $r$ 是 $\mathbb F_P$ 的**原根** (primitive root)。很不幸地，目前沒有高效率的演算法可以找 $r$${sidenote`高斯證明這樣的元素一定存在。另外，一共有 $\phi(\phi(P))$ 個這樣子的元素。 [這頁](https://cp-algorithms.com/algebra/primitive-root.html)有詳細的說明和演算法。`}，只能一個一個試；還好 $P$ 通常是給定的，所以 $r$ 允許預先算好。找到 $r$ 之後，所有 $\mathbb{F}_P$ 的非零元素都能表示成 $r$ 的冪次。再來是用 $r$ 表示 $\mathbb{F}_P$ 的 $2^k$ 階單位根，也就是某個數 $x$ 滿足 $x^{2^k} = 1$。顯然這個數是 $r^c$，因為 $(r^c)^{2^k} = r^{c \cdot 2^k} = r^{P-1} = 1 \mod P$。這樣子 $r^c$ 就扮演了前文中 $w_n$ 的角色： 1. 對於任意 $x \in \left[ 1, 2^k \right)$，$(r^c)^x \ne 1 \mod P$：因為 $1 \le cx < c \cdot 2^k = P-1$。 2. 對於相異 $\alpha, \beta \in \left[ 0, 2^k \right)$，$r^{c\alpha} \ne r^{c\beta}$：不妨假設 $\alpha < \beta$。假設兩者相等，那就有 $r^{c(\alpha + 2^k - \alpha)} = r^{c(\beta + 2^k - \alpha)}$，得 $r^{c(\beta - \alpha)} = 1 \bmod P$，但由前述，這是不會發生的。事實上，$r^c$ 模 $P$ 生成的乘法群同構於 $(\mathbb Z_{2^k})_+$。把反向變換的 "$/n$" 改成「乘上 $n$ 的乘法反元素」，然後就可以愉快地在有限體下做 FFT 了。 `33:T697,md`### 實例考慮在 $\mathbb{F}_{97}$ 下計算兩個多項式的乘積：$$ \begin{aligned} f(x) &= x^3 + 4x^2 + 5x + 7, \\ g(x) &= 2x^3 + 8x^2 + 16x + 31. \end{aligned}$$ 首先分解 $97 - 1 = 96 = 3 \cdot 2^5$。因為 $5^{96/2} = 96 \ne 1$ 和 $5^{96/3} = 35 \ne 1$，可以證明 $5$ 是一個原根。因為乘積會有 $7$ 項係數，我們要取 $2(3+1) = 8$ 次單位根，也就是 $c = 96/8 = 12$，所以取 $w = 5^c = 64$。預先算好所有 $w$ 的冪次，表記如下： ${(() => { const iter = [].map.bind([...new Array(8).keys()]) return smallTable([iter(x => tex`w^${x}`), iter(x => tex`${Math.pow(64, x) % 97}`)]) })()} #### NTT

NTT 計算過程

$f(x) = (4x^2+7) + x(x^2+5)$， 1. $4w^2+7 = 95, w^2+5 = 27$，故 $f(w) = 95 + 64 \cdot 27 = 77, f(w^5) = 95 - 64 \cdot 27 = 16,$ 2. $4w^4+7 = 3, w^4+5 = 4$，故 $f(w^2) = 3 + 22 \cdot 4 = 91, f(w^6) = 3 - 22 \cdot 4 = 12,$ 3. $4w^6+7 = 16, w^6+5 = 80$，故 $f(w^3) = 16 + 50 \cdot 80 = 39, f(w^7) = 16 - 50 \cdot 80 = 90,$ 4. $4w^8+7 = 11, w^8+5 = 6$，故 $f(w^4) = 11 + 96 \cdot 6 = 5, f(w^8) = 11 - 96 \cdot 6 = 17.$

所以 $f(x)$ 經過數論轉換後是 $(17, 77, 91, 39, 5, 16, 12, 90)$。對 $g(x)$ 進行同樣操作，得到 $(57, 70, 40, 7, 21, 53, 6, 91)$。 #### Convolution 將對應項相乘得到 $(96, 55, 51, 79, 8, 72, 72, 42)$。 #### Inverse NTT 將上面的 $w$ 換成反元素 $w^{-1} = 47$，用這個值對上一步結果進行 NTT 得到 $(87, 2, 43, 28, 76, 31, 16, 0)$，再乘上 $8^{-1} = 85$ 得到 $(23, 73, 66, 52, 58, 16, 2, 0)$，表示 $f(x)g(x) = 2x^6 + 16x^5 + 58x^4 + 52x^3 + 66x^2 + 73x + 23$。 `34:[] 2e:{"title":"數論轉換 (NTT)","cells":[{"id":0,"value":"md`# 數論轉換 (NTT)`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6,"value":"md`\n\n\n> Andy Pan ()\n> Feb 11, 2020\n\n\n\n>\n> *A reading report. Trying to write on Observable.*\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":314,"value":"toc('h2,h3')","pinned":false,"hidden":false,"mode":"ojs"},{"id":47,"value":"md`## 前情提要\n一般想要計算兩個 $n$ 次多項式 $A(x), B(x)$ 的乘法需要 $O(n^2)$ 的時間乘遍所有項，但是我們從高中數學又知道：「一個單變數的 $n$ 次多項式可以由它在 $n+1$ 個點上的取值來唯一決定。」利用這個觀察來計算多項式乘法，我們可以有個大膽的想法：\n\n1. 選取一些點 $(X_n) = (x_0, x_1, \\cdots, x_n)$.\n2. 分別計算 $A(x), B(x)$ 在 $X$ 上的值，用 $(n+1)$-tuple 來表示：$(\\alpha_0, \\alpha_1, \\cdots, \\alpha_n)$, $(\\beta_0, \\beta_1, \\cdots, \\beta_n)$.\n3. 將對應的值乘起來：$(\\gamma_0, \\gamma_1, \\cdots, \\gamma_n) = (\\alpha_i \\times \\beta_i)_{i=0 .. n}$.\n4. 將這個 $(n+1)$-tuple 轉換成對應的多項式。${sidenote`這樣其實不夠，因為 $C(x)$ 會是 $2n$ 次多項式，至少需要 $2n+1$ 個點才能唯一確定。各取 $2(n+1)$ 個點總行了吧。為了簡化起見，下文討論都敘述取 $n$ 個點的值，實務上可能會差一個常數倍。`}\n\n我們可以透過**快速傅立葉變換** (Fast Fourier Transform, FFT) 以及其反向變換 (Inverse FFT)，分別加速 2. 和 4. 兩步，讓整個演算法在 $O(n \\log n)$ 的時間下完成，但運算過程中會牽涉到複數和浮點數。如果定義域是有限體，可以使用衍生自此的**數論轉換** (Number Theoretic Transform, NTT)，同樣是一種快速計算摺積 (convolution) 的演算法，因為每一步操作都是精確的，許多操作甚至能查表完成，大大減少運算時間。\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":53,"value":"md`## 複習離散傅立葉變換\n\n定義 $$w_n = e^{2 \\pi i / n},$$\n\n數列 $\\{a_t\\}_n$ 的離散傅立葉變換 (Discrete Fourier Transform, DFT) 為數列 $\\{y_t\\}_n$，其中 $$y_t = \\sum_{k=0}^{n-1} a_t w_n^{-kt}.$$\n\n這裡想做的事情可以看作是 DFT 的一類特例：對於某個多項式，$a_t$ 是其 $x^t$ 項係數。\n\n### 幾個基本性質\n\n因為 \\${w_n}^0, {w_n}^1, \\cdots, {w_n}^{n-1}$ 是 $x^n = 1$ 的 $n$ 個根，在複數平面上依序為單位圓內接正 $n$ 邊形的頂點。由其代數和幾何關係不難得到以下觀察：\n\n1. 對於所有整數 $n, d, k \\ge 0$，\\${w_{dn}}^{dk} = {w_n}^k$.\n2. 對於所有偶數 $n > 0$，\\${w_n}^{n/2} = w_2 = -1$.\n3. 對於所有偶數 $n > 0$，\\$({w_n}^k)^2 = ({w_n}^{k + n/2})^2 = {w_{n/2}}^k$.\n4. 對於所有整數 $n, k \\ge 0$，\\${w_n}^{k + n/2} = -{w_n}^k$.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":102,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":467,"value":"md`\n> **Input**: An array A[i] of $n = 2^\\gamma$ elements.\n> **Output**: An array of the same size sampled on $\\omega^k$, $k = [0..n-1]$.\n> \\`\\`\\`plain\n> do BIT-REVERSE:\n> for i = 0 .. n-1:\n> if i < BR(i)：\n> swap(A[i], A[BR(i)]) // why this is correct?\n> \n> s = 2\n> while s < n: // each stage we merge a pair of (s/2) elements\n> ω = 1, ω_step = e^(2 pi / s)\n> for j = 0 .. (s/2 - 1):\n> for k = j .. (n/2-1) step m:\n> u = A[k] // A^even\n> t = A[k + s/2] * ω // A^odd\n> A[k] = u + t // }\n> A[k + s/2] = u - t // }- the so-called \"butterfly operation\"\n> ω *= ω_step\n> s *= 2\n> return A\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":157,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":233,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":497,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":13,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":17,"value":"md`# References\n\n* [Tutorial on FFT/NTT — The tough made simple. ( Part 1 )](https://codeforces.com/blog/entry/43499)\n * [*3. Divide & Conquer: FFT*: MIT OpenCourseWare](https://www.youtube.com/watch?v=iTMn0Kt18tg)\n* [Tutorial on FFT/NTT — The tough made simple. ( Part 2 )](https://codeforces.com/blog/entry/48798)\n* [演算法筆記：Waves](http://www.csie.ntnu.edu.tw/~u91029/Wave.html#4)\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"md`---`","pinned":false,"hidden":false,"mode":"ojs"},{"id":528,"value":"smallTable = xss => {\n const el = html``\n xss.forEach(xs => {\n const row = html``\n xs.forEach(x => row.appendChild(html``))\n el.appendChild(row)\n })\n return el\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":394,"value":"","pinned":true,"hidden":false,"mode":"html"},{"id":35,"value":"mdit.set({breaks: true})","pinned":true,"hidden":false,"mode":"ojs"},{"id":383,"value":"sidenote = Sidenote(md_)","pinned":true,"hidden":false,"mode":"ojs"},{"id":386,"value":"md = SidenoteWrapper(md_)","pinned":true,"hidden":false,"mode":"ojs"},{"id":311,"value":"import {toc} from '@bryangingechen/toc'","pinned":false,"hidden":false,"mode":"ojs"},{"id":364,"value":"import {Sidenote, SidenoteWrapper} from '@andy0130tw/sidenotes'","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"import {$r, md as md_, mdit} from '@andy0130tw/my-markdown-it'","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}

${x}