多項式

空白菌發表於2024-03-09

目錄
  • 多項式
    • 多項式基礎
      • 數域的定義
      • 多項式的定義與基本性質
      • 多項式帶餘式除法
      • 形式冪級數的定義
      • 冪級數的導數和不定積分
      • 常見冪級數展開
    • 多項式插值
      • 多項式插值的定義
      • 多項式插值的方法
        • 拉格朗日插值法
        • 重心拉格朗日插值法
    • 加法卷積
      • 加法卷積的定義
      • 加法卷積的變換
        • 快速傅立葉變換（FFT）
        • 快速數論變換（NTT）
        • 任意模數NTT（MTT）
    • 多項式初等函式
      • 多項式牛頓迭代
      • 多項式求逆
      • 多項式除法&取模
      • 多項式開方
      • 多項式對數函式
      • 多項式指數函式
      • 多項式冪函式
      • 多項式三角函式
      • 多項式反三角函式
    • 位運算卷積
      • 位運算卷積的定義
      • 快速沃爾什變換（FWT）
        • 異或卷積
        • 或卷積
        • 與卷積
    • 子集卷積
      • 子集卷積的定義
      • 快速莫比烏斯變換（FMT）
        • 並卷積
        • 交卷積
        • 子集卷積
    • 多項式分治
      • 多項式多點求值
      • 多項式快速插值
      • 分治加法卷積
    • 多項式線性齊次遞推
    • 多項式平移

多項式

多項式基礎

數域的定義

定義 複數集的子集 \(K\) ，滿足 \(0,1 \in K\) 且元素的和差積商（除數不為 \(0\) ）都屬於 \(K\) ，那麼稱 \(K\) 是一個數域。

關於群環域的詳細定義可以看看抽代，這裡只提及必須知道的。

例如有理數集、複數集、模素數 \(m\) 剩餘系都是數域，但整數集不是數域。

多項式的定義與基本性質

定義 設 \(a_0, a_1, \cdots ,a_n\) 是數域 \(K\) 上的數，其中 \(n\) 為非負整數，那麼稱 \(f(x) = \displaystyle \sum_{i = 0}^n a_i x^i\) 是數域 \(K\) 上的一元多項式，其中 \(a_ix^i\) 是 \(f(x)\) 的 \(i\) 次項，\(a_i\) 則是 \(f(x)\) 的 \(i\) 次項係數。

此外，也可以用 \([x^i]f(x)\) 表示多項式 \(f(x)\) 的 \(i\) 次項係數。

注意，這裡的 \(x\) 是形式上的記號，而非真正的數。換句話說，單獨寫出係數序列也能代表一個多項式， \(x\) 的次數更多時候只是用來區分系數。

一元多項式環的定義 數域 \(K\) 上多項式的全體，稱為一元多項式環，記作 \(K[x]\) ，而 \(K\) 稱為 \(K[x]\) 的係數域。

次數的定義 對於多項式 \(f(x)\) ，其係數非零的最高次項的次數稱為多項式的次數，記作 \(\partial(f(x))\) 。

相等的定義 若兩個多項式對應次數的各系數均相等，那麼這兩個多項式相等。

零多項式的定義 係數全為 \(0\) 的多項式，記為 \(0\) ，其次數不考慮。

多項式加法的定義 對於兩個多項式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ，則 \(\displaystyle f(x) + g(x) = \sum_{i=0}^{\max(n,m)} (a_i + b_i)x^i\) 。

多項式乘法的定義 對於兩個多項式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ，則 \(\displaystyle f(x) \cdot g(x) = \sum_{i=0}^n \sum_{j=0}^m a_ib_jx^{i+j}\) 。

多項式乘法有一個等價形式，\(\displaystyle f(x) \cdot g(x) = \sum_{k=0}^{n+m} x^{k} \sum_{i = 0}^k a_ib_{k-i}\) ，即求兩個多項式係數的加法卷積（下標之和相等的位置的值的乘積之和），這是今後許多問題可以轉化為多項式計算的關鍵。

多項式複合的定義 對於兩個多項式 \(\displaystyle f(x) = \sum_{i = 0}^n a_i x^i, g(x) = \sum_{i = 0}^m b_i x^i\) ，則 \(\displaystyle f(g(x)) = a_0 + \sum_{i=1}^n a_ig^i(x)\) 。

性質1 數域 \(K\) 上的兩個多項式經過加、減、乘等運算後，所得結果仍然是數域 \(K\) 上的多項式。

性質2 對於兩個多項式 \(f(x),g(x)\) ，滿足加法結合律、加法交換律、乘法結合律、乘法交換律、乘法對加法分配律、乘法消去律。

性質3 任意 \(n+1\) 個不同的取樣點，就可以唯一確定一個 \(n\) 次多項式。

多項式帶餘式除法

定理1（帶餘式除法） 在一元多項式環 \(K[x]\) 中，任意兩個多項式 \(A(x),B(x)\) 且 \(B(x) \neq 0\) ，一定存在唯一的兩個多項式 \(Q(x),R(x)\) 滿足 \(A(x) = Q(x)B(x) + R(x)\) ，其中 \(\partial(R(x)) < \partial(B(x))\) 或 \(R(x) = 0\) ，稱 \(Q(x)\) 為 \(A(x)\) 除以 \(B(x)\) 的商， \(R(x)\) 為 \(A(x)\) 除以 \(B(x)\) 的餘式。

大部分數論整除同餘的性質都可以類似地應用到多項式上，後面就不展開了。

形式冪級數的定義

