目录
多项式
多项式基础
数域的定义
定义 复数集的子集 \(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}\) 。
单位根具有许多有用的性质：
- 互异性：若 \(i \neq j\) ，则 \(\omega_n^i \neq \omega_n^j\) 。
- 折半律：\(\omega_{2n}^{2i} = \omega_{n}^{i}\) 。
- 周期律：\(\omega_n^{i + n} = \omega_n^i\) 。
- 半周期律： \(\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\) 项多项式的变换过程：
- 初始层：\(\{a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7\}\) 。
- 第一层：\(\{a_0, a_2, a_4, a_6\},\{a_1, a_3, a_5, a_7\}\) 。
- 第二层：\(\{a_0, a_4\},\{a_2, a_6\},\{a_1, a_5\}, \{a_3, a_7\}\) 。
- 第三层：\(\{a_0\}, \{a_4\},\{a_2\}, \{a_6\},\{a_1\}, \{a_5\}, \{a_3\}, \{a_7\}\) 。
我们需要从下往上迭代，那么就需要知道最后一层的顺序。
容易知道，\(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);
}
};