FFT & NTT

22717 高翊恩

關於講師

OJ id: pring、pringDeSu
慢慢改掉喜歡壓一行的習慣了
左閉右開

先備知識

a.k.a. 數學

向量

沒事今天沒有幾何

\(\bold{v}^T=\left[\begin{array}{cccc}2&3&-5&7\end{array}\right]\)

\(\bold{v}=\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]\)

向量的線性組合

把對應位置相乘後加總

\(\bold{v}^T=\left[\begin{array}{cccc}2&3&-5&7\end{array}\right]\)

\(\bold{u}^T=\left[\begin{array}{cccc}4&2&0&-1\end{array}\right]\)

\(\bold{v}^T\cdot\bold{u}^T=2\cdot4+3\cdot2+(-5)\cdot0+7\cdot(-1)=7\)

矩陣

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\end{array}\right]\)

矩陣乘法

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\\3&-2&0&4\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}9\\-24\\-15\\28\end{array}\right]\)

矩陣乘法

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\\3&-2&0&4\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}9\\-24\\-15\\28\end{array}\right]\)

矩陣乘法

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\\3&-2&0&4\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}9\\-24\\-15\\28\end{array}\right]\)

矩陣乘法

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\\3&-2&0&4\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}9\\-24\\-15\\28\end{array}\right]\)

矩陣乘法

\(\left[\begin{array}{cccc}2&4&0&-1\\3&5&2&-5\\-4&6&5&0\\3&-2&0&4\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}9\\-24\\-15\\28\end{array}\right]\)

多項式

\(A(x)=2+3x-5x^2+7x^3\)

\(\bold{A}=\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]\)

帶數字進去

\(A(x)=2+3x-5x^2+7x^3\)

\(A(-3)=2+3\cdot(-3)-5\cdot(-3)^2+7\cdot(-3)^3=-5\)