定義 設 \(a_0, a_1, \cdots ,a_n\) 是數域 \(K\) 上的數，那麼稱 \(f(x) = \displaystyle \sum_{i = 0}^{\infin} a_i x^i\) 是數域 \(K\) 上的形式冪級數，簡稱冪級數。

形式冪級數環的定義 數域 \(K\) 上形式冪級數的全體，稱為形式冪級數環，記作 \(K[[x]]\) 。

冪級數可以看作是一元多項式的擴充套件，其具有更多良好的性質，如形式導數和形式不定積分等。

在模意義下，冪級數可等價為有限項的多項式，因此通常會把多項式擴充套件到冪級數上，藉由冪級數的諸多性質得到許多有用的結論，例如模意義下多項式的初等函式。

冪級數的導數和不定積分

注意，極限在環上可能並不存在，但依然可以在形式上的定義導數與積分。

形式導數 設形式冪級數 \(\displaystyle f(x) = \sum_{i = 0}^{\infin}a_ix^i\) ，其形式導數 \(\displaystyle f'(x) = \sum_{i = 1}^{\infin}ia_ix^{i-1}\) 。

此外，我們也可將 \(f(x)\) 的 \(t\) 階導數記作 \(f^{(t)}(x)\) 。

形式不定積分 設形式冪級數 \(\displaystyle f(x) = \sum_{i = 0}^{\infin}a_ix^i\) ，其形式不定積分 \(\displaystyle \int f(x) \text{d} x = \sum_{i = 1}^{\infin}ia_ix^{i-1} + C\) 。

其他的基本求導法則皆可適用，就不再展開了。

常見冪級數展開

\[\begin{aligned} e^x &= \sum_{i = 0}^{\infin} \frac{1}{i!}x^i \\ \sin x &= \sum_{i = 0}^{\infin} \frac{(-1)^{i}}{(2i+1)!}x^{2i+1} \\ \cos x &= \sum_{i = 0}^{\infin} \frac{ (-1)^{i}}{(2i)!}x^{2i} \\ \frac{1}{1+x} &= \sum_{i = 0}^{\infin} (-1)^ix^i \\ (1+x)^{\alpha} &= \sum_{i = 0}^{\infin} \frac{\alpha^{\underline i}}{i!}x^i \\ \ln(1+x) &= \sum_{i = 1}^{\infin} \frac{(-1)^{i-1}}{i}x^i \\ \arctan x &= \sum_{i = 0}^{\infin} \frac{(-1)^{i}}{2i+1}x^{2i+1} \\ \end{aligned} \]

多項式插值

多項式插值的定義

定義 給定 \(n\) 個點 \((x_1,y_1), \cdots, (x_n, y_n)\) ，其中 \(x_i\) 互不相同，求這些點確定的 \(n-1\) 次多項式函式 \(f(x)\) 的過程，稱為多項式插值。

多項式插值的方法

拉格朗日插值法

考慮構造 \(n\) 個輔助函式 \(\displaystyle g_i(x) = y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j}\) ，顯然 \(g_i(x)\) 滿足 \(g_i(x_i) = y_i\) 且 \(g_i(x_j) = 0, j \neq i\) 。

因此我們令 \(\displaystyle f(x) = \sum_{i = 1}^n g_i(x) = \sum_{i = 1}^n y_i \prod_{j \neq i} \frac{x-x_j}{x_i-x_j}\) 即可唯一確定所求多項式，此為拉格朗日插值公式。

其中，若 \(x_i = i\) ，可以預處理階乘以及 \(x - x_j\) 的前字尾積，將公式化簡為 \(O(n)\) 。

若我們只需要求出在 \(x = k\) 處的值，那麼代入即可。

若要求出具體的係數則設計多項式運算的模擬。

這裡只給出求單點值的程式碼。

時間複雜度 \(O(n^2)\)

空間複雜度 \(O(n)\)

int lagrange_poly(const vector<pair<int, int>> &point, int x) {
    int n = point.size() - 1;
    int ans = 0;
    for (int i = 1;i <= n;i++) {
        int res1 = point[i].second;
        int res2 = 1;
        for (int j = 1;j <= n;j++) {
            if (i == j) continue;
            res1 = 1LL * res1 * (x - point[j].first + P) % P;
            res2 = 1LL * res2 * (point[i].first - point[j].first + P) % P;
        }
        (ans += 1LL * res1 * qpow(res2, P - 2) % P) %= P;
    }
    return ans;
}

重心拉格朗日插值法

若插值點會新增 \(q\) 次，每次重新計算 \(f(k)\) 都是 \(O(n^2)\) ，我們需要對公式做些調整。

\(\displaystyle f(x) = \sum_{i = 1}^n y_i \prod_{j \neq i} \frac{x-x_j}{x_i-x_j} = \prod_{i=1}^n (x - x_i) \sum_{i = 1}^n \frac{y_i}{(x-x_i)\prod_{j \neq i} (x_i-x_j)}\) 。

我們設 \(\displaystyle A = \prod_{i=1}^n (x - x_i) , B(i)=\displaystyle \prod_{j \neq i} (x_i-x_j)\) 。

每次新增插值點時 \(O(1)\) 更新 \(A\) ， \(O(n)\) 更新 \(B(i)\) 後，即可在 \(O(n)\) 內得到新的 \(\displaystyle f(k) = A \sum_{i = 1}^n \frac{y_i}{(k-x_i)B(i)}\) 。

程式碼可借鑑的拉格朗日插值法做修改，這裡就不給出了。

時間複雜度 \(O(n^2 + qn)\)

空間複雜度 \(O(n)\)

加法卷積

加法卷積的定義

