多项式

目录

• 多项式
  • 多项式基础
    • 数域的定义
    • 多项式的定义与基本性质
    • 多项式带余式除法
    • 形式幂级数的定义
    • 幂级数的导数和不定积分
    • 常见幂级数展开
  • 多项式插值
    • 多项式插值的定义
    • 多项式插值的方法
      • 拉格朗日插值法
      • 重心拉格朗日插值法
  • 加法卷积
    • 加法卷积的定义
    • 加法卷积的变换
      • 快速傅里叶变换（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）

并卷积

交卷积

子集卷积

多项式分治

多项式多点求值

多项式快速插值

分治加法卷积

多项式线性齐次递推

多项式平移