\(\left[\begin{array}{cccc}(-3)^0&(-3)^1&(-3)^2&(-3)^3\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=-5\)

長得很像做內積

帶一堆數字進去

\(A(x)=2+3x-5x^2+7x^3\)

\(\left[\begin{array}{cccc}(-3)^0&(-3)^1&(-3)^2&(-3)^3\\2^0&2^1&2^2&2^3\\0^0&0^1&0^2&0^3\end{array}\right]\left[\begin{array}{c}2\\3\\-5\\7\end{array}\right]=\left[\begin{array}{c}-5\\25\\2\end{array}\right]\)

複數（平面）

\(\operatorname{cis}(\theta)=e^{i\theta}=\cos\theta+i\sin\theta\)

事先說好

\(1\)的\(n\)次方根中，由\(x\)軸正向開始逆時針掃，掃到的第\(k\)個解，稱為\(\omega_n^k\)

\(\omega_n^k=\operatorname{cis}(2\pi\cdot\frac{k}{n})\)

多項式的特性

先從問題開始

給你兩個\(n\)次多項式，請輸出它們相乘的結果

複雜度限制\(\Omicron(n^2)\)

先從問題開始

給你兩個\(n\)次多項式，請輸出它們相乘的結果

複雜度限制\(\Omicron(n\log n)\)

換個問題

有兩個未知的多項式\(A(x)\)和\(B(x)\)和一個長度為\(n\)的已知陣列\(x\)，已知

\(\left\{\begin{array}{l}A(x_i)=\alpha_i\\B(x_i)=\beta_i\end{array}\right.,1\leq i\leq n\)

令

\(c_i=(A*B)(x_i)\ ,1\leq i\leq n\)

求\(c\)

\(c_i=\alpha_i\cdot\beta_i\)

hehe

回到這個問題

給你兩個\(n\)次多項式，請輸出它們相乘的結果

複雜度限制\(\Omicron(n\log n)\)

跳躍性思考

挑一堆數字\(x_i\)給\(A\)和\(B\)，就會有\(\alpha_i\)和\(\beta_i\)

令\(c_i=\alpha_i\cdot\beta_i\)，就可以得到\((A*B)(x_i)\)

這相當於是\((A*B)(x)\)函數圖形上的一堆點，想辦法做一些奇怪的數學操作就可以把它變回來了

\(A(x)=x-1\)

\(B(x)=-x+2\)

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((0,2)\)

\((1,1)\)

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((0,2)\)

\((1,1)\)

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((0,2)\)

\((1,1)\)

這是\(A*B\)嗎？

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((2,1)\)

\((0,2)\)

\((1,1)\)

\((2,0)\)

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((2,1)\)

\((0,2)\)

\((1,1)\)

\((2,0)\)

\(A(x)=x-1\)

\(B(x)=-x+2\)

\((0,-1)\)

\((1,0)\)

\((2,1)\)

\((0,2)\)

\((1,1)\)

\((2,0)\)

一些問題

兩個次數為\((n-1)\)的多項式相乘，至少要帶多少個數字進去？

至少要帶幾個點才能決定唯一一個次數為\(n\)的多項式？

大統整

\(A(x),B(x)\)

\(\alpha,\beta\)

\(c\)

\((A*B)(x)\)

帶值進去

兩兩相乘

想辦法弄回來

實作這個東西

1. 帶值進去

令\(\bold{M}=\left[\begin{array}{ccccc}x_0^0&x_0^1&x_0^2&\cdots&x_0^{n-1}\\x_1^0&x_1^1&x_1^2&\cdots&x_1^{n-1}\\x_2^0&x_2^1&x_2^2&\cdots&x_2^{n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\x_{n-1}^0&x_{n-1}^1&x_{n-1}^2&\cdots&x_{n-1}^{n-1}\end{array}\right]\)

則\(\left\{\begin{array}{l}\bold{MA}=\bold{\alpha}\\\bold{MB}=\bold{\beta}\end{array}\right.\)

2. 兩兩相乘

就兩兩相乘阿

3. 想辦法變回來

我們可以找出\(\bold{M}^{-1}\)

則\(\bold{M}^{-1}\bold{c}=\bold{A}*\bold{B}\)

大統整 · 改

\(A(x),B(x)\)

\(\alpha,\beta\)

\(c\)

\((A*B)(x)\)

乘一個矩陣

兩兩相乘

乘一個反矩陣

複雜度分析

\(A(x),B(x)\)

\(\alpha,\beta\)

\(c\)

\((A*B)(x)\)

\(\Omicron(?)\)

\(\Omicron(n)\)

\(\Omicron(?)\)

如果我們可以讓「乘一個矩陣」在\(\Omicron(n\log n)\)的時間做完，

就有辦法在\(\Omicron(n\log n)\)的時間內將兩個多項式相乘！

傅立葉變換

DFT

要帶一堆數字

理論上我們要帶\(2n\)個數字進去

多項式帶值是\(\Omicron(n)\)

有沒有更快的做法？

於是我們就想到了 Divide and Conquer

要帶哪些數字

帶\(\omega_N^0,\omega_N^1,\omega_N^2,\dots,\omega_N^{N-1}\)這\(N\)個數字！

※溫馨提醒：\(N\)要是一個比\(2n\)大的「\(2\)的冪次」

所以說會有這個東西

\left[\begin{array}{ccccc}\omega_N^0&\omega_N^0&\omega_N^0&\cdots&\omega_N^0\\\omega_N^0&\omega_N^1&\omega_N^2&\cdots&\omega_N^{N-1}\\\omega_N^0&\omega_N^2&\omega_N^4&\cdots&\omega_N^{2(N-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\\omega_N^0&\omega_N^{N-1}&\omega_N^{2(N-1)}&\cdots&\omega_N^{(N-1)(N-1)}\end{array}\right] \left[\begin{array}{c}a_0\\a_1\\a_2\\\vdots\\a_{N-1}\end{array}\right]= \left[\begin{array}{c}y_0\\y_1\\y_2\\\vdots\\y_{N-1}\end{array}\right]

所以說會有這個東西

\(\bold{\Omega}_N\bold{A}=\operatorname{DFT}(\bold{A})=\bold{Y}\)

我們的目標是找出\(\bold{Y}\)

其實\(\bold{\Omega}_N\bold{A}\)可以\(\Omicron(N\log N)\)算完

成功的話就是 Fast Fourier Transform

當我們要算\(\bold{\Omega}_N\bold{A}\dots\)

若\(N=1\)，則\(\bold{\Omega}_N\bold{A}=\bold{A}\)

當我們要算\(\bold{\Omega}_N\bold{A}\dots\)

否則：

把\(\bold{A}\)拆成兩部分\(\bold{A}_0,\bold{A}_1\)
對它們兩個做\(\operatorname{DFT}\)，得到\(\bold{Y}_0,\bold{Y}_1\)
想辦法用\(\bold{Y}_0,\bold{Y}_1\)湊出\(\bold{Y}\)

\bold{A}=\left[\begin{array}{c}a_0\\a_1\\a_2\\a_3\\a_4\\a_5\\\vdots\\a_{N-2}\\a_{N-1}\end{array}\right]\Rightarrow\bold{A}_0=\left[\begin{array}{c}a_0\\a_2\\a_4\\\vdots\\a_{N-2}\end{array}\right],\bold{A}_1=\left[\begin{array}{c}a_1\\a_3\\a_5\\\vdots\\a_{N-1}\end{array}\right]

\operatorname{DFT}(\bold{A}_0)=\bold{Y}_0=\left[\begin{array}{c}y_{0,0}\\y_{0,1}\\y_{0,2}\\\vdots\\y_{0,n-1}\end{array}\right]; \operatorname{DFT}(\bold{A}_1)=\bold{Y}_1=\left[\begin{array}{c}y_{1,0}\\y_{1,1}\\y_{1,2}\\\vdots\\y_{1,n-1}\end{array}\right]

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

分析時間

y_{0,k}=\bold{A}_0(\omega_n^k)

y_{1,k}=\bold{A}_1(\omega_n^k)

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(x)\\ =&\bold{A}_0(x^2)+x\bold{A}_1(x^2)\\ \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(x)\\ =&\bold{A}_0(x^2)+x\bold{A}_1(x^2)\\ \end{array}

x=\omega_N^k

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&\bold{A}_0(\omega_N^{2k})+\omega_N^k\bold{A}_1(\omega_N^{2k})\\ \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&\bold{A}_0(\omega_N^{2k})+\omega_N^k\bold{A}_1(\omega_N^{2k})\\ \end{array}

\omega_N^{2k}=\omega_{2n}^{2k}=\omega_n^k

（若\(0\leq k<n\)）

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&\bold{A}_0(\omega_n^k)+\omega_N^k\bold{A}_1(\omega_n^k)\\ \end{array}

（若\(0\leq k<n\)）

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

到這裡已經好一半了

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2(k+n)})+\omega_N^{k+n}\bold{A}_1(\omega_N^{2(k+n)}) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2(k+n)})+\omega_N^{k+n}\bold{A}_1(\omega_N^{2(k+n)}) \end{array}

\omega_N^{2(k+n)}=\omega_N^{2k+N}=\omega_N^{2k}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2k})+\omega_N^{k+n}\bold{A}_1(\omega_N^{2k}) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2k})+(\omega_N^k\cdot\omega_N^n)\bold{A}_1(\omega_N^{2k}) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2k})+(\omega_N^k\cdot(-1))\bold{A}_1(\omega_N^{2k}) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2k})-\omega_N^k\bold{A}_1(\omega_N^{2k}) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_n^k)-\omega_N^k\bold{A}_1(\omega_n^k) \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&y_{0,k}-\omega_N^k y_{1,k} \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&y_{0,k}-\omega_N^k y_{1,k} \end{array}

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&y_{0,k}-\omega_N^k y_{1,k} \end{array}