定義 對於兩個序列 \(f,g\) ，它們的加法卷積序列 \(h\) 滿足 \(\displaystyle h_k = \sum_{i + j = k} f_ig_j\) 。

把序列當作多項式係數理解， \(h\) 其實就是 \(f,g\) 表示的多項式的乘積的係數，因此可以用加法卷積的演算法加速多項式乘積，下面也會用多項式的角度介紹加法卷積的演算法。

加法卷積的變換

快速傅立葉變換（FFT）

多項式在係數域直接加法卷積是 \(O(n^2)\) 的，但如果我們在若干個不同位置取兩個多項式的點值（取樣點），容易發現這些點值相乘後得到的新點值就落在所求的多項式上，最後只要把點值變換回係數，就得到目標多項式。

換句話說，係數域的加法卷積可以變換為點值域的對應相乘，而對應相乘這個過程是 \(O(n)\) 的，現在我們只需要一個能夠快速在係數域和點值域之間變換演算法即可。

這也是大多數變換加速卷積的核心思路，即找到一個快速的可逆變換，使得兩個序列變換後的對應位置做運算的結果，恰好為兩個序列卷積的變換，最後逆變換回去就是兩個序列的卷積。

接下來就是加法卷積的需要的變換，離散傅立葉變換。

離散傅立葉變換（DFT）

首先，點值域的點不能隨便取，我們要取 \(n\) 次單位根 \(\omega_n\) 的 \(n\) 個冪 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n-1}\) ，\(n\) 要大於等於多項式的項數。為了方便，我們通常需要令 \(n\) 為 \(2\) 的冪。

\(n\) 次單位根等價於將複平面單位圓弧劃分為 \(n\) 等分，其中第 \(k\) 份即 \(\omega_n^k = \cos \dfrac{2k\pi}{n} + \text{i}\sin \dfrac{2k\pi}{n}\) 。

單位根具有許多有用的性質：

  1. 互異性：若 \(i \neq j\) ，則 \(\omega_n^i \neq \omega_n^j\) 。
  2. 折半律：\(\omega_{2n}^{2i} = \omega_{n}^{i}\) 。
  3. 週期律：\(\omega_n^{i + n} = \omega_n^i\) 。
  4. 半週期律： \(\omega_n^{i + \frac{n}{2}} = -\omega_n^i\) 。

互異性保證了 \(n\) 個點一定互不相同，接下來考慮如何求值。

設 \(\displaystyle f(x) = \sum_{i = 0}^{n-1} a_i x^i\) ，那麼顯然有 \(f(x) = a_0 + a_2x^2 + \cdots + a_{n-1}x^{n-1} + x(a_1 + a_3x^2 + \cdots + a_nx^n)\) 。

設 \(f_1(x) = a_0 + a_2x + \cdots a_{n-1}x^{n-1} ,f_2(x) = a_1 + a_3x + \cdots + a_nx^n\) ，那麼有 \(f(x) = f_1(x^2) + xf_2(x^2)\) 。

當 \(i \in [0, \dfrac{n}{2} - 1]\) 時，我們代入單位根 \(\omega_n^i\) 可得

