[總結]多項式類數學相關（定理&證實&板子）

目錄

• 寫在前面
• 前置技能
  • 多項式相關
    • 多項式的係數表示
    • 多項式的點值表示
  • 複數相關
    • 複數的意義
    • 複數的基本運算
    • 單位根
  • 代碼相關
• 多項式乘法
• 快速傅里葉變換
  • DFT
  • IDFT
  • 算法實現
    • 遞歸實現
    • 迭代實現
• 快速數論變換
  • 原根
  • 算法實現
  • 模數任意的解決方案
• 應用
  • 快速卷積
  • 多項式求逆
    • 基本概念
    • 求解方法
    • 算法實現
  • 求第二類斯特林數
    • 第二類斯特林數
    • \(\text{NTT}\) 優化
• 快速沃爾什變換
  • \(xor\) 卷積
  • 結論（三種卷積求法）
    • 正向 \(\text{tf}\)
    • 逆向 \(\text{utf}\)
  • 算法實現

寫在前面

1. 因爲上一篇總結的版面限制，特開此文來記錄 \(\text{OI}\) 中多項式類數學相關的問題。
2. 該文啓發於Miskcoo的博客，甚至一些地方直接引用，在此特別說明；若文章中出現錯誤，煩請告知。
3. 感謝你的造訪。

前置技能

多項式相關

形同 \(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)\) 。

若是可以找到一種有效的方法幫助咱們在多項式的點值表示和係數表示之間轉換，咱們就能夠快速地計算多項式的乘法了，快速傅里葉變換就能夠作到這一點。

快速傅里葉變換

由剛纔的設想，咱們只需按序進行下述三次操做，便可下降多項式乘法的複雜度：

1. 將 \(A(x),B(x)\) 的係數表達轉爲點值表達；
2. 將 \(A(x),B(x)\) 的點值表達相乘；
3. 將獲得的結果轉爲係數表達。

咱們已知 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*}\]

而對於這個式子容易發現

1. 當 \(i=j\) 時， \(\omega_n^{k(j-i)}=\omega_n^{0}=1\) ，故 \(e_{ij}=\sum_{k=0}^{n-1}1=n\) 。
2. 當 \(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\) 的什麼性質。

1. \(\omega_n^n=1\) ，有這樣一個初值，方便後面的計算；
2. \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) 是互不相同的，這樣帶入計算出來的點值才能夠用來還原出係數；
3. \(\omega_n^2=\omega_{\frac{n}{2}}, \omega_n^{\frac{n}{2}+k}=-\omega_n^k\) ，這使得在按照指數奇偶分類以後可以把帶入的值也減半使得問題規模可以減半。
4. \[\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)\) 。

1. 當 \(n=1\) 時，容易發現 \(A(x)\equiv c\pmod x\) 此時 \(c\) 是一個常數，故 \(A^{-1}(x)\) 也是一個常數，即 \(c^{-1}\) 。
2. 當 \(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}\) 。

遞推式

1. \(\begin{Bmatrix}i\\0\end{Bmatrix}=0,i\in\mathbb{N_+}\) ， \(\begin{Bmatrix}0\\0\end{Bmatrix}=1\)
2. \(\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}\)

1. \(xor\) 卷積： \(\text{tf}(A)=<A_0+A_1,A_0-A_1>\)
2. \(and\) 卷積： \(\text{tf}(A)=<A_0+A_1,A_1>\)
3. \(or\) 卷積： \(\text{tf}(A)=<A_0,A_0+A_1>\)

逆向 \(\text{utf}\)

1. \(xor\) 卷積： \(\text{utf}(A)=\left<\frac{A_0+A_1}{2},\frac{A_0-A_1}{2}\right>\)
2. \(and\) 卷積： \(\text{utf}(A)=<A_0-A_1,A_1>\)
3. \(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;
            }
}