\(y_{0,k},y_{1,k}皆已知\)

萬歲！

想通的話可以看看Code

typedef complex<double> cd;

cd cis(double theta) {
    return complex<double>(cos(theta), sin(theta));
}

cd omega(int n, int k) {
    return cis(acos(-1) * 2 / n * k);
}

vector<cd> FFT(vector<cd> &A) {
    int N = A.size();
    // N is the power of 2
    if (N == 1) return vector<cd>(1, A[0]);

    int n = N / 2;
    vector<cd> A0(n), A1(n);
    for (int i = 0; i < n; i++) {
        A0[i] = A[i * 2];
        A1[i] = A[i * 2 + 1];
    }

    vector<cd> Y0 = FFT(A0), Y1 = FFT(A1);
    vector<cd> Y(N);
    for (int k = 0; k < n; k++) {
        Y[k] = Y0[k] + omega(N, k) * Y1[k];
        Y[k + n] = Y0[k] - omega(N, k) * Y1[k];
    }
    return Y;
}

不回傳陣列的Code

typedef complex<double> cd;

cd cis(double theta) {
    return complex<double>(cos(theta), sin(theta));
}

cd omega(int n, int k) {
    return cis(acos(-1) * 2 / n * k);
}

void FFT(vector<cd> &A) {
    int N = A.size();
    // N is the power of 2
    if (N == 1) return;

    int n = N / 2;
    vector<cd> A0(n), A1(n);
    for (int i = 0; i < n; i++) {
        A0[i] = A[i * 2];
        A1[i] = A[i * 2 + 1];
    }

    FFT(A0), FFT(A1);
    for (int k = 0; k < n; k++) {
        A[k] = A0[k] + omega(N, k) * A1[k];
        A[k + n] = A0[k] - omega(N, k) * A1[k];
    }
}