\[\begin{aligned} f(\omega_n^i) &= f_1((\omega_n^i)^2) +\omega_n^if_2((\omega_n^i)^2) \\ &= f_1(\omega_n^{2i}) +\omega_n^if_2(\omega_n^{2i}) \\ &= f_1(\omega_{\frac{n}{2}}^{i}) +\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

我們代入單位根 \(\omega_n^{i + \frac{n}{2}}\) 可得

\[\begin{aligned} f(\omega_n^{i + \frac{n}{2}}) &= f_1((\omega_n^{i + \frac{n}{2}})^2) +\omega_n^{i + \frac{n}{2}}f_2((\omega_n^{i + \frac{n}{2}})^2) \\ &= f_1(\omega_n^{2i}) -\omega_n^if_2(\omega_n^{2i}) \\ &= f_1(\omega_{\frac{n}{2}}^{i}) -\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

注意到 \(f_1(\omega_{\frac{n}{2}}^{i})\) 和 \(f_2(\omega_{\frac{n}{2}}^{i})\) 正是子問題的答案。

因此一個大小為 \(n\) 的問題，變成兩個大小為 \(\dfrac{n}{2}\) 的子問題外加 \(O(n)\) 複雜度計算，遞迴下去總體複雜度是 \(O(n\log n)\) 的。（如果隨便取點，問題規模不會折半，也就不能快速了）

於是，多項式係數域到點值域的快速正變換就找到了。

位逆序置換

遞迴的常數較大，我們希望改為迭代，考慮 \(2^3\) 項多項式的變換過程：

  1. 初始層：\(\{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\) 。
  2. 第一層：\(\{a_0, a_2, a_4, a_6\},\{a_1, a_3, a_5, a_7\}\) 。
  3. 第二層：\(\{a_0, a_4\},\{a_2, a_6\},\{a_1, a_5\}, \{a_3, a_7\}\) 。
  4. 第三層：\(\{a_0\}, \{a_4\},\{a_2\}, \{a_6\},\{a_1\}, \{a_5\}, \{a_3\}, \{a_7\}\) 。

我們需要從下往上迭代，那麼就需要知道最後一層的順序。

容易知道，\(a_i\) 最後會出現在 \(a_{rev_i}\) ，其中 \(rev_i\) 表示 \(i\) 的二進位制逆序，例如 \(110\) 的逆序就是 \(011\) 。

根據 \(rev\) 的定義，我們可以遞推它在 \(n\) 項多項式的情況：

\[\begin{aligned} rev_i = \dfrac{rev_{\frac{i}{2}}}{2} + [2 \not\mid i] \cdot 2^{\log n - 1} \end{aligned} \]

因此，我們預處理 \(rev\) 後，將對應係數置換即可迭代。

蝶形運算最佳化

在上面，當我們求出 \(f_1(x)\) 和 \(f_2(x)\) 的各 \(\dfrac{n}{2}\) 個點值後，設 \(i \in [0, \dfrac{n}{2}-1]\) ，那麼

\[\begin{aligned} f(\omega_n^i) &= f_1(\omega_{\frac{n}{2}}^{i}) +\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ f(\omega_n^{i + \frac{n}{2}}) &= f_1(\omega_{\frac{n}{2}}^{i}) -\omega_n^if_2(\omega_{\frac{n}{2}}^{i}) \\ \end{aligned} \]

注意到 \(f(\omega_n^i)\) 和 $ f(\omega_n^{i + \frac{n}{2}})$ 分別在 \(i\) 和 \(i + \frac{n}{2}\) 位置上，而 \(f_1(\omega_{\frac{n}{2}}^{i})\) 和 \(f_2(\omega_{\frac{n}{2}}^{i})\) 也恰好在 \(i\) 和 \(i + \frac{n}{2}\) 位置上，因此我們不需要額外空間，只需要原地覆蓋即可。

離散傅立葉逆變換（IDFT）

現在，我們嘗試推導一下逆變換。

我們定義點值多項式 \(\displaystyle F(x) = \sum_{i = 0}^{n-1} f(\omega_n^{i})x^i\) ，即 \(f(x)\) 點值當作係數的多項式。

我們代入 \(x = \omega_n^k\) ，那麼 \(\displaystyle F(\omega_n^k) = \sum_{i = 0}^{n-1} \omega_n^{ik}\sum_{j = 0}^{n-1} a_j\omega_n^{ij} = \sum_{j = 0}^{n-1} a_j\sum_{i = 0}^{n-1} \omega_n^{i(k+j)}\) 。

構造輔助多項式 \(\displaystyle G(x) = \sum_{i = 0}^{n-1} x^i\) ，那麼 \(\displaystyle F(\omega_n^k) = \sum_{j=0}^{n-1}a_jG(\omega_n^{j+k})\) 。

考慮 \(G(\omega_n^i)\) ，當 \(i = 0\) 時 \(G(\omega_n^i) = n\) ，否則單位根兩兩配對 \(G(\omega_n^i) = 0\) 。

因此 \(F(\omega_n^k) = na_{n-k}\) ，特別地當 \(k = 0\) 時 \(F(\omega_n^0) = na_0\) ，所以我們只需要對點值多項式進行一次DFT，隨後將 \([1,n-1]\) 項反轉，最後對所有項除以 \(n\) 即可還原多項式。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

const db PI = acos(-1.0);

vector<int> rev;
vector<complex<db>> Wn = { {0, 0}, {1, 0} }; // 0, w20, w40, w41, w80, w81, w82, w83, ...
void FFT(vector<complex<db>> &A, bool is_inv) {
    int n = A.size();

    if (rev.size() != n) {
        int k = __builtin_ctz(n) - 1;
        rev.resize(n);
        for (int i = 0;i < n;i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
    }
    for (int i = 0;i < n;i++) if (i < rev[i]) swap(A[i], A[rev[i]]);

    if (Wn.size() < n) {
        int k = Wn.size();
        Wn.resize(n);
        while (k < n) {
            complex<db> w(cos(PI / k), sin(PI / k));
            for (int i = k >> 1;i < k;i++) {
                Wn[i << 1] = Wn[i];
                Wn[i << 1 | 1] = Wn[i] * w;
            }
            k <<= 1;
        }
    }

    for (int k = 1;k < n; k <<= 1) {
        for (int i = 0;i < n;i += k << 1) {
            for (int j = 0;j < k;j++) {
                complex<db> u = A[i + j];
                complex<db> v = A[i + k + j] * Wn[k + j];
                A[i + j] = u + v;
                A[i + k + j] = u - v;
            }
        }
    }

    if (is_inv) {
        reverse(A.begin() + 1, A.end());
        for (int i = 0;i < n;i++) A[i] /= n;
    }
}

vector<complex<db>> poly_mul(vector<complex<db>> A, vector<complex<db>> B) {
    if (A.empty() || B.empty()) return vector<complex<db>>();
    int n = 1, sz = A.size() + B.size() - 1;
    while (n < sz) n <<= 1;
    A.resize(n);
    B.resize(n);
    FFT(A, 0);
    FFT(B, 0);
    for (int i = 0;i < n;i++) A[i] *= B[i];
    FFT(A, 1);
    A.resize(sz);
    return A;
}

快速數論變換（NTT）

考慮在素域內的多項式變換，我們發現原根剛好能代替單位根。

考慮一個素數 \(P\) ，表達為 \(P = r \cdot 2^k + 1\) ，其原根 \(G\) 的階為 \(\varphi(P) = P-1 = r \cdot 2^k\) ，當多項式含有 \(n = 2^{k'}\) 項時，我們令 \(G_n = G^{\frac{P-1}{n}}\) 那麼 \(G_n\) 等價為 \(n\) 次單位根。

注意到，一個素數能支援的多項式長度為 \(2^k\) ，因此 \(k\) 越大越好，不過常見的 \(10^9 + 7\) 就比較雞肋，因為 \(k = 1\) 。

NTT其他部分和FTT完全一致。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

const int P = 998244353, G = 3;

int qpow(int a, ll k) {
    int ans = 1;
    while (k) {
        if (k & 1) ans = 1LL * ans * a % P;
        k >>= 1;
        a = 1LL * a * a % P;
    }
    return ans;
}

std::vector<int> rev;
std::vector<int> Wn = { 0,1 }; // 0, w20, w40, w41, w80, w81, w82, w83, ...
/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
    explicit Poly(int n = 0, int val = 0) : std::vector<int>(n, val) {}
    explicit Poly(const std::vector<int> &src) : std::vector<int>(src) {}
    explicit Poly(const std::initializer_list<int> &src) : std::vector<int>(src) {}

    Poly modx(int k) const {
        assert(k >= 0);
        if (k > size()) {
            Poly X = *this;
            X.resize(k);
            return X;
        }
        else return Poly(std::vector<int>(begin(), begin() + k));
    }
    Poly mulx(int k) const {
        assert(k >= 0 || -k < size());
        Poly X = *this;
        if (k >= 0) X.insert(X.begin(), k, 0);
        else X.erase(X.begin(), X.begin() + (-k));
        return X;

    }
    Poly derive(int k = 0) const {
        if (k == 0) k = std::max((int)size() - 1, 0);
        Poly X(k);
        for (int i = 1;i < std::min((int)size(), k + 1);i++) X[i - 1] = 1LL * i * (*this)[i] % P;
        return X;
    }
    Poly integral(int k = 0) const {
        if (k == 0) k = size() + 1;
        Poly X(k);
        for (int i = 0;i < std::min((int)size(), k - 1);i++)  X[i + 1] = 1LL * qpow(i + 1, P - 2) * (*this)[i] % P;
        return X;
    }
    
    Poly &operator+=(const Poly &X) {
        resize(std::max(size(), X.size()));
        for (int i = 0;i < size();i++) ((*this)[i] += X[i]) %= P;
        return *this;
    }
    Poly &operator-=(const Poly &X) {
        resize(std::max(size(), X.size()));
        for (int i = 0;i < size();i++) ((*this)[i] -= X[i] - P) %= P;
        return *this;
    }
    Poly &operator*=(int x) {
        for (int i = 0;i < size();i++) (*this)[i] = 1LL * (*this)[i] * x % P;
        return *this;
    }
    Poly &operator/=(int x) {
        int val = qpow(x, P - 2);
        for (int i = 0;i < size();i++) (*this)[i] = 1LL * (*this)[i] * val % P;
        return *this;
    }
    friend Poly operator+(Poly A, const Poly &B) { return A += B; }
    friend Poly operator-(Poly A, const Poly &B) { return A -= B; }
    friend Poly operator*(Poly A, int x) { return A *= x; }
    friend Poly operator*(int x, Poly A) { return A *= x; }
    friend Poly operator/(Poly A, int x) { return A /= x; }
    friend Poly operator-(Poly A) { return (P - 1) * A; }
    friend std::ostream &operator<<(std::ostream &os, const Poly &X) {
        for (int i = 0;i < X.size();i++) os << X[i] << ' ';
        return os;
    }

    static void NTT(Poly &X, bool is_inv) {
        int n = X.size();

        if (rev.size() != n) {
            int k = __builtin_ctz(n) - 1;
            rev.resize(n);
            for (int i = 0;i < n;i++) rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
        }
        for (int i = 0;i < n;i++) if (i < rev[i]) std::swap(X[i], X[rev[i]]);

        if (Wn.size() < n) {
            int k = __builtin_ctz(Wn.size());
            Wn.resize(n);
            while (1 << k < n) {
                int w = qpow(G, P - 1 >> k + 1);
                for (int i = 1 << k - 1;i < 1 << k;i++) {
                    Wn[i << 1] = Wn[i];
                    Wn[i << 1 | 1] = 1LL * Wn[i] * w % P;
                }
                k++;
            }
        }

        for (int k = 1;k < n; k <<= 1) {
            for (int i = 0;i < n;i += k << 1) {
                for (int j = 0;j < k;j++) {
                    int u = X[i + j];
                    int v = 1LL * X[i + k + j] * Wn[k + j] % P;
                    X[i + j] = (u + v) % P;
                    X[i + k + j] = (u - v + P) % P;
                }
            }
        }

        if (is_inv) {
            std::reverse(X.begin() + 1, X.end());
            int inv = qpow(n, P - 2);
            for (int i = 0;i < n;i++) X[i] = 1LL * X[i] * inv % P;
        }
    }
    Poly &operator*=(Poly X) {
        if (empty() || X.empty()) return *this = Poly();
        int n = 1, sz = size() + X.size() - 1;
        while (n < sz) n <<= 1;
        resize(n);
        X.resize(n);
        NTT(*this, 0);
        NTT(X, 0);
        for (int i = 0;i < n;i++) (*this)[i] = 1LL * (*this)[i] * X[i] % P;
        NTT(*this, 1);
        resize(sz);
        return *this;
    }
    friend Poly operator*(Poly A, const Poly &B) { return A *= B; }
};

常用NTT模數：

\(r\cdot 2^k + 1\) \(r\) \(k\) \(g\)
5767169 11 19 3
7340033 7 20 3
23068673 11 21 3
104857601 25 22 3
167772161 5 25 3
469762049 7 26 3
998244353 119 23 3
1004535809 479 21 3
2013265921 15 27 31
2281701377 17 27 3
3221225473 3 30 5
75161927681 35 31 3
77309411329 9 33 7
206158430209 3 36 22
2061584302081 15 37 7
2748779069441 5 39 3
6597069766657 3 41 5
39582418599937 9 42 5
79164837199873 9 43 5
263882790666241 15 44 7
1231453023109121 35 45 3
1337006139375617 19 46 3
3799912185593857 27 47 5
4222124650659841 15 48 19
7881299347898369 7 50 6
31525197391593473 7 52 3
180143985094819841 5 55 6
1945555039024054273 27 56 5
4179340454199820289 29 57 3

任意模數NTT（MTT）

暫時不學。

多項式初等函式

初等函式 公式 方法 備註
乘法 \(f(x) \cdot g(x)\) NTT/FTT
求逆 \(f^{-1}(x) \equiv f^{-1}_0(x)(2 - f^{-1}_0(x)f(x)) \pmod{x^n}\) 牛頓迭代、乘法 常數項逆元存在
整除 \(\left[\dfrac{f(x)}{g(x)} \right]_R \equiv f_R(x)g_R^{-1}(x) \pmod{x^{n-m+1}}\) 求逆 除式非零
取模 \(f(x) \bmod g(x) = f(x) - g(x)\left[\dfrac{f(x)}{g(x)} \right]\) 整除 除式非零
開方 \(\sqrt{f(x)} \equiv \dfrac{1}{2} \left(\left( \sqrt{f(x)} \right)_0 + f(x)\left( \sqrt{f(x)} \right)_0^{-1} \right) \pmod{x^n}\) 牛頓迭代、求逆 首非零項是偶次項，且二次剩餘存在
對數函式 \(\displaystyle \ln f(x) \equiv \int f'(x)f^{-1}(x) \text{d}x \pmod{x^n}\) 求逆 常數項為 \(1\)
指數函式 \(\text{e}^{f(x)} \equiv \left(\text{e}^{f(x)}\right)_0 \left(1-\ln \left(\text{e}^{f(x)}\right)_0 + f(x) \right) \pmod{x^n}\) 牛頓迭代、對數函式 常數項為 \(0\)
冪函式 \(f^k(x) \equiv e^{k \ln f(x)} \pmod{x^n}\) 指數函式
三角函式 尤拉公式 指數函式 常數項為 \(0\)
反三角函式 求導積分 開方 常數項為 \(0\)

多項式牛頓迭代

給定多項式 \(g(x)\) ，求 \(f(x)\) ，滿足 \(g(f(x)) \equiv 0 \pmod{x^n}\) 。

考慮倍增法。

當 \(n = 1\) 時， \([x^0]g(f(x)) = 0\) 需要單獨解出。

假設在模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 時的 \(f(x)\) 的解是 \(f_0(x)\) ，那麼我們對 \(g(f(x))\) 在 \(f_0(x)\) 處泰勒展開有

\[\displaystyle \sum_{i=0}^{\infin} \dfrac{g^{(i)}(f_0(x))}{i!}(f(x) - f_0(x))^i \equiv 0 \pmod {x^n} \]

顯然，當 \(i \geq 2\) 時， \((f(x) - f_0(x))^i \equiv 0 \pmod{x^n}\) ，因此有

\[\displaystyle \sum_{i=0}^{\infin} \dfrac{g^{(i)}(f_0(x))}{i!}(f(x) - f_0(x))^i \equiv g(f_0(x)) + g'(f_0(x))(f(x) - f_0(x))) \equiv 0 \pmod {x^n} \]

最後化簡得

\[f(x) \equiv f_0(x) - \frac{g(f_0(x))}{g'(f_0(x))} \pmod{x^n} \]

這就是模意義下的牛頓迭代，每次都能倍增多項式有效係數，一些關鍵多項式初等函式需要由此推導。

模 \(x^n\) 是因為精確解實際上大機率是冪級數，但大部分時候我們的操作只需要前幾項，就能保證覆蓋所有有意義的部分，因此冪級數是不必要的，求出模意義下的就夠了。

多項式求逆

給定多項式 \(f(x)\) ，求 \(f^{-1}(x)\) ，滿足 \(f(x)f^{-1}(x) \equiv 1 \pmod{x^n}\) 。

設 \(g(f^{-1}(x)) = \dfrac{1}{f^{-1}(x)} - f(x) \equiv 0 \pmod{x^n}\) 。

當 \(n = 1\) 時，\([x^0]f^{-1}(x) = ([x^0]f(x))^{-1}\) ，因此需要保證常數項逆元存在。

假設模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解為 \(f_0(x)\) 。

根據牛頓迭代

\[f^{-1}(x) \equiv f^{-1}_0(x) - \dfrac{g(f^{-1}_0(x))}{g'(f^{-1}_0(x))} \equiv f^{-1}_0(x) - \dfrac{\dfrac{1}{f^{-1}_0(x)} - f(x)}{-\dfrac{1}{f^{-2}_0(x)}} \equiv f^{-1}_0(x)(2 - f^{-1}_0(x)f(x)) \pmod{x^n} \]

因此，我們可以用這個公式迭代出多項式的逆，複雜度由主定理 \(T(n) = T(n/2) + O(n \log n) = O(n \log n)\) 。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly inv(int n = 0) const {
        assert(size() && (*this)[0] > 0); // assert [x^0]f(x) inverse exists
        if (n == 0) n = size();
        Poly X{ qpow((*this)[0], P - 2) };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (X * (Poly{ 2 } - X * modx(k))).modx(k);
        }
        return X.modx(n);
    }
};

多項式除法&取模

給定多項式 \(f(x),g(x)\) ，求 \(q(x),r(x)\) ，滿足 \(f(x) = q(x)g(x) + r(x)\) 。

其中 \(\partial(q(x)) = \partial(f(x)) - \partial(g(x)), \partial(r(x)) < \partial(g(x))\) 。

設 \(n = \partial(f(x)), m = \partial(g(x))\) ，不妨設 \(\partial(r(x)) = m-1\)。

因為存在餘式，我們不能直接使用模 \(x^m\) 的多項式求逆。

設 \(f_R(x) = x^nf\left( \dfrac{1}{x} \right)\) ，其實就是將係數反轉。

我們對原式變形

\[d\begin{aligned} f(x) &= q(x)g(x) + r(x)\\ f\left( \dfrac{1}{x} \right) &= q\left( \dfrac{1}{x} \right)g\left( \dfrac{1}{x} \right) + r\left( \dfrac{1}{x} \right) \\ x^nf\left( \dfrac{1}{x} \right) &= x^nq\left( \dfrac{1}{x} \right)g\left( \dfrac{1}{x} \right) + x^nr\left( \dfrac{1}{x} \right) \\ f_R(x) &= g_R(x)q_R(x) + x^{n - m + 1} r_R(x) \end{aligned} \]

有 \(\partial(q_R(x)) = \partial(q(x)) = n-m < n-m+1\) ，因此在模 \(x^{n-m+1}\) 下 \(q_R(x)\) 是不會被影響的，而 \(x^{n-m+1}r_R(x)\) 會被模掉。

所以有 \(f_R(x) \cdot g^{-1}_R(x) \equiv q_R(x) \pmod{x^{n-m+1}}\) 。

求出 \(q_R(x)\) 後，反轉系數就是 \(q(x)\) ，最後 \(r(x) = f(x) - q(x)g(x)\) 。

實現上注意處理除式的後導 \(0\) ，會導致結果出錯，雖然一般題目不需要這個處理。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
    Poly &operator/=(Poly X) {
        while (X.size() && X.back() == 0) X.pop_back();
        assert(X.size()); // assert X(x) is not 0-polynomial
        if (size() < X.size()) return *this = Poly();
        std::reverse(begin(), end());
        std::reverse(X.begin(), X.end());
        *this = (modx(size() - X.size() + 1) * X.inv(size() - X.size() + 1)).modx(size() - X.size() + 1);
        std::reverse(begin(), end());
        return *this;
    }
    Poly &operator%=(Poly X) {
        while (X.size() && X.back() == 0) X.pop_back();
        return *this = (*this - *this / X * X).modx(X.size() - 1);
    }
    friend Poly operator/(Poly A, const Poly &B) { return A /= B; }
    friend Poly operator%(Poly A, const Poly &B) { return A %= B; }
};

多項式開方

給定多項式 \(f(x)\) ，求 \(\sqrt{f(x)} \bmod x^n\) 。

設 \(g(\sqrt{f(x)}) = \left( \sqrt{f(x)} \right)^2 - f(x) \equiv 0 \pmod {x^n}\) 。

當 \(n = 1\) 時， \([x^0]\sqrt{f(x)} = \sqrt{[x^0]f(x)}\) ，因此需要保證常數項二次剩餘存在。

假設模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解為 \(f_0(x)\) 。

根據牛頓迭代

\[\sqrt{f(x)} \equiv \left( \sqrt{f(x)} \right)_0 - \frac{\left( \sqrt{f(x)} \right)_0^2 - f(x)}{2\left( \sqrt{f(x)} \right)_0} \equiv \dfrac{1}{2} \left(\left( \sqrt{f(x)} \right)_0 + f(x)\left( \sqrt{f(x)} \right)_0^{-1} \right) \pmod{x^n} \]

程式碼沒實現求二次剩餘，目前只能對常數項為 \(1\) 的開方。

特別地，出現前導 \(0\) 時，前導 \(0\) 個數 \(cnt\) 是偶數（即第一個非零項是偶次）則多項式可以整體除以 \(x^{cnt}\) 再開方，最後乘 \(x^{cnt/2}\) ，否則無解。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

struct Poly : public std::vector<int> {
    Poly sqrt(int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        assert(!(cnt & 1) && (*this)[cnt] == 1); // assert cnt is even and [x^cnt]f(x) exists 2-residue 
        Poly X{ 1 };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (P + 1 >> 1) * (X + mulx(-cnt).modx(k) * X.inv(k)).modx(k);
        }
        return X.mulx(cnt >> 1).modx(n);
    }
};

多項式對數函式

給定多項式 \(f(x)\) ，求 \(\ln f(x) \bmod x^n\) 。

求導再積分後， \(\displaystyle \ln f(x) \equiv \int \frac{f'(x)}{f(x)} \text{d}x \pmod {x^n}\) 。

注意根據 \(\ln\) 的定義， \(f(x)\) 的常數項必須為 \(1\) ，否則對數無意義。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly log(int n = 0) const {
        assert(size() && (*this)[0] == 1); // assert [x^0]f(x) = 1
        if (n == 0) n = size();
        return (derive(n) * inv(n)).integral(n);
    }
};

多項式指數函式

給定多項式 \(f(x)\) ，求 \(e^{f(x)} \bmod x^n\) 。

設 \(g(e^{f(x)}) = \ln e^{f(x)} - f(x) \equiv 0 \pmod{x^n}\) 。

當 \(n = 1\) 時，\([x^0]e^{f(x)} = e^{[x^0]f(x)}\) ，因此需要保證常數項為 \(0\) ，否則無意義。

假設模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 的解為 \(f_0(x)\) 。

根據牛頓迭代

\[e^{f(x)} \equiv \left( e^{f(x)} \right)_0 - \frac{\ln \left( e^{f(x)} \right)_0 - f(x)}{\frac{1}{\left( e^{f(x)} \right)_0}} \equiv \left( e^{f(x)} \right)_0 \left( 1 - \ln \left( e^{f(x)} \right)_0 + f(x) \right) \pmod{x^n} \]

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly exp(int n = 0) const {
        assert(empty() || (*this)[0] == 0); // assert [x^0]f(x) = 0
        if (n == 0) n = size();
        Poly X{ 1 };
        int k = 1;
        while (k < n) {
            k <<= 1;
            X = (X * (Poly{ 1 } - X.log(k) + modx(k))).modx(k);
        }
        return X.modx(n);
    }
};

多項式冪函式

給定多項式 \(f(x)\) ，求 \(f^k(x) \bmod x^n\) 。

顯然有 \(f^k(x) \equiv e^{k \ln f(x)} \pmod{x^n}\) 。

指數並非真正的指數，而是多項式函式的自變數，因此指數上的 \(k\) 也是對 \(P\) 取模。

當 \([x^0]f(x) \neq 1\) 時，我們找到第一個非零項 \([x^{cnt}]f(x)\) ，多項式可以整體除以 \([x^{cnt}]f(x) \cdot x^{cnt}\) 再用上面的公式，最後乘 \(([x^{cnt}]f(x))^k \cdot x^{k \cdot cnt}\) ，其中 \([x^{cnt}]f(x)\) 的 \(k\) 次方要模 \(\varphi(P)\) ，因為他是真正意義的指數。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly pow(int k = 0, int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        if (1LL * k * cnt >= n) return Poly(n);
        int k1 = k % P, k2 = k % (P - 1);
        return ((k1 * (mulx(-cnt).modx(n) / (*this)[cnt]).log(n)).exp(n) * qpow((*this)[cnt], k2)).mulx(cnt * k1).modx(n);
    }
    Poly pow(const std::string &s, int n = 0) const {
        if (n == 0) n = size();
        int cnt = 0;
        while (cnt < size() && (*this)[cnt] == 0) cnt++;
        if (cnt == size()) return Poly(n);
        int k1 = 0, k2 = 0;
        for (auto ch : s) {
            if ((1LL * 10 * k1 + ch - '0') * cnt >= n) return Poly(n);
            k1 = (1LL * 10 * k1 % P + ch - '0') % P;
            k2 = (1LL * 10 * k2 % (P - 1) + ch - '0') % (P - 1);
        }
        return ((k1 * (mulx(-cnt).modx(n) / (*this)[cnt]).log(n)).exp(n) * qpow((*this)[cnt], k2)).mulx(cnt * k1).modx(n);
    }
};

多項式三角函式

給定多項式 \(f(x)\) ，求 \(\sin f(x) \bmod x^n\) 和 \(\cos f(x) \bmod x^n\) 。

考慮尤拉公式 \(e^{\text{i}\theta} = \cos \theta + \text{i} \sin \theta\) 。

因此有 \(\sin \theta = \dfrac{e^{\text{i} \theta} - e^{-\text{i} \theta}}{2 \text{i}} , \cos \theta = \dfrac{e^{\text{i} \theta} + e^{-\text{i} \theta}}{2}\) 。

令 \(\theta = f(x)\) 即可，其中 \(\text{i}\) 在模 \(P\) 下等價於 \(g^{\frac{P-1}{4}}\) ，注意 \(f(x)\) 常數項必須為 \(0\) ，否則無意義。

此外，用這兩個函式可以推匯出其他三角函式（但有些無法在 \(0\) 處展開，且在其他位置展開都存在超越數，就不存在於這個模多項式體系了），這裡就不一一列舉了。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly sin(int n = 0) const {
        if (n == 0) n = size();
        return ((I * modx(n)).exp(n) - (((P - I) % P) * modx(n)).exp(n)).modx(n) * qpow(2LL * I % P, P - 2);
    }
    Poly cos(int n = 0) const {
        if (n == 0) n = size();
        return ((I * modx(n)).exp(n) + (((P - I) % P) * modx(n)).exp(n)).modx(n) * qpow(2, P - 2);
    }
};

多項式反三角函式

給定多項式 \(f(x)\) ，求 \(\arcsin f(x) \bmod x^n\) 和 \(\arctan f(x) \bmod x^n\) 。

考慮求導再積分回去， \(\displaystyle \arcsin f(x) = \int \frac{f'(x)}{\sqrt{1 - f^2(x)}} \text{d}x, \arctan f(x) = \int \frac{f'(x)}{1 + f^2(x)} \text{d}x\) 。

為什麼沒有 \(\arccos\) ？因為他的多項式常數是超越數，在這個體系無意義。上面能求的積分出來的常數是 \(0\) 。

時間複雜度 \(O(n \log n)\)

空間複雜度 \(O(n)\)

/// 有限域多項式板子（部分）
struct Poly : public std::vector<int> {
	Poly asin(int n = 0) const {
        if (n == 0) n = size();
        return (derive(n) * (Poly{ 1 } - (modx(n) * modx(n))).sqrt(n).inv(n)).integral(n);
    }
    Poly atan(int n = 0) const {
        if (n == 0) n = size();
        return (derive(n) * (Poly{ 1 } + (modx(n) * modx(n))).inv(n)).integral(n);
    }
};

位運算卷積

位運算卷積的定義

快速沃爾什變換（FWT）

異或卷積

或卷積

與卷積

子集卷積

子集卷積的定義

快速莫比烏斯變換（FMT）

並卷積

交卷積

子集卷積

多項式分治

多項式多點求值

多項式快速插值

分治加法卷積

多項式線性齊次遞推

多項式平移

相關文章