目錄
寫在前面
- 由於上一篇總結的版面限制，特開此文來記錄 \(\text{OI}\) 中多項式類數學相關的問題。
- 該文啟發於Miskcoo的博客，甚至一些地方直接引用，在此特別說明；若文章中出現錯誤，煩請告知。
- 感謝你的造訪。
前置技能
多項式相關
形同 \(P(X)=a_0+a_1X+a_2X^2+\cdots+a_nX^n\) 的形式冪級數 \(P(X)\) 稱為多項式。其中 \(\{a_i|i\in[0,n]\}\) 為多項式的系數； \(n\) 表示多項式的次數。
多項式的系數表示
對於 \(n\) 次多項式 \(P(X)\) 的系數 \(\{a_i|i\in[0,n]\}\) ，我們記向量 \(\vec a=(a_0,a_1,\cdots,a_n)\) 為 \(P(x)\) 的系數表示。
多項式的點值表示
對於一個 \(n\) 次多項式 \(P(X)\) ，我們由函數和方程的思想，對於函數 \(P(X)\) 的圖像，我們只要確定了該圖像上的 \(n+1\) 個點，那么 \(P(X)\) 是確定且唯一的。我們記這 \(n+1\) 個點組成的集合為 \(\{\left(x_i, P(x_i)\right)| i\in[0,n]\}\) 。那么該集合為多項式 \(P(x)\) 的點值表示。
一個多項式有不同的點值表示，但任意一個點值都能找出其唯一對應的多項式。
多項式的系數表示和點值表示能夠互相轉換。
復數相關
記 \(i=\sqrt{-1}\) 。我們把形如 \(a+bi\) （ \(a,b\) 均為實數）的數稱為復數。其中 \(a\) 為虛數的實部， \(b\) 為虛數的虛部。
復數的意義
對於任意一個復數 \(a+bi\) ，都能用復平面上的唯一一個向量 \((a,b)\) 來表示。
復平面：即復數平面，由實軸作為 \(x\) 軸，虛軸作為 \(y\) 軸構成。
下文中出現的向量 \((a,b)\) 來代指虛數 \(a+bi\) 。
共軛復數：對於復數 \(z=a+bi\) ，對於另外一個復數 \(z'=a-bi\) ，我們稱 \(z'\) 為 \(z\) 的共軛復數，記做 \(\overline{z}\) 。容易發現，一個復數與其共軛復數的實部相等，虛部互為相反數。
復數的輻角：我們可以將復數 \(z\) 寫成 \(z=r\times(\cos\theta+i\sin\theta)\) ，其中 \(r\) 為復數 \(z\) 的模長， \(\theta\) 為復數 \(z\) 的輻角。一個復數有多個輻角，這些值相差 \(2\pi\) 。我們將 \(\theta\in[-\pi,\pi)\) 的 \(\theta\) 叫做輻角的主值。指數形式： \(z=r\times(\cos\theta+i\sin\theta)=re^{i\theta}\) 。
復數的基本運算
對於復數 \((a,b)\) 和 \((c,d)\) 。
滿足加法法則 \((a,b)\pm(c,d)=(a\pm c,c\pm d)\) 。
滿足乘法法則 \((a,b)\times(c,d)=(ac-bd,bc+ad)\) 。
對於除法，只需分子分母同乘分母的共軛復數，將分母實數化，分子做復數乘法即可。
復數乘法的幾何意義：模長相乘，幅角相加。
證明：
對於兩個復數 \(z_1=r_1\times(\cos\theta_1+i\sin\theta_1),z_2=r_2\times(\cos\theta_2+i\sin\theta_2)\) 。
\[\begin{aligned}z_1\times z_2&=r_1\times(\cos\theta_1+i\sin\theta_1)\times r_2\times(\cos\theta_2+i\sin\theta_2)\\&=r_1r_2\times((\cos\theta_1\cos\theta_2-\sin\theta_1\sin\theta_2)+i(\sin\theta_1\cos\theta_2+\cos\theta_1\sin\theta_2))\\&=r_1r_2\times(\cos(\theta_1+\theta_2)+i\sin(\theta_1+\theta_2))\end{aligned}\]
得證。
單位根
\(n\) 次單位根是指能夠滿足方程 \(z^n=1\) 的復數，這些復數一共有 \(n\) 個它們都分布在復平面的單位圓上，並且構成一個正 \(n\) 邊形，它們把單位圓等分成 \(n\) 個部分。
根據復數乘法相當於模長相乘，幅角相加就可以知道， \(n\) 次單位根的模長一定是 \(1\) ，幅角的 \(n\) 倍是 \(0\) 。
這樣， \(n\) 次單位根也就是
\[e^{\frac{2\pi ki}{n}}, k = 0, 1, 2, \cdots, n - 1\]
再根據歐拉公式
\[e^{\theta i}=\cos\theta + i\sin\theta\]
就可以知道 \(n\) 次單位根的算術表示
如果記 \(\omega_n=e^{\frac{2\pi i}{n}}\) ，那么 \(n\) 次單位根就是 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 。
代碼相關
給出此文出現的代碼中的一些宏定義
#define dob complex<double>
const double pi = acos(-1.0);
const int mod = 998244353;多項式乘法
給定兩個多項式 \(A(x),B(x)\)
\[A(x) = \sum_{i=0}^na_ix^i = a_nx^n+a_{n-1}x^{n-1}+\cdots+a_1x+a_0 \\ B(x) = \sum_{i=0}^nb_ix^i = b_nx^n+b_{n-1}x^{n-1}+\cdots+b_1x+b_0\]
將這兩個多項式相乘得到 \(C(x) = \sum_{i=0}^{2n}c_ix^i\) ，在這里
\[c_i=\sum_{j+k=i,0\leq j,k\leq n}a_jb_k\]
如果一個個去算 \(c_i\) 的話，要花費 \(O(n^2)\) 的時間才可以完成，但是，這是在系數表示下計算的，如果轉換成點值表示，知道了 \(A(x),B(x)\) 的點值表示后，由於點數是 \(O(n)\) ，就可以直接將其相乘，在 \(O(n)\) 的時間內得到 \(C(x)\) 的點值表示。
由於 \(C(x)\) 的次數為 \(2n\) ，所以我們可以在 \(A(x),B(x)\) 上取 \(2n+1\) 個點，便於唯一確定 \(C(x)\) 。
如果能夠找到一種有效的方法幫助我們在多項式的點值表示和系數表示之間轉換，我們就可以快速地計算多項式的乘法了，快速傅里葉變換就可以做到這一點。
快速傅里葉變換
由剛才的設想，我們只需按序進行下述三次操作，即可降低多項式乘法的復雜度：
- 將 \(A(x),B(x)\) 的系數表達轉為點值表達；
- 將 \(A(x),B(x)\) 的點值表達相乘；
- 將得到的結果轉為系數表達。
我們已知 2. 是 \(O(n)\) 。剩下的 1.和3. 可以用快速傅里葉變換來實現 \(O(n\log_2 n)\) 的變換。
DFT
\(\text{DFT}\) 是實現上述的 1. 過程的。它是一個基於分治策略的算法。
具體的思想是將多項式的 \(n\) 個系數通過變換變成 \(n\) 個點值。
對於一個多項式 \(A(x)=\sum_{i=0}^{n-1}a_ix^i\) ，我們先將 \(n\times 2\) ，理由上面說了，因為為了確定 \(C(x)\) 需要多取一些點。為了方便之后的處理，我們繼續將 \(n\) 增大使得 \(n=2^m\) ，找到這個最小的 \(n\) ，並將不足的項的系數變為 \(0\) 。
接着將 \(n\) 個 \(n\) 次單位根 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 代入 \(A(x)\)
\[A(\omega_n^k) = \sum_{i=0}^{n-1}a_i\omega_n^{ki} , k\in[0,n)\]
我們發現，這樣計算出多項式 \(A(x)\) 的點值表示 \(\{\left(\omega_n^k, A(\omega_n^k)\right)| k\in[0,n)\}\) 的復雜度仍是 \(O(n^2)\) 的，如何將這個復雜度優化就是整個算法的關鍵。具體優化，就要利用單位根的性質了。
第一步，是將點集中每一項按照指數的奇偶分類：
\[\begin{eqnarray*} A(\omega_n^k) &=& \sum_{i=0}^{n-1}a_i\omega_n^{ki} \\ &=& \sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_n^{2ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_n^{2ki} \\ \end{eqnarray*}\]
但是這樣分類后沒什么用，因為對於每個系數 \(a_i\) ，它被計算的次數依舊是 \(n\) 次，因為對於每個 \(k\) 都要與 \(a_i\) 乘一次。
我們在回到上面這個式子，試着將它變形
\[A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\left(\omega_n^{ki}\right)^2+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\left(\omega_n^{ki}\right)^2\]
注意到的是
\[\omega_n^2=\left(e^{\frac{2\pi i}{n}}\right)^2=e^{\frac{2\pi i}{\frac{n}{2}}}=\omega_{\frac{n}{2}}\]
這個等式可以直接由上述式子推出，但是我們可以去直觀感受一下為什么會這樣。
對於一個復數 \(z\) ，它的平方 \(z^2\) 依舊滿足復數乘法運算的規律：模長相乘，幅角相加。
那么對於一個單位根 \(\omega_n=e^{\frac{2\pi i}{n}}\) ，它的平方模長不變，輻角翻倍。
我們想到一個復數有多個輻角，並且差為 \(2\pi\) 。既然如此，那么對於 \(n\) 次單位根 \(\omega_n\) ，與另一個 \(n\) 次單位根 \(\omega_n^{\frac{n}{2}}\) 輻角差 \(\pi\) 。平方后輻角翻倍那么相差 \(2\pi\) ，所以說這兩個 \(n\) 次單位根的平方是相同的。
同時我們可以得到
\[\omega_n^{\frac{n}{2}+k} = \omega_n^{\frac{n}{2}}\cdot \omega_n^k = -\omega_n^k\]
因為 \(\omega_n^{\frac{n}{2}}=e^{\pi i}=-1\) ，得證。
即原式變成了
\[A(\omega_n^k)=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\]
不過這時是 \(k<\frac{n}{2}\) 時才成立的。由之前得到的式子 \(\omega_n^{\frac{n}{2}+k} = \omega_n^{\frac{n}{2}}\cdot \omega_n^k = -\omega_n^k\) 其余部分則應該滿足：
\[\begin{eqnarray*} A(\omega_n^{k+\frac{n}{2}}) &=& \sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^{k+\frac{n}{2}}\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \\ &=&\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki} \end{eqnarray*}\]
這樣對於每個 \(a_i\) ，代入的值變成了 \(1, \omega_{\frac{n}{2}}, \omega_{\frac{n}{2}}^2, \cdots, \omega_{\frac{n}{2}}^{\frac{n}{2}-1}\) ，問題變成了兩個規模減半的子問題，只要遞歸下去計算就可以了，復雜度是 \(O(n\log_2 n)\) 。
IDFT
剛才說了， \(\text{DFT}\) 是將系數表示轉為點值表示。而我們現在需要解決的是將一個多項式的點值表示轉為系數表示。 \(\text{IDFT}\) 則可實現這一過程。 \(\text{IDFT}\) 是 \(\text{DFT}\) 的逆過程。
考慮到一個本質的問題，如何基礎地將點值表示轉為系數表示，顯然就是解如下的一個方程組：
\[\begin{equation*} \left\{ \begin{array}{ccccccccc} a_0(\omega_n^0)^{0}&+&\cdots&+&a_{n-2}(\omega_n^0)^{n-2}&+&a_{n-1}(\omega_n^0)^{n-1}&=&A(\omega_n^0) \\ a_0(\omega_n^1)^{0}&+&\cdots&+&a_{n-2}(\omega_n^1)^{n-2}&+&a_{n-1}(\omega_n^1)^{n-1}&=&A(\omega_n^1) \\ \vdots & & \vdots & &\vdots& & \vdots & & \vdots\\ a_0(\omega_n^{n-1})^{0}&+&\cdots&+&a_{n-2}(\omega_n^{n-1})^{n-2}&+&a_{n-1}(\omega_n^{n-1})^{n-1}&=&A(\omega_n^{n-1}) \end{array} \right. \end{equation*}\]
寫成矩乘形式就是：
\[\begin{equation} \begin{bmatrix} (\omega_n^0)^0 & (\omega_n^0)^1 & \cdots & (\omega_n^0)^{n-1} \\ (\omega_n^1)^0 & (\omega_n^1)^1 & \cdots & (\omega_n^1)^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{n-1})^0 & (\omega_n^{n-1})^1 & \cdots & (\omega_n^{n-1})^{n-1} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation}\]
左上的第一個矩陣為 \(\mathbf V\) 。現在考慮構造下面這個矩陣 \(d_{ij}=\omega_n^{-ij}\)
\[\begin{equation*} \mathbf D = \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \end{equation*}\]
記 \(\mathbf E=\mathbf D \cdot \mathbf V\) ，容易發現
\[\begin{eqnarray*} e_{ij} &=& \sum_{k=0}^{n-1} d_{ik} v_{kj} \\ &=& \sum_{k=0}^{n-1} \omega_n^{-ik}\omega_n^{kj} \\ &=& \sum_{k=0}^{n-1} \omega_n^{k(j-i)} \end{eqnarray*}\]
而對於這個式子容易發現
- 當 \(i=j\) 時， \(\omega_n^{k(j-i)}=\omega_n^{0}=1\) ，故 \(e_{ij}=\sum_{k=0}^{n-1}1=n\) 。
- 當 \(i\neq j\) 時， \(e_{ij}=\sum_{k=0}^{n-1}\omega_n^{k(j-i)}=\frac{\omega_n^0(1-\omega_n^{n(j-i)})}{1-\omega_n^{j-i}}\) ，由於 \(\omega_n^n=1\) ，故 \(e_{ij}=0\) 。
所以單位矩陣 \(\mathbf I_n=\frac{1}{n}\mathbf E\) ，由於 \(\mathbf I_n=\mathbf V\cdot\mathbf V^{-1},\mathbf E=\mathbf D \cdot \mathbf V\) ， \(\mathbf V\cdot\mathbf V^{-1}=\frac{1}{n}\mathbf D \cdot \mathbf V\) ， \(\mathbf V^{-1}=\frac{1}{n}\mathbf D\) 。其中 \(\mathbf V^{-1}\) 為 \(\mathbf V\) 的逆矩陣。
在上述 \((1)\) 式中等式左右兩邊左乘 \(\frac{1}{n}\mathbf D\) 。能得到
\[\begin{equation*} \begin{bmatrix} a_0 \\ a_1 \\ \vdots \\ a_{n-1} \end{bmatrix} = \frac{1}{n} \begin{bmatrix} (\omega_n^{-0})^0 & (\omega_n^{-0})^1 & \cdots & (\omega_n^{-0})^{n-1} \\ (\omega_n^{-1})^0 & (\omega_n^{-1})^1 & \cdots & (\omega_n^{-1})^{n-1} \\ \vdots & \vdots & \ddots & \vdots \\ (\omega_n^{-(n-1)})^0 & (\omega_n^{-(n-1)})^1 & \cdots & (\omega_n^{-(n-1)})^{n-1} \end{bmatrix} \begin{bmatrix} A(\omega_n^0) \\ A(\omega_n^1) \\ \vdots \\ A(\omega_n^{n-1}) \end{bmatrix} \end{equation*}\]
即 \(a_k=\frac{1}{n}\sum_{i=0}^{n-1}=\omega_n^{-ki}A(\omega_n^i)\) ，這樣， \(\text{IDFT}\) 就相當於把 \(\text{DFT}\) 過程中的 \(\omega_n^i\) 換成 \(\omega_n^{-i}\) ，然后做一次 \(\text{DFT}\) ，之后結果除以 \(n\) 就可以了。
算法實現
遞歸實現
由上面的說明想必很容易模擬上述過程，來實現遞歸的 \(\text{DFT}\) 和 \(\text{IDFT}\) 。
為了方便閱讀，我們在這重寫上面的式子，對於 \(\text{DFT}\) ：
\[\begin{aligned}A(\omega_n^k)&=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}+\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\\ A(\omega_n^{k+\frac{n}{2}}) &=\sum_{i=0}^{\frac{n}{2}-1}a_{2i}\omega_{\frac{n}{2}}^{ki}-\omega_n^k\sum_{i=0}^{\frac{n}{2}-1}a_{2i+1}\omega_{\frac{n}{2}}^{ki}\end{aligned}\]
void FFT(dob *A, int len, int o) {
// len 為當前遞歸區間長度； o 為識別因子，若 o=1 表示進行 DFT ，若為 -1 表示進行 IDFT
if(len == 1) return;
dob wn(cos(2.0*pi/len), sin(2.0*pi*o/len)), w(1,0), t;
//注意此處 wn 的初值
dob A0[len>>1], A1[len>>1];
for (int i = 0; i < (len>>1); i++) A0[i] = A[i<<1], A1[i] = A[i<<1|1];
FFT(A0, len>>1, o); FFT(A1, len>>1, o);
for(int i = 0; i < (len>>1); i++, w *= wn) {
t = w*A1[i];
A[i] = A0[i]+t;
A[i+(len>>1)] = A0[i]-t;
}
}遞歸來實現 \(\text{FFT}\) 的顯著的優點就是直觀簡潔，幾乎就是按照式子模擬即可；不過缺點是常數巨大，在實戰上毫不適用。
迭代實現
假設現在有 \(16\) 個數要進行 \(\text{DFT}\) 來看看遞歸的過程。
（圖片轉自Miskcoo）
在 \(\text{Step1} \rightarrow \text{Step2}\) 的過程中，按照奇偶分類，二進制位中最后一位相同的被分到同一組；
在 \(\text{Step2} \rightarrow \text{Step3}\) 的過程中，仍然按照奇偶，只不過不是按照數字的奇偶性，而是下標的奇偶性，二進制位中最后兩位相同的才被分到同一組；
在 \(\text{Step3} \rightarrow \text{Step4}\) 的過程中，二進制位中最后三位相同的數字被分在同一組；
我們將所有數的二進制翻轉，容易發現每次分在同一組內的都是前綴相同的，並且同一組內的數值是連續的。
由於迭代實現，我們可以先將原來的 \(A\) 數組的位置預處理成 \(B\) 數組——遞歸最下面一層（最后一步）的 \(A\) 的位置。
假設 reverse(i) 是將二進制位反轉的操作，那么 \(A\) 和 \(B\) 之間有這樣的關系 B[reverse(i)]=A[i] ，也就是說， B[i+1]=A[reverse(reverse(i)+1)] ， \(B\) 中第 \(i\) 項的下一項就是將 \(i\) 反轉后加 \(1\) 再反轉回來 \(A\) 中的那一項，所以現在要模擬的就是從高位開始的二進制加法。
我們可以處理一個數組 \(R\) ： R[i]=reverse(i) ，即 B[R[i]]=A[i] 。
對於 \(R\) 的預處理，可以 \(O(n)\) 遞推出，其中 \(L\) 表示 \(n=2^L\) ：
for (int i = 0; i < n; i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));顯然滿足 R[R[i]]=i ，所以對於一組 \(i,R_i\) ，只要交換一次位置即可。
既然已經處理最后的位置 \(B\) 。那么我們就可以枚舉區間長度迭代計算了。
void FFT(dob *A, int o) {
for (int i = 0; i < n; i++) if (i < R[i]) swap(A[i], A[R[i]]);
for (int i = 1; i < n; i <<= 1) {
//枚舉區間長度
dob wn(cos(pi/i), sin(pi*o/i)), x, y;
for(int j = 0; j < n; j += (i<<1)) {
//枚舉區間左端點
dob w(1, 0);
for(int k = 0; k < i; k++, w *= wn) {
//枚舉 k
x = A[j+k]; y = w*A[j+i+k];
A[j+k] = x+y;
A[j+k+i] = x-y;
}
}
}
}快速數論變換
有些題目要求答案要對一個質數取模，用 \(\text{FFT}\) 難免會有缺陷了；畢竟在復數域計算難免有精度損失。
首先來看 \(\text{FFT}\) 中能在 \(O(n\log_2n)\) 時間內變換用到了單位根 \(\omega\) 的什么性質。
- \(\omega_n^n=1\) ，有這樣一個初值，方便后面的計算；
- \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 是互不相同的，這樣帶入計算出來的點值才可以用來還原出系數；
- \(\omega_n^2=\omega_{\frac{n}{2}}, \omega_n^{\frac{n}{2}+k}=-\omega_n^k\) ，這使得在按照指數奇偶分類之后能夠把帶入的值也減半使得問題規模能夠減半。
- \[\sum_{k=0}^{n-1} (\omega_n^{j-i})^k = \begin{eqnarray*} \left\{ \begin{aligned}0, ~~~&i \neq j\\ n, ~~~&i = j \end{aligned} \right. \end{eqnarray*}\]
這樣保證了能夠使用相同的方法進行逆變換得到系數表示。
如果我們能夠在數論域上找到這樣的類似單位根的東西。就可以用相同的思想來進行簡化運算。
這樣我們引出了原根的概念。
原根
根據費馬小定理我們知道，對於一個素數 \(p\) ，有下面這樣的關系
\[a^{p-1} \equiv 1 \pmod p\]
這樣類比單位根，可以滿足“周期性”。
對於一個數 \(g\) 滿足 \(g^0, g^1, \cdots, g^{p-2} \pmod p\) 互不相同，那么稱 \(g\) 是 \(p\) 的原根。
令 \(n=2^k\) ，我們取素數 \(p = r\cdot 2^k + 1\) （滿足這樣形式的素數叫做費馬素數，朴素的 \(\text{NTT}\) 是只在費馬素數下適用的。），然后取 \(p\) 的原根 \(g\) ，然后我們令 \(g_n\equiv g^r\pmod p\) ，這樣就能滿足 \(g_n^0, g_n^1, \cdots, g_n^{n-1} \pmod p\) 互不相同，並且 \(g_n^n\equiv 1\pmod p\) 。
首先 \(g_n^n\equiv 1\pmod p\) 是顯然的，因為 \(2^{rn}\geq r\cdot 2^n\) ，就會滿足 \(r\cdot 2^n\mid 2^{r+n}\) ；
而對於 \(g_n^0, g_n^1, \cdots, g_n^{n-1} \pmod p\) 互不相同的證明，其實等價於證明 \(0r,1r,\cdots,(n-1)r \pmod{p-1}\) 互不相同，這個在之前的博客中有講，就不贅述了。
於是這樣就滿足了上面所說的 1.2. 。
由於 \(p\) 是質數且 \(g_n^n \equiv 1 \pmod p\) ，所以 \(g_n^{\frac{n}{2}} \equiv \pm 1 \pmod p\) ，又由於 2. ，所以 \(g_n^{\frac{n}{2}} \equiv -1 \pmod p\) 。直接滿足了性質 3. 。
對於 4. ，可以用相同的方法代入計算。不再贅述。
我們這樣就找到了這樣一系列滿足條件的數 \(g_n\) 。
對於一個模數，需要去找到它的原根，似乎比較麻煩，但Miskcoo提供了一些常用的數值。查表就好了。
算法實現
實現和 \(\text{FFT}\) 類似，我們依舊用迭代實現。
void NTT(int *A, int o) {
for (int i = 0; i < n; i++) if (i < R[i]) swap(A[i], A[R[i]]);
for (int i = 1; i < n; i <<= 1) {
//枚舉區間長度
int gn = quick_pow(3, (mod-1)/(i<<1)), x, y;
if (o == -1) gn = quick_pow(gn, mod-2);
for (int j = 0; j < n; j += (i<<1)) {
//枚舉區間左端點
int g = 1;
for (int k = 0; k < i; k++, g = 1ll*g*gn%mod) {
//枚舉 k
x = A[j+k], y = 1ll*g*A[j+k+i]%mod;
A[j+k] = (x+y)%mod;
A[j+k+i] = (x-y+mod)%mod;
}
}
}
}模數任意的解決方案
傳說中的 \(\text{MTT}\) （快速毛爺爺變換），具體思想是取幾個乘積大於 \(n(mod-1)^2\) 的費馬素數作為模數，分別求出結果后用 \(crt\) 合並就好了。
應用
快速卷積
對於兩個定義在 \(\mathbb{N}\) 上的函數 \(f(n),g(n)\) ，定義 \(f\) 和 \(g\) 卷積為 \(f\otimes g\)
\[(f \otimes g)(n) = \sum_{i=0}^n f(i)g(n-i)\]
容易發現兩個數論函數的卷積其實就是兩個多項式的乘積。
可以用 \(\text{FFT}\) 或 \(\text{NTT}\) 優化。
更一般地，對於 \(h(k) = \sum_{i=0}^n f(i)g(i+k)\) ，我們可以設 \(f'(x)=f(n-x)\) ，顯然 \(h(k) = \sum_{i=0}^n f'(n-i)g(i+k)\) ，也是一個卷積的形式。
多項式求逆
基本概念
對於多項式 \(A(x),B(x)\) ，存在唯一的 \(Q(x),R(x)\) 滿足 \(A(x)=B(x)Q(x)+R(x)\) ，其中 \(R\) 的次數小於 \(B\) 的次數，我們稱 \(Q(x)\) 為 \(B(x)\) 除 \(A(x)\) 的商， \(R(x)\) 為 \(B(x)\) 除 \(A(x)\) 的余數，可以記作
\[A(x) \equiv R(x) \pmod {B(x)}\]
對於一個多項式 \(A(x)\) ，如果存在 \(B(x)\) 滿足 \(B\) 的次數小於等於 \(A\) 的次數，並且
\[A(x)B(x) \equiv 1 \pmod {x^n}\]
那么稱 \(B(x)\) 為 \(A(x)\) 在 \(\mod x^n\) 意義下的逆元，記作 \(A^{−1}(x)\) 。
對於多項式 \(A(x)\) ，對於 \(A(x)\mod x^n\) 直觀的解釋是提出 \(A(x)\) 中的 \(x^0\sim x^{n-1}\) 次項。
求解方法
考慮如何求 \(A^{-1}(x)\) 。
- 當 \(n=1\) 時，容易發現 \(A(x)\equiv c\pmod x\) 此時 \(c\) 是一個常數，故 \(A^{-1}(x)\) 也是一個常數，即 \(c^{-1}\) 。
- 當 \(n>1\) 時，記 \(B(x)=A^{-1}(x)\) ，此時應該滿足
\[\begin{equation}A(x)B(x)\equiv 1\pmod {x^n}\end{equation}\]
假設在 \(\mod x^{\lceil \frac{n}{2} \rceil}\) 意義下 \(A(x)\) 的逆元是 \(B'(x)\) 並且我們已經求出，那么
\[\begin{equation}A(x)B'(x) \equiv 1 \pmod {x^{\lceil \frac{n}{2} \rceil}} \end{equation}\]
再將 \((2)\) 放在 \(\mod x^{\lceil \frac{n}{2} \rceil}\) 意義下
\[\begin{equation}A(x)B(x)\equiv 1\pmod {x^{\lceil \frac{n}{2} \rceil}}\end{equation}\]
由 \((3)-(4)\) 得到
\[B(x) - B'(x) \equiv 0 \pmod {x^{\lceil \frac{n}{2} \rceil}}\]
兩邊同時平方
\[B^2(x) - 2B'(x)B(x) + B'^2(x) \equiv 0 \pmod {x^n}\]
解釋一下平方后為什么模的 \(x^{\lceil \frac{n}{2} \rceil}\) 也會平方。
因為，左邊多項式在 \(\mod x^n\) 意義下為 \(0\) ，那么就說明其 \(0\) 到 \(n-1\) 次項系數都為 \(0\) ，平方了之后，對於第 \(0\leq i\leq 2n-1\) 項，其系數 \(a_i\) 為 \(\sum_{j=0}^i a_ja_{i-j}\) ，很明顯 \(j\) 和 \(i-j\) 之間必然有一個值小於 \(n\) ，因此 \(a_i\) 必然是 \(0\) ，也就是說平方后在 \(\mod x^{2n}\) 意義下仍然為 \(0\) 。
這時我們只要在等式兩邊同乘上 \(A(x)\) ，移項得
\[B(x) \equiv 2B'(x) - A(x)B'^2(x) \pmod {x^n}\]
這樣就可以得到 \(\mod x^n\) 意義下的逆元了，利用 \(\text{FFT}\) 加速之后可以做到在 \(O(n\log_2n)\) 時間內解決當前問題。
由主定理，最后總的時間復雜度也就是 \(O(n\log_2n)\) 的。
順便一提，由這個過程可以看出，一個多項式有沒有逆元完全取決於其常數項是否有逆元。
算法實現
#include <bits/stdc++.h>
using namespace std;
const int mod = 998244353;
const int N = 4e5;
int n, a[N+5], b[N+5], L, R[N+5], len, tmp[N+5];
int quick_pow(int a, int b) {
int ans = 1;
while (b) {
if (b&1) ans = 1ll*ans*a%mod;
b >>= 1, a = 1ll*a*a%mod;
}
return ans;
}
void NTT(int *A, int o) {
for (int i = 0; i < len; i++) if (i < R[i]) swap(A[i], A[R[i]]);
for (int i = 1; i < len; i <<= 1) {
int gn = quick_pow(3, (mod-1)/(i<<1)), x, y;
if (o == -1) gn = quick_pow(gn, mod-2);
for (int j = 0; j < len; j += (i<<1)) {
int g = 1;
for (int k = 0; k < i; k++, g = 1ll*g*gn%mod) {
x = A[j+k], y = 1ll*A[j+k+i]*g%mod;
A[j+k] = (x+y)%mod;
A[j+k+i] = (x-y+mod)%mod;
}
}
}
}
void poly_inv(int *A, int *B, int deg) {
// deg 表示多項式的度
if (deg == 1) {B[0] = quick_pow(A[0], mod-2); return; }
poly_inv(A, B, (deg+1)>>1);
for (L = 0, len = 1; len <= (deg<<1); len <<= 1) ++L;
// A*B 的度為 deg^2
for (int i = 0; i < len; i++) R[i] = (R[i>>1]>>1)|((i&1)<<(L-1));
for (int i = 0; i < deg; i++) tmp[i] = A[i];
for (int i = deg; i < len; i++) tmp[i] = 0;
for (int i = (deg+1)>>1; i < len; i++) B[i] = 0;
//注意將高次項系數補為 0
NTT(tmp, 1), NTT(B, 1);
for (int i = 0; i < len; i++)
B[i] = 1ll*B[i]*((2ll-1ll*tmp[i]*B[i]%mod+mod)%mod)%mod;
NTT(B, -1); int inv = quick_pow(len, mod-2);
for (int i = 0; i < len; i++) B[i] = 1ll*B[i]*inv%mod;
}
void work() {
scanf("%d", &n); for (int i = 0; i < n; i++) scanf("%d", &a[i]);
poly_inv(a, b, n); for (int i = 0; i < n; i++) printf("%d ", b[i]);
}
int main() {work(); return 0; }求第二類斯特林數
第二類斯特林數
定義
將 \(n\) 個有區別的球放入 \(m\) 個無區別的盒子中非空的方案數，記為 \(S(n,m)\) 或 \(\begin{Bmatrix}n\\m\end{Bmatrix}\) 。
遞推式
- \(\begin{Bmatrix}i\\0\end{Bmatrix}=0,i\in\mathbb{N_+}\) ， \(\begin{Bmatrix}0\\0\end{Bmatrix}=1\)
- \(\begin{Bmatrix}n\\m\end{Bmatrix}=\begin{Bmatrix}n-1\\m-1\end{Bmatrix}+m\begin{Bmatrix}n-1\\m\end{Bmatrix}\)
① 邊界情況顯然；
② 考慮第 \(n\) 個球如何放：放在之前的盒子里面，則共 \(m\) 方法；或新開一個盒子。
通項公式
\[S(n,m)=\frac{1}{m!} \sum _{k=0}^m (-1)^k{m\choose k}(m-k)^n\]
證明的話大致就是容斥原理， \(k\) 枚舉有多少個集合是空的，每種情況有 \(m\choose k\) 種空集情況，\(n\) 個元素可以放進非空的 \(m-k\) 個集合中。這樣求出來的答案是有序的，所以我們除以 \(m!\) 使得其變為無序。
\(\text{NTT}\) 優化
由第二類斯特林數的通項公式，我們將組合數拆開，得到
\[S(n,m)=\sum _{k=0}^m \frac{(-1)^k}{k!}\frac{(m-k)^n}{(m-k)!}\]
記多項式
\[\begin{aligned}C(x)&=\sum_{i=0}^\infty S(n,i)x^i\\A(x)&=\sum_{i=0}^\infty \frac{(-1)^i}{i!}x^i\\B(x)&=\sum_{i=0}^\infty\frac{i^n}{i!}x^i\end{aligned}\]
那么 \(C(x)=A(x)B(x)\) 。 \(\text{NTT}\) 優化即可。
快速沃爾什變換
我們回到多項式乘法的概念：我們已知多項式 \(A(x),B(x)\) 。要求多項式 \(C(x)=A(x)B(x)\) 。其中
\[c_i=\sum_{j+k=i}a_jb_k\]
似乎求和式下面的 j+k=i 比較單調。我們試着將 + 換成其他的符號。
考慮如何求
\[c_i=\sum_{j\oplus k=i}a_jb_k\]
其中 \(\oplus\) 為按位運算符號。包括常用的 xor and or ，即“按位異或”，“按位與”，“按位或”。
由於這三種卷積具有相似性，這里僅舉 xor 為例，其余兩種可以類比推出。
\(xor\) 卷積
注意下文中的 \(\oplus\) 表示“按位異或”。
我們類比 \(\text{FFT}\) 的過程。對於 \(C(x)=A(x)B(x)\) ，那么滿足
\[\text{DFT}(A(x))_i\times\text{DFT}(B(x))_i=\text{DFT}(C(x))_i\]
其中 \(\text{DFT}\) 表示多項式系數表示轉點值表示的過程。
我們考慮是否能也構造一個變換 \(\text{tf}\) ，使得 \(C(x)=A(x)\oplus B(x)\) 滿足
\[\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\]
由於是“按位異或”，我們不妨先考慮只含一位的情況，此時多項式只有兩項分別為 \(0\) 和 \(1\) 。記 \(A_0=0,A_1=1\) 。 \(B,C\) 類似。那么
\[\begin{aligned}\text{tf}(A)&=<A_0+A_1,A_0-A_1>\\\text{tf}(B)&=<B_0+B_1,B_0-B_1>\\\text{tf}(C)&=<C_0+C_1,C_0-C_1>\end{aligned}\]
至於為什么是這種形式，可以結合式子
\[\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\]
得出 \[\begin{cases}C_0&=A_0B_0+A_1B_1\\C_1&=A_0B_1+A_1B_0\end{cases}\]
顯然這是滿足異或卷積的表達式的，成立。
推廣到一般的情況，當 \(A\) 的長度為 \(2^k\) 時，我們記 \(A_0\) 為前 \(2^{k-1}\) 位， \(A_1\) 為后 \(2^{k-1}\) 位。用數學歸納法證明。容易發現 \(A_0\) 中每一項的最高位為 \(0\) ， \(A_1\) 中每一項的最高位為 \(1\) 。
只要證明滿足
\[\begin{aligned}\text{tf}(A)&=<\text{tf}(A_0)+\text{tf}(A_1),\text{tf}(A_0)-\text{tf}(A_1)>\\\text{tf}(B)&=<\text{tf}(B_0)+\text{tf}(B_1),\text{tf}(B_0)-\text{tf}(B_1)>\\\text{tf}(C)&=<\text{tf}(C_0)+\text{tf}(C_1),\text{tf}(C_0)-\text{tf}(C_1)>\end{aligned}\]
時， \(\text{tf}(A(x))_i\times\text{tf}(B(x))_i=\text{tf}(C(x))_i\) 依舊成立即可。
我們依舊暴力拆開，得到
\[\begin{cases}\text{tf}(C_0)&=\text{tf}(A_0)\text{tf}(B_0)+\text{tf}(A_1)\text{tf}(B_1)\\\text{tf}(C_1)&=\text{tf}(A_0)\text{tf}(B_1)+\text{tf}(A_1)\text{tf}(B_0)\end{cases}\]
由於異或每一位是獨立，而這里如果我們把 \(C\) 按照最高位為 \(0\) 或 \(1\) 分成兩部分，最高位的異或和其它位不相關。而 \(\text{tf}\) 已經將除最高位外的其他所有位處理好了。顯然是滿足條件的。注意： \(\text{tf}\times\text{tf}\) 是按位乘的，而不是卷積。
似乎我們已經做好了類似“系數轉點值”的過程；考慮逆過程 \(\text{utf}\) 。依舊用同樣的方法驗證，先考慮只含一位的情況。
\[\text{utf}(A)=\left<\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}\right>\]
由於 \(\text{uft}\) 是多項式，注意是“卷積”的形式。
得到 \(\begin{cases}C_0&=A_0B_0\\C_1&=A_1B_1\end{cases}\) 。顯然滿足。
接着就直接數歸來證即可。方法類似之前證明 \(\text{tf}\) 的過程。
結論（三種卷積求法）
正向 \(\text{tf}\)
- \(xor\) 卷積： \(\text{tf}(A)=<A_0+A_1,A_0-A_1>\)
- \(and\) 卷積： \(\text{tf}(A)=<A_0+A_1,A_1>\)
- \(or\) 卷積： \(\text{tf}(A)=<A_0,A_0+A_1>\)
逆向 \(\text{utf}\)
- \(xor\) 卷積： \(\text{utf}(A)=\left<\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}\right>\)
- \(and\) 卷積： \(\text{utf}(A)=<A_0-A_1,A_1>\)
- \(or\) 卷積： \(\text{utf}(A)=<A_0,A_1-A_0>\)
算法實現
這里僅提供 \(xor\) 卷積的模板，其它情況類似。
值得注意的是 \(\text{FWT}\) 與 \(\text{NTT},\text{FFT}\) 並不完全相似。只是用了類比的思想得到 \(\text{FWT}\) 的 \(\text{tf}\) 和 \(\text{utf}\) 過程，本質上是不同的。
所以說代碼實現也只是借用了迭代的思想等，一些操作如“交換初始位置”，“逆變換后除以 \(n\) ”，是不需要的。
void FWT(int *A, int o) {
for (int i = 1; i < n; i <<= 1)
for (int j = 0; j < n; j += (i<<1))
for (int k = 0; k < i; k++) {
int x = A[k+j], y = A[k+j+i];
A[k+j] = (x+y)%mod, A[k+j+i] = (x-y+mod)%mod;
if (o == -1)
A[k+j] = 1ll*A[k+j]*inv2%mod,
A[k+j+i] = 1ll*A[k+j+i]*inv2%mod;
}
}