FFT with No Recursion

\(a_0\)

\(a_1\)

\(a_2\)

\(a_3\)

\(a_4\)

\(a_5\)

\(a_6\)

\(a_7\)

]

\(a_0\)

\(a_1\)

\(a_2\)

\(a_3\)

\(a_4\)

\(a_5\)

\(a_6\)

\(a_7\)

[

]

[

]

\(a_0\)

\(a_4\)

\(a_2\)

\(a_6\)

\(a_1\)

\(a_5\)

\(a_3\)

\(a_7\)

[

]

\(a_0\)

\(a_4\)

\(a_2\)

\(a_6\)

\(a_1\)

\(a_5\)

\(a_3\)

\(a_7\)

[

]

[

]

\(y_0\)

\(y_1\)

\(y_2\)

\(y_3\)

\(y_4\)

\(y_5\)

\(y_6\)

\(y_7\)

]

[

\(a_0\)

\(a_1\)

\(a_2\)

\(a_3\)

\(a_4\)

\(a_5\)

\(a_6\)

\(a_7\)

]

[

\(a_0\)

\(a_4\)

\(a_2\)

\(a_6\)

\(a_1\)

\(a_5\)

\(a_3\)

\(a_7\)

[

]

[

]

\(y_0\)

\(y_1\)

\(y_2\)

\(y_3\)

\(y_4\)

\(y_5\)

\(y_6\)

\(y_7\)

]

[

\(a_0\)

\(a_1\)

\(a_2\)

\(a_3\)

\(a_4\)

\(a_5\)

\(a_6\)

\(a_7\)

]

[

\(a_0\)

\(a_4\)

\(a_2\)

\(a_6\)

\(a_1\)

\(a_5\)

\(a_3\)

\(a_7\)

[

]

[

]

\(y_0\)

\(y_1\)

\(y_2\)

\(y_3\)

\(y_4\)

\(y_5\)

\(y_6\)

\(y_7\)

]

[

奇怪排序 \(\longrightarrow\) 一段一段合併

範例

\([\)

\(4\)

\(0\)

\(2\)

\(-5\)

\(-2\)

\(6\)

\(3\)

\(-1\)

\(]\)

範例

\([\)

\(4\)

\(0\)

\(2\)

\(-5\)

\(-2\)

\(6\)

\(3\)

\(-1\)

\(]\)

\(0\)

\(1\)

\(2\)

\(3\)

\(4\)

\(5\)

\(6\)

\(7\)

範例

\([\)

\(4\)

\(0\)

\(2\)

\(-5\)

\(-2\)

\(6\)

\(3\)

\(-1\)

\(]\)

\(0\)

\(1\)

\(2\)

\(3\)

\(4\)

\(5\)

\(6\)

\(7\)

範例

\([\)

\(4\)

\(0\)

\(2\)

\(-5\)

\(-2\)

\(6\)

\(3\)

\(-1\)

\(]\)

範例

\([\)

\(4\)

\(0\)

\(2\)

\(-5\)

\(-2\)

\(6\)

\(3\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(0\)

\(2\)

\(-5\)

\(6\)

\(3\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(0\)

\(2\)

\(-5\)

\(6\)

\(3\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(0\)

\(5\)

\(-5\)

\(6\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(0\)

\(5\)

\(-5\)

\(6\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(6\)

\(5\)

\(-5\)

\(6\)

\(-6\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(6\)

\(5\)

\(-5\)

\(6\)

\(-6\)

\(-1\)

\(]\)

範例

\([\)

\(2\)

\(6\)

\(5\)

\(-6\)

\(6\)

\(-6\)

\(-1\)

\(-4\)

\(]\)

範例

\([\)

\(2\)

\(6\)

\(5\)

\(-6\)

\(6\)

\(-6\)

\(-1\)

\(-4\)

\(]\)

範例

\([\)

\(2\)

\(6\)

\(5\)

\(-6\)

\(6\)

\(-6\)

\(-1\)

\(-4\)

\(]\)

範例

\([\)

\(6\)

\(-6\)

\(-4\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

範例

\([\)

\(6\)

\(-6\)

\(-4\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

範例

\([\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

\(0\)

\(-6-4i\)

\(12\)

\(-6+4i\)

範例

\([\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

\(0\)

\(-6-4i\)

\(12\)

\(-6+4i\)

範例

\([\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

\(0\)

\(-6-4i\)

\(12\)

\(-6+4i\)

範例

\([\)

\(]\)

\(6-i\)

\(6+i\)

\(-3\)

\(7\)

\(0\)

\(-6-4i\)

\(12\)

\(-6+4i\)

自己做

Code結構

奇怪排序
做\(\operatorname{lg}N\)次操作
第\(i\)次操作會有\(2^{\operatorname{lg}N-i+1}\)個長度為\(2^{i-1}\)子陣列，將他們兩兩合併

Code結構

void SORT(int *a, int N);
void CONQUER(pii y0_range, pii y1_range);

void FFT(int *a, int N) {
    
    SORT(a, N);

    for (int w = 1; w < N; w <<= 1) {
        // w: the width of subarrays that hasn't been conquered

        for (int s = 0; s < N; s += w * 2) {
            // s: the start position of the first array to be conquered

            CONQUER({s, s + w}, {s + w, s + w * 2});
        }
    }
}

SORT

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

SORT

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

SORT

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

SORT

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

把index反過來排

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

\(000\)

\(001\)

\(010\)

\(011\)

\(100\)

\(101\)

\(110\)

\(111\)

速度優化

不喜歡\(\operatorname{cis}\)

相信我們都很會做普通的 CONQUER 了

0,0	0,1

1,0	1,1

0	1	2	3

相信我們都很會做普通的 CONQUER 了

0,0	0,1

1,0	1,1

0	1	2	3

一次做一堆 CONQUERS ?

0,0	0,1

1,0	1,1

0	1	2	3

0,0	0,1

1,0	1,1

0	1	2	3

所有紅色的地方都要搭配\(\omega_4^0\)

所有綠色的地方都要搭配\(\omega_4^1\)

0,0	0,1

1,0	1,1

0	1	2	3

0,0	0,1

1,0	1,1

0	1	2	3

比起分開 Conquer

不如依序把\(\omega\)丟進去算！

0,0	0,1

1,0	1,1

0	1	2	3

0,0	0,1

1,0	1,1

0	1	2	3

Code結構

奇怪排序
做\(\operatorname{lg}N\)次操作
第\(i\)次操作會有\(2^{\operatorname{lg}N-i+1}\)個長度為\(2^{i-1}\)的子陣列，將他們兩兩合併

Code結構

奇怪排序
做\(\operatorname{lg}N\)次操作
第\(i\)次操作會有\(2^{\operatorname{lg}N-i+1}\)個長度為\(2^{i-1}\)的子陣列，將他們兩兩合併

\(\Rightarrow\) 每次操作中的每一組子陣列的長度都相同，用剛剛的方式壓縮\(\omega\)的運算量！

Code結構

cd OMEGA(int n, int k);
void SORT(cd *a, int N);
// void CONQUER(pii y0_range, pii y1_range);
void SMALL_CONQUER(cd &L, cd &R, cd omega);

void FFT(cd *a, int N) {
    
    SORT(a, int N);

    for (int w = 1; w < N; w <<= 1) {
        // w: the width of subarrays that hasn't been conquered

        int omega_n = w * 2;

        for (int omega_k = 0; omega_k < w; omega_k++) {
            cd omega = OMEGA(omega_n, omega_k);

            for (int s = 0; s < N; s += w * 2) {
                // s: the start position of the first array to be conquered

                SMALL_CONQUER(a[s + omega_k], a[s + w + omega_k], omega);

            }
        }
    }
}

Code

typedef long long ll;
typedef pair<int, int> pii;
typedef complex<double> cd;

inline cd cis(double theta) {
    return cd(cos(theta), sin(theta));
}

inline cd OMEGA(int n, int k) {
    return cis(acos(-1) * 2 * k / n);
}

inline int REVERSE(int x, int N) {
    int ans = 0;
    for (int i = 1; i < N; i <<= 1) {
        ans <<= 1;
        if (i & x) ans |= 1;
    }
    return ans;
}

void FFT(cd *a, int N) {
    for (int i = 0; i < N; i++) {
        int r = REVERSE(i, N);
        if (i < r) swap(a[i], a[r]);
    }
    /* processing */
    for (int w = 1; w < N; w <<= 1) {
        int omega_n = w * 2;
        for (int omega_k = 0; omega_k < w; omega_k++) {
            cd omega = OMEGA(omega_n, omega_k);
            for (int s = 0; s < N; s += w << 1) {
                cd &L = a[s + omega_k], &R = a[s + omega_k + w];
                cd l = L, r = R * omega;
                L = l + r;
                R = l - r;
            }
        }
    }
}

inverse FFT

大統整

\(A(x),B(x)\)

\(\alpha,\beta\)

\(c\)

\((A*B)(x)\)

\(\operatorname{DFT}\)

兩兩相乘

\(\operatorname{DFT}^{-1}\)

\bold{\Omega}_N=\left[\begin{array}{ccccc}\omega_N^0&\omega_N^0&\omega_N^0&\cdots&\omega_N^0\\\omega_N^0&\omega_N^1&\omega_N^2&\cdots&\omega_N^{N-1}\\\omega_N^0&\omega_N^2&\omega_N^4&\cdots&\omega_N^{2(N-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\\omega_N^0&\omega_N^{N-1}&\omega_N^{2(N-1)}&\cdots&\omega_N^{(N-1)(N-1)}\end{array}\right]

\bold{\Omega}_N^{-1}=\left\lbrack\begin{array}{ccccc}\omega_N^0&\omega_N^0&\omega_N^0&\cdots&\omega_N^0\\\omega_N^0&\omega_N^{-1}&\omega_N^{-2}&\cdots&\omega_N^{-(N-1)}\\\omega_N^0&\omega_N^{-2}&\omega_N^{-4}&\cdots&\omega_N^{-2(N-1)}\\\vdots&\vdots&\vdots&\ddots&\vdots\\\omega_N^0&\omega_N^{-(N-1)}&\omega_N^{-2(N-1)}&\cdots&\omega_N^{-(N-1)(N-1)}\end{array}\right]\cdot\frac{1}{N}

時間不夠了

請讓我們省略證明

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&y_{0,k}-\omega_N^k y_{1,k} \end{array}

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&y_{0,k}-\omega_N^k y_{1,k} \end{array}

\begin{array}{cl}&\bold{A}(\omega_N^{-k})=y_k\\ =&y_{0,k}+\omega_N^{-k}y_{1,k}\\ \end{array}

\begin{array}{cl}&\bold{A}(\omega_N^{-(k+n)})=y_{k+n}\\ =&y_{0,k}-\omega_N^{-k}y_{1,k} \end{array}

if (inverse) ...

\(\omega_N^k\longrightarrow\omega_N^{-k}\)
最後記得每個元素\(\times\frac{1}{N}\)

Code

typedef long long ll;
typedef pair<int, int> pii;
typedef complex<double> cd;

inline cd cis(double theta) {
    return cd(cos(theta), sin(theta));
}

inline cd OMEGA(int n, int k) {
    return cis(acos(-1) * 2 * k / n);
}

inline int REVERSE(int x, int N) {
    int ans = 0;
    for (int i = 1; i < N; i <<= 1) {
        ans <<= 1;
        if (i & x) ans |= 1;
    }
    return ans;
}

void FFT(cd *a, int N, bool inv) {
    for (int i = 0; i < N; i++) {
        int r = REVERSE(i, N);
        if (i < r) swap(a[i], a[r]);
    }
    /* processing */
    for (int w = 1; w < N; w <<= 1) {
        int omega_n = w * 2;
        for (int omega_k = 0; omega_k < w; omega_k++) {
            cd omega = OMEGA(omega_n, (inv ? -1 : 1) * omega_k);
            for (int s = 0; s < N; s += w << 1) {
                cd &L = a[s + omega_k], &R = a[s + omega_k + w];
                cd l = L, r = R * omega;
                L = l + r;
                R = l - r;
            }
        }
    }
    if (inv) for (int i = 0; i < N; i++) a[i] /= N;
}

NTT

為什麼FFT複雜度是好的

\operatorname{DFT}(\bold{A})=\bold{Y}=\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\y_4\\y_5\\\vdots\\y_{N-2}\\y_{N-1}\end{array}\right]

\begin{array}{cl}&\bold{A}(\omega_N^k)=y_k\\ =&y_{0,k}+\omega_N^ky_{1,k}\\ \end{array}

（若\(0\leq k<n\)）

\begin{array}{cl}&\bold{A}(\omega_N^{k+n})=y_{k+n}\\ =&\bold{A}_0(\omega_N^{2(k+n)})+\omega_N^{k+n}\bold{A}_1(\omega_N^{2(k+n)}) \end{array}

\omega_N^{2(k+n)}=\omega_N^{2k+N}=\omega_N^{2k}

替代方案

誠徵：一個定義\(\omega_N^k\)的方式使得

\(\omega_N^{k_1}\omega_N^{k_2}=\omega_N^{k_1+k_2}\)

\(\omega_N^{k+N}=\omega_N^k\)

\(\omega_N^0,\omega_N^1,\dots\omega_N^{N-1}\)皆不相同

直接破梗

剩餘系

#define int long long
struct Rest {
    int val;
    const static int mod;
    // mod must be a prime
    Rest() {
        val = 0;
    }
    Rest(int _x) {
        val = _x % mod;
    }
    Rest operator*(Rest other) {
        return Rest(val * other.val % mod);
    }
    Rest pow(int x) {
        Rest now(1);
        for (int i = 1 << 30; i > 0; i >>= 1) {
            now = now * now;
            if (i & x) now = operator*(now);
        }
        return now;
    }
    Rest inv() {
        return pow(mod - 2);
    }
}

剩餘系

令\(\omega_6^k=3^k\%7\)

k	0	1	2	3	4	5
omega	1	3	2	6	4	5

剩餘系

令\(\omega_6^k=3^k\%7\)

k	0	1	2	3	4	5
omega	1	3	2	6	4	5

我們要的是\(N=2^x\)

剩餘系

令\(\omega_{8388608}^k=15311432^k\%998244353\)

我們要的是\(N=2^x\)

剩餘系

令\(\omega_{2^{23}}^k=(3^{119})^k\%998244353\)

我們要的是\(N=2^x\)

剩餘系

令\(\omega_{2^{23}}^k=(3^{119})^k\%998244353\)

我們要的是\(N=2^x\)

它是好的

剩餘系

令\(\omega_{2^{23}}^k=(3^{119})^k\%998244353\)

喔我現在只要用\(\omega_4^k\)

它是好的

剩餘系

令\(\omega_{2^{23}}^k=(3^{119})^k\%998244353\)

喔我現在只要用\(\omega_4^k\)

\(\omega_4^k=\omega_{4\cdot2^{21}}^{k\cdot2^{21}}\)

剩餘系

令\(\omega_{2^{23}}^k=(3^{119})^k\%998244353\)

當\(N\leq2^{23}\)，且運算後的數字不超過\(998244353\)時可以使用這套系統

FFT\(\longrightarrow\)NTT

只是換掉\(\omega\)的系統而已
沒有浮點數誤差
值域有點小\(\longrightarrow\)中國剩餘

Code

#define int long long
struct Rest {
    int val, mod;
    Rest() {
        val = 0;
        mod = 0;
    }
    Rest(int _val, int _mod) {
        val = _val % _mod;
        mod = _mod;
    }
    Rest operator+(Rest other) {
        return Rest((val + other.val) % mod, mod);
    }
    Rest operator-(Rest other) {
        return Rest((val - other.val + mod) % mod, mod);
    }
    Rest operator*(Rest other) {
        return Rest(val * other.val % mod, mod);
    }
    Rest pow(int x) {
        Rest ans(1, mod);
        for (int i = 1 << 30; i > 0; i >>= 1) {
            ans = ans * ans;
            if (i & x) ans = operator*(ans);
        }
        return ans;
    }
    Rest inv() {
        return pow(mod - 2);
    }
};

int REVERSE(int id, int N) {
    int ans = 0;
    for (int i = 1; i < N; i <<= 1) {
        ans <<= 1;
        if (i & id) ans |= 1;
    }
    return ans;
}

struct NTTBASE {
    Rest g, omegaB, omegaB_1;
    int mul, upLG;
    NTTBASE() {
        g = Rest();
        mul = 0;
        omegaB = Rest();
        omegaB_1 = Rest();
    }
    NTTBASE(int _g, int _mul, int _upLG, int _mod) {
        g = Rest(_g, _mod);
        mul = _mul;
        upLG = _upLG;
        omegaB = g.pow(mul);
        omegaB_1 = omegaB.inv();
    }
    Rest OMEGA(int _lgn, int _k) {
        Rest ans = omegaB;
        for (int i = _lgn; i < upLG; i++) {
            ans = ans * ans;
        }
        return ans.pow(_k);
    }
    void operator()(Rest *a, int N, bool inv) {
        for (int i = 0; i < N; i++) {
            int rev = REVERSE(i, N);
            if (i < rev) swap(a[i], a[rev]);
        }
        for (int w = 1, wlg = 0; w < N; w <<= 1, wlg++) {
            Rest omegaBase = OMEGA(wlg + 1, 1), omega(1, g.mod);
            if (inv) omegaBase = omegaBase.inv();
            for (int omega_k = 0; omega_k < w; omega_k++) {
                for (int s = 0; s < N; s += w * 2) {
                    Rest &L = a[s + omega_k], &R = a[s + w + omega_k];
                    Rest l = L, r = R * omega;
                    L = l + r;
                    R = l - r;
                }
                omega = omega * omegaBase;
            }
        }
        if (inv) {
            Rest temp = Rest(N, g.mod).inv();
            for (int i = 0; i < N; i++) a[i] = a[i] * temp;
        }
    }
};

NTTBASE NTT(3, 119, 23, MOD);

一些題目

Apples and Bananas

\bold{A}=\left[\begin{array}{ccccc}0&2&4&6&8\end{array}\right]\\ \bold{B}=\left[\begin{array}{ccccc}1&3&5&7&9\end{array}\right]\\ \bold{C}=\bold{A}*\bold{B}=\left[\begin{array}{ccccc}c_0&c_1&c_2&\cdots&c_8\end{array}\right]\\ c_3=?