（v4 更新）OI+ACM 笔记：0x30 数学知识
0x31 数学知识 多项式
拉格朗日插值
多项式的插值：给出 \(n + 1\) 个点 \((x_1, y_1), (x_2, y_2), \cdots, (x_{n + 1}, y_{n + 1})\)，其中横坐标两两不同。求一个至多 \(n\) 次的多项式，使得这 \(n + 1\) 个点都在 \(f(x)\) 上。
拉格朗日插值：
构造方法：考虑构造 \(n + 1\) 个函数 \(f_1(x), \cdots, f_n(x), f_{n + 1}(x)\)，其中第 \(i\) 个函数 \(f_i(x)\) 的图像经过 \((x_i, y_i)\)，对于其他的 \(j \neq i\) 经过 \((x_j, 0)\)。则我们要求的 \(f(x)\) 即为这 \(n + 1\) 个函数的和 \(\sum_{i = 1}^{n + 1} f_i(x)\)。
尝试构造 \(f_i(x) = a_i \prod\limits_{j \neq i}(x - x_j)\)，其中 \(a_i\) 是常数。
这样构造可以满足，对于其他的 \(j \neq i\) 均有 \(f_i(x_j) = 0\)。代入 \(f_i(x_i) = y_i\)，解得 \(a_i = \frac{y_i}{\prod_{j \neq i}(x_i - x_j)}\)。
回代即可完成构造，不难发现 \(f(x)\) 是一个至多 \(n\) 次的多项式。并且可以证明拉格朗日插值是唯一的插值多项式。
朴素实现的时间复杂度是 \(\mathcal{O}(n^2)\)，可以使用多项式优化到 \(\mathcal{O}(n \log^2 n)\)（详见多项式快速插值）。
0x31 拉格朗日插值.cpp：
/* 拉格朗日插值 */
int lagrange(const std::vector< std::pair<int, int> > &seq, int k) {
int n = seq.size();
int ans = 0;
for (int i = 0; i < n; i ++) {
int p = 1, q = 1; // 先求出分子与分母
for (int j = 0; j < n; j ++) {
if (i == j) continue;
mul(p, k - seq[j].first);
mul(q, seq[i].first - seq[j].first);
}
ans = (ans + 1ll * p * qpow(q, mod - 2, mod) % mod * seq[i].second) % mod;
}
return (ans + mod) % mod;
}
拉格朗日插值求系数：相当于是直接将整个多项式求出来。运用到了暴力除以单项式的技巧。
- 记系数部分 \(a_i = \frac{y_i}{\prod_{j \neq i}(x_i - x_j)}\)。记多项式 \(g(x) = \prod_j(x - x_j)\)，记多项式 \(h_i(x) = \prod_{j \neq i} (x - x_j)\)。
- 先暴力多项式乘法 \(\mathcal{O}(n^2)\) 处理出 \(g(x)\)，考虑由 \(g(x)\) 得到 \(h_i(x)\)，相当于是多项式 \(g(x)\) 除以一个单项式 \((x - x_i)\)，仍然可以暴力 \(\mathcal{O}(n)\) 除以单项式。具体地，由于 \((x - x_i)h_i(x) = g(x)\)，则对于 \(g(x)\) 的 \(j(j \geq 1)\) 次项系数，有 \([x^{j - 1}]h_i - x_i[x^j]h_i = [x^j]g\)，即
- 求出了 \(a_i\) 与 \(h_i(x)\) 之后，\(f(x) = \sum_{i = 1}^{n + 1} a_ih_i(x)\)。时间复杂度 \(\mathcal{O}(n^2)\)。
0x31 拉格朗日插值求系数.cpp：
/* 拉格朗日插值求系数 */
std::vector<int> lagrange(const std::vector< std::pair<int, int> > &seq) {
int n = seq.size();
/* 求系数 a[] */
std::vector<int> a(n);
for (int i = 0; i < n; i ++) {
int q = 1;
for (int j = 0; j < n; j ++) {
if (i == j) continue;
mul(q, seq[i].first - seq[j].first);
}
if (q < 0) q += mod;
a[i] = 1ll * seq[i].second * qpow(q, mod - 2, mod) % mod;
}
/* 求总的多项式 g */
std::vector<int> g(n + 1);
g[0] = 1;
for (int i = 0; i < n; i ++) {
for (int j = i + 1; j >= 1; j --) {
g[j] = (1ll * g[j] * (mod - seq[i].first) + g[j - 1]) % mod;
}
g[0] = 1ll * g[0] * (mod - seq[i].first) % mod;
}
std::vector<int> f(n); // 答案多项式 f
for (int i = 0; i < n; i ++) {
std::vector<int> h(n); // 求分的多项式 h_i
int inv = qpow(mod - seq[i].first, mod - 2, mod);
h[0] = 1ll * g[0] * inv % mod;
for (int j = 1; j < n; j ++) {
h[j] = 1ll * (g[j] - h[j - 1] + mod) * inv % mod;
}
for (int j = 0; j < n; j ++) {
f[j] = (f[j] + 1ll * a[i] * h[j]) % mod;
}
}
return f;
}
int func(const std::vector<int> &f, int k) {
int n = f.size();
int y = 0;
for (int i = n - 1; i >= 0; i --) {
y = (1ll * y * k + f[i]) % mod;
}
return y;
}
连续点值的拉格朗日插值：给出 \(n + 1\) 个连续点值 \(f(1), f(2), \cdots, f(n + 1)\)，代入拉格朗日插值公式，得
现在考虑给定一个自变量 \(x\)，求出因变量 \(f(x)\)。考虑右边的乘积项怎么求。
- 分子：记 \(\mathrm{pre}_i = \prod_{1 \leq j \leq i}(x - j), \mathrm{suf}_i = \prod_{i\leq j \leq n + 1}(x - j)\)，则分子显然为 \(\mathrm{pre}_{i - 1} \times \mathrm{suf}_{i + 1}\)。
- 分母：对于 \(j < i\) 的部分，为 \(i - 1\) 到 \(1\) 乘起来；对于 \(j > i\) 的部分，为 \(-1\) 到 \(-(n + 1 - i)\) 乘起来。故分母为 \((i - 1)! (n + 1 - i)!(-1)^{n + 1 - i}\)。
故最后的式子为
求出单个点值的时间为 \(\mathcal{O}(n)\)。
若给出 \(n + 1\) 个连续点值 \(f(h), f(h + 1), \cdots, f(h + n)\)，分母的求法是完全不变的，只需将分子改成 \(x - x_j\) 相乘即可。
或者可以考虑将坐标系向左平移 \(h - 1\)。
0x31 连续点值的拉格朗日插值.cpp：
/* 连续点值的拉格朗日插值 */
int lagrange_conti(const std::vector<int> &y, int k) {
int n = y.size() - 1; // y 的下标从 1 开始
if (k <= n) {
return y[k];
}
std::vector<int> pre(n + 2), suf(n + 2);
pre[0] = 1;
for (int i = 1; i <= n; i ++) {
pre[i] = 1ll * pre[i - 1] * (k - i) % mod;
}
suf[n + 1] = 1;
for (int i = n; i >= 1; i --) {
suf[i] = 1ll * suf[i + 1] * (k - i) % mod;
}
std::vector<int> fact(n + 1), facv(n + 1);
fact[0] = 1;
for (int i = 1; i <= n; i ++) {
fact[i] = 1ll * fact[i - 1] * i % mod;
}
facv[n] = qpow(fact[n], mod - 2, mod);
for (int i = n - 1; i >= 0; i --) {
facv[i] = 1ll * facv[i + 1] * (i + 1) % mod;
}
int ans = 0;
for (int i = 1; i <= n; i ++) {
int cur = y[i];
cur = 1ll * cur * pre[i - 1] % mod * suf[i + 1] % mod;
cur = 1ll * cur * facv[i - 1] % mod * facv[n - i] % mod;
if ((n - i) & 1) {
dec(ans, cur);
} else {
add(ans, cur);
}
}
return ans;
}
FFT（快速傅里叶变换）
单位根：\(n\) 次单位根，即为 \(x^n = 1\) 的所有复数解。单位圆上辐角为 \(0, \frac{2\pi}{n}, \frac{4\pi}{n}, \cdots, \frac{2(n - 1)\pi}{n}\) 的复数都是单位根，共 \(n\) 个。这些单位根依次表示为
单位根的简单性质：欧拉公式：\(e^{ix} = \cos(x) + i \sin(x)\)。这说明 \(x \in \R\) 时，\(e^{ix}\) 可以描述一个单位圆上辐角为 \(x\) 的复数。
- \(\omega_n^1 = e^{\frac{2\pi i}{n}} = \cos(\frac{2\pi}{n}) + i \sin(\frac{2\pi}{n})\)。
- \(\omega_n^i \times \omega_n^j = \omega_n^{i + j}\)。
- \(\left(\omega_n^i\right)^j = \omega_n^{ij}\)。
- \(\omega_{2n}^{2k} = \omega_n^k\)。
- \(\omega_{2n}^k = -\omega_{2n}^{k + n}\)。
补充界：在进行 FFT 前，需要对其最高位补 \(0\)，将其补成一个最高次数为 \(2^P - 1\) 的多项式。方便后续的奇偶性分组。
DFT
DFT 的思想："系数表示法" 转 "点值表示法"。取自变量 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n - 1}\)，计算因变量 \(F(\omega_n^0), F(\omega_n^1), \cdots, F(\omega_n^{n - 1})\)。
记 \(n = 2^P\)，考虑一个 \(2^P - 1\) 次的多项式 \(F(x)\)
将每一项按照次数的奇偶性，划分成两部分（由于 \(n\) 是 \(2\) 的若干次幂，每次分出来的奇偶部分的次数一定相等）。
记
则
考虑计算因变量 \(F(\omega_n^0), F(\omega_n^1), \cdots, F(\omega_n^{n - 1})\)。
- 对于 \(k < n / 2\)，代入 \(\omega_n^k\) 得
- 对于 \(k < n / 2\)，代入 \(\omega_n^{k + n/2}\) 得
注意到 \(F(\omega_n^k), F(\omega_n^{k + n / 2})\) 等于 \(F_L(\omega_{n / 2}^k) \pm \omega_n^kF_R(\omega_{n / 2}^k)\)。
如果我们知道 \(F_L(x), F_R(x)\) 分别在 \(\omega_{n / 2}^0, \omega_{n / 2}^1, \cdots, \omega_{n / 2}^{n / 2 - 1}\) 处的点值表示，就可以 \(\mathcal{O}(n)\) 求出 \(F(x)\) 在 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n - 1}\) 处的点值表示。
这是一个分治的过程，不断分治直到仅剩一个项即可。
时间复杂度 \(\mathcal{O}(n \log n)\)。
IDFT
结论：将 DFT 中的 \(\omega_n^1\) 换成 \(\omega_n^{-1}\)，做完以后除以 \(n\) 即可。
IDFT 的思想："点值表示法" 转 "系数表示法"。已知因变量 \(F(\omega_n^0), F(\omega_n^1), \cdots, F(\omega_n^{n - 1})\)，求出自变量 \(\omega_n^0, \omega_n^1, \cdots, \omega_n^{n - 1}\)。
记 \(n = 2^P\)，对于一个 \(2^{P} - 1\) 的多项式 \(F(x)\)，设其经过 DFT 之后，得到的点值序列 \(g_0, g_1, \cdots, g_{n - 1}\)，可以证明：
证明：将 \(g_i\) 代入该式子，得
后面括号里的和式，即为单位根反演的经典结论 \([n \mid k] = \frac{1}{n}\sum_{i = 0}^{n - 1}\omega_{n}^{ik}\)。具体地，分类讨论：
- 当 \(j = k\) 时，和式为 \(\frac{1}{n}\sum_{i = 0}^{n - 1}\omega_n^0 = 1\)。
- 当 \(j \neq k\) 时，和式为
故 \(a_k = \frac{1}{n} \sum_{i = 0}^{n - 1}g_i \left( \omega_n^{-k} \right)^i\) 成立。
进一步，观察 IDFT 的式子，相当于将 DFT 中的 \(\omega_n^1\) 换成 \(\omega_n^{-1}\)，做完以后再除以 \(n\)。
再进一步，发现 \(F(\omega_n^{-i}) = F(\omega_n^{n - i})\)，于是 IDFT 可以看作将 DFT 后得到的点值进行翻转（注意 \(F(\omega_n^0)\) 这一项不参与翻转），然后再除以 \(n\)。
套用 DFT 的板子即可。
*DFT 与 IDFT 的矩阵角度
对于 DFT 前的系数表示法 \(a_0, a_1, \cdots, a_{n - 1}\) 以及 DFT 后的点值表示法 \(g_0, g_1, \cdots, g_{n - 1}\)。可以写成一个系数向量左乘一个转移矩阵得到点值向量的式子：
其中转移矩阵（从 \(0\) 开始的）第 \(i\) 行第 \(j\) 列的系数即为因变量 \(\omega_n^i\) 的 \(j\) 次方 \(\left(\omega_{n}^i\right)^j\)。这恰好是范德蒙德矩阵！
如何从点值向量得到系数向量？只需让上式等号两边都左乘一个转移矩阵的逆即可。结论：范德蒙德矩阵的逆，为各个系数取倒数，再除以矩阵大小 \(n\)。故有
这也解释了 IDFT 的结论：将 DFT 中的 \(\omega_n^1\) 换成 \(\omega_n^{-1}\)，做完以后除以 \(n\)。
FFT 倍增实现
位逆序置换：
以 \(n = 8\) 为例，模拟分治的过程：
从分治的角度，每次将 \(x_i\) 按照奇偶性分类，偶则分到左边，奇则分到右边。直到只剩下一个项为止。
从倍增的角度，我们希望对于每个项 \(x_i\)，求出其经过不断分治以后，只剩下一个项时，在哪一个位置。然后根据变换后的位置，不断地倍增合并点值。
将 \(x\) 用二进制数表示，变换后的位置即为该二进制数的翻转（例如 011000 变成 000110）。证明：每次根据奇偶性将 \(x_i\) 进行分类，相当于是由最低位决定最高位。以此类推，次低位决定次高位，第三低位决定第三高位 ...
设 \(\mathrm{rev}(x)\) 表示二进制数 \(x\) 的翻转。首先有 \(\mathrm{rev}(0) = 0\)，考虑从小到大递推，有递推式
蝶形运算：已知 \(F_L(\omega_{n / 2}^k), F_R(\omega_{n / 2}^k)\)。每次更新，需要用 \(F_L(\omega_{n / 2}^k) \pm \omega_n^kF_R(\omega_{n / 2}^k)\) 求出 \(F(\omega_n^k), F(\omega_n^{k + n / 2})\)。使用位逆序置换后
- 更新前：\(F_L(\omega_{n / 2}^k)\) 储存在数组下标为 \(k\) 的位置，\(F_R(\omega_{n / 2}^k)\) 储存在数组下标为 \(k + n / 2\) 的位置。
- 更新后：\(F(\omega_n^k)\) 储存在数组下标为 \(k\) 的位置，\(F(\omega_n^{k + n / 2})\) 储存在数组下标为 \(k + n / 2\) 的位置。
因此可以直接在数组下标为 \(k\) 和 \(k + n / 2\) 的位置进行覆写，而不用开额外的数组保存值。
FFT 多项式乘法
对于两个多项式 \(A, B\)，现在要求出多项式 \(C = A \times B\)。
考虑对 \(A, B\) 进行 DFT，经过 DFT 之后得到一系列关于 \(A, B\) 的点值，将关于 \(A, B\) 的点值对应位置相乘即可得到 \(C\) 的点值，再对 \(C\) 的点值进行 IDFT 即可得到 \(C\) 的系数。
核心思想：\(A, B\) 的系数转 \(A, B\) 的点值，\(A, B\) 的点值点积转 \(C\) 的点值，\(C\) 的点值转 \(C\) 的系数。
时间复杂度 \(\mathcal{O}((n + m) \log (n + m))\)。
0x31 FFT.cpp：
const double pi = acos(-1);
using comp = std::complex<double>;
std::vector<int> rev;
void dft(std::vector<comp> &a) {
int n = a.size();
for (int i = 0; i < n; i ++) {
if (i < rev[i]) {
std::swap(a[i], a[rev[i]]);
}
}
for (int k = 1; k < n; k <<= 1) {
comp omega(cos(pi / k), sin(pi / k)); // 2k 阶单位根
for (int i = 0; i < n; i += (k << 1)) {
comp x(1, 0);
for (int j = 0; j < k; j ++, x *= omega) {
comp u = a[i + j], v = x * a[i + j + k];
a[i + j] = u + v, a[i + j + k] = u - v;
}
}
}
}
void idft(std::vector<comp> &a) {
int n = a.size();
std::reverse(a.begin() + 1, a.end());
dft(a);
for (int i = 0; i < n; i ++) {
a[i] /= n;
}
}
struct poly : public std::vector<comp> {
poly() : std::vector<comp>() {}
explicit constexpr poly(int n) : std::vector<comp>(n) {}
explicit constexpr poly(const std::vector<comp> &a) : std::vector<comp>(a) {}
constexpr poly(const std::initializer_list<comp> &a) : std::vector<comp>(a) {}
template <class InputIt, class = std::_RequireInputIter<InputIt>>
explicit constexpr poly(InputIt st, InputIt ed) : std::vector<comp>(st, ed) {}
constexpr friend poly operator * (poly a, poly b) {
int tot = a.size() + b.size() - 1;
if (tot < 128) {
poly c(tot);
for (int i = 0; i < a.size(); i ++) {
for (int j = 0; j < b.size(); j ++) {
c[i + j] += a[i] * b[j];
}
}
return c;
}
int L = 1, P = 0;
while (L < tot) L <<= 1, P ++;
rev.resize(L);
for (int i = 0; i < L; i ++) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (P - 1));
}
a.resize(L), b.resize(L);
dft(a), dft(b);
for (int i = 0; i < L; i ++) {
a[i] *= b[i];
}
idft(a);
a.resize(tot);
return a;
}
};
NTT（快速数论变换）
回顾一些内容：
阶：设整数 \(m > 0\)，满足 \(a^x \equiv 1 \pmod m\) 的最小正整数 \(x\)，称作 \(a\) 模 \(m\) 的阶。记作 \(\mathrm{ord}_m(a)\)。
阶的存在条件：当且仅当 \(\gcd(a, m) = 1\) 时，\(a\) 模 \(m\) 的阶存在。
阶的简单性质：
- 若 \(a\perp m\)，则 \(a, a^2, \cdots, a^{\mathrm{ord}_m(a)}\) 模 \(m\) 两两不同。
- 若 \(a\perp m\)，则 \(\mathrm{ord}_m(a) \mid \varphi(m)\)。
- 若 \(a \perp m\)，则
原根：设整数 \(m > 0\)，若 \(\mathrm{ord}_m(g) = \varphi(m)\)，则称 \(g\) 为模 \(m\) 的原根。
在模意义下， 我们希望找到一个东西来代替单位根。主要是要满足 FFT 过程中，用到的单位根满足的三个性质：
- \(\omega_n\) 的阶为 \(n\)。
- \(\omega_{2n}^{2k} = \omega_n^k\)。
- \(\omega_{2n}^k = -\omega_{2n}^{n + k}\)。
对于一个质数 \(m\)，考虑 \(m\) 的一个原根 \(g\)。此时 \(g\) 的阶为 \(m - 1\)，显然原根 \(g\) 不能直接代替 \(n\) 阶单位根。
考虑 \(g\) 的若干次幂。\(g^k\) 的阶为 \(\frac{m - 1}{\gcd(m - 1, k)}\)，必为 \(m - 1\) 的约数。且当 \(n \mid (m - 1)\) 时，可以构造 \(g^{\frac{m - 1}{n}}\) 的阶恰好等于 \(n\)。
尝试使用 \(g^{\frac{m - 1}{n}}\) 代替 \(n\) 阶单位根，性质 1, 2 的证明显然，性质 3 相当于是要证明 \(g^{\frac{m - 1}{2}} \equiv -1 \pmod m\)，可以使用二次剩余中的欧拉判别法来证明（原根必定不是二次剩余）。
于是当 \(n \mid (m - 1)\) 时，我们完全可以使用 \(g^{\frac{m - 1}{n}}\) 来代替 \(n\) 阶单位根 \(\omega_n\)。
现在还要满足 \(n \mid (m - 1)\) 才可以代替，对模数 \(m\) 的取值有特殊的要求。注意到 FFT 过程中的 \(n\) 都是 \(2\) 的若干次幂，所以 \(m - 1\) 只需要包含大量的因子 \(2\) 即可。例如 \(998244353 = 2^{23} \times 7 \times 17 + 1\)，就是一个经典的 NTT 模数。
常用的 NTT 模数：
998244353：原根为 \(3\)。1004535809：原根为 \(3\)。469762049：原根为 \(3\)。167772161：原根为 \(3\)。924844033：原根为 \(5\)。4179340454199820289（\(4 \times 10^{18}\) 级别，以时间换精度）：原根为 \(3\)。
0x31 NTT.cpp：
// 预处理单位根
constexpr auto findRootPower() {
int g = 3; // 模数对应的原根，需根据实际问题调整
std::array<int, 32> w{};
for (int k = 1, idx = 0; (mod - 1) % (k << 1) == 0; k <<= 1) {
w[idx ++] = qpow(g, (mod - 1) / (k << 1), mod); // 2k 阶单位根
}
return w;
}
constexpr auto RootPower = findRootPower();
std::vector<int> rev;
void dft(std::vector<int> &a) {
int n = a.size();
for (int i = 0; i < n; i ++) {
if (i < rev[i]) {
std::swap(a[i], a[rev[i]]);
}
}
for (int k = 1, idx = 0; k < n; k <<= 1) {
int omega = RootPower[idx ++];
for (int i = 0; i < n; i += (k << 1)) {
int x = 1;
for (int j = 0; j < k; j ++, mul(x, omega)) {
int u = a[i + j], v = 1ll * x * a[i + j + k] % mod;
add(a[i + j] = u, v), dec(a[i + j + k] = u, v);
}
}
}
}
void idft(std::vector<int> &a) {
int n = a.size(), inv = qpow(n, mod - 2, mod);
std::reverse(a.begin() + 1, a.end());
dft(a);
for (int i = 0; i < n; i ++) {
mul(a[i], inv);
}
}
struct poly : public std::vector<int> {
poly() : std::vector<int>() {}
explicit constexpr poly(int n) : std::vector<int>(n) {}
explicit constexpr poly(const std::vector<int> &a) : std::vector<int>(a) {}
constexpr poly(const std::initializer_list<int> &a) : std::vector<int>(a) {}
template <class InputIt, class = std::_RequireInputIter<InputIt>>
explicit constexpr poly(InputIt st, InputIt ed) : std::vector<int>(st, ed) {}
constexpr friend poly operator * (poly a, poly b) {
int tot = a.size() + b.size() - 1;
if (tot < 128) {
poly c(tot);
for (int i = 0; i < a.size(); i ++) {
for (int j = 0; j < b.size(); j ++) {
c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % mod;
}
}
return c;
}
int L = 1, P = 0;
while (L < tot) L <<= 1, P ++;
rev.resize(L);
for (int i = 0; i < L; i ++) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (P - 1));
}
a.resize(L), b.resize(L);
dft(a), dft(b);
for (int i = 0; i < L; i ++) {
mul(a[i], b[i]);
}
idft(a);
a.resize(tot);
return a;
}
};
MTT（任意模数 NTT）
当模数不为常用的 NTT 模数时，我们就需要使用一些技巧使得多项式乘法可以继续进行。
首先不可以简单地使用 FFT，因为 double 的有效精度只有 15~16 位，而对于常见的数据范围 \(n, m \leq 10^5\) 以及 \(a_i, b_i \leq 10^9\)，卷积结果的最大值接近 \(10^5 \times 10^9 \times 10^9 = 10^{23}\)，一定会掉精度。
三模 NTT
由于常见的数据范围下，卷积结果的最大值接近 \(10^{23}\)。于是我们考虑选择三个 \(10^9\) 级别的常用 NTT 模数，用三个模数分别做一次多项式乘法，将得到的三个结果使用 CRT 合并。
常用的 NTT 模数选法如下，因为这三个 NTT 模数均有一个原根为 \(3\)。
CRT 过程中需要使用 __int128。
0x31 MTT（三模 NTT）.cpp：
// 预处理单位根
constexpr auto findRootPower(int mod) {
int g = 3; // 模数对应的原根
std::array<int, 32> w{};
for (int k = 1, idx = 0; (mod - 1) % (k << 1) == 0; k <<= 1) {
w[idx ++] = qpow(g, (mod - 1) / (k << 1), mod); // 2k 阶单位根
}
return w;
}
std::vector<int> rev;
template <const int mod>
struct mulSpace {
std::array<int, 32> RootPower;
mulSpace() { RootPower = findRootPower(mod); }
inline void add(int &x, const int &y) { x += y; if (x >= mod) x -= mod; }
inline void dec(int &x, const int &y) { x -= y; if (x < 0) x += mod; }
inline void mul(int &x, const int &y) { x = 1ll * x * y % mod; }
void dft(std::vector<int> &a) {
int n = a.size();
for (int i = 0; i < n; i ++) {
if (i < rev[i]) {
std::swap(a[i], a[rev[i]]);
}
}
for (int k = 1, idx = 0; k < n; k <<= 1) {
int omega = RootPower[idx ++];
for (int i = 0; i < n; i += (k << 1)) {
int x = 1;
for (int j = 0; j < k; j ++, mul(x, omega)) {
int u = a[i + j], v = 1ll * x * a[i + j + k] % mod;
add(a[i + j] = u, v), dec(a[i + j + k] = u, v);
}
}
}
}
void idft(std::vector<int> &a) {
int n = a.size(), inv = qpow(n, mod - 2, mod);
std::reverse(a.begin() + 1, a.end());
dft(a);
for (int i = 0; i < n; i ++) {
mul(a[i], inv);
}
}
std::vector<int> mul(std::vector<int> a, std::vector<int> b, int L) {
for (int &x : a) x %= mod;
for (int &x : b) x %= mod;
a.resize(L), b.resize(L);
dft(a), dft(b);
for (int i = 0; i < L; i ++) {
mul(a[i], b[i]);
}
idft(a);
return a;
}
};
const int mod1 = 998244353, mod2 = 1004535809, mod3 = 469762049;
mulSpace<mod1> m1;
mulSpace<mod2> m2;
mulSpace<mod3> m3;
/* （完整地）得出 CRT 的解 */
int CRT(int x1, int x2, int x3) {
static const s128 M = (s128)mod1 * mod2 * mod3;
static const s64 M1 = M / mod1, M2 = M / mod2, M3 = M / mod3;
static const int
V1 = qpow(M1 % mod1, mod1 - 2, mod1),
V2 = qpow(M2 % mod2, mod2 - 2, mod2),
V3 = qpow(M3 % mod3, mod3 - 2, mod3);
s128 x = 0;
x += (s128)x1 * M1 * V1;
x += (s128)x2 * M2 * V2;
x += (s128)x3 * M3 * V3;
x %= M; // 此时的 x 即为 CRT 的解
return x % mod;
}
struct poly : public std::vector<int> {
poly() : std::vector<int>() {}
explicit constexpr poly(int n) : std::vector<int>(n) {}
explicit constexpr poly(const std::vector<int> &a) : std::vector<int>(a) {}
constexpr poly(const std::initializer_list<int> &a) : std::vector<int>(a) {}
template <class InputIt, class = std::_RequireInputIter<InputIt>>
explicit constexpr poly(InputIt st, InputIt ed) : std::vector<int>(st, ed) {}
constexpr friend poly operator * (poly a, poly b) {
int tot = a.size() + b.size() - 1;
if (tot < 128) {
poly c(tot);
for (int i = 0; i < a.size(); i ++) {
for (int j = 0; j < b.size(); j ++) {
c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % mod;
}
}
return c;
}
int L = 1, P = 0;
while (L < tot) L <<= 1, P ++;
rev.resize(L);
for (int i = 0; i < L; i ++) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (P - 1));
}
std::vector<int> c1 = m1.mul(a, b, L), c2 = m2.mul(a, b, L), c3 = m3.mul(a, b, L);
poly c(tot);
for (int i = 0; i < tot; i ++) {
c[i] = CRT(c1[i], c2[i], c3[i]);
}
return c;
}
};
Trick：翻转卷积
给出序列 \(a_0, \cdots, a_n\) 与序列 \(b_0, \cdots, b_m\)，对于每一个 \(k\) 需要求出 \(c_k = \sum_{i} a_ib_{i - k}\)，即下标差恒为 \(k\) 的卷积式（相当于 \(a_ib_j\) 贡献给了 \(c_{i - j}\)）。
将序列 \(b\) 翻转得到 \(b'\)，其中 \(b'_{m - j} = b_j\)。对序列 \(a, b'\) 进行普通卷积得到序列 \(c'\)，此时 \(a_ib_j\) 贡献给了 \(c'_{m + i - j}\)，于是有 \(c_k = c'_{k + m}\)。
Trick：积卷积
给出序列 \(a_1, \cdots, a_{p - 1}\) 与序列 \(b_1, \cdots, b_{p - 1}\)，对于每一个 \(k \in [1, p)\) 需要求出 \(c_k = \sum_{i\cdot j \bmod p = k} a_ib_j\)。
当模数为较小的质数 \(p\) 时，可以使用 "离散对数"，将模 \(p\) 意义下的乘法转换为模 \(p - 1\) 意义下的加法。
Trick：多项式点乘向量
给出多项式 \(A(x) = \sum_{i = 0}^n a_ix^i\) 与向量 \(\mathbf{B} = (b_0, \cdots, b_n)\)，现在我们要求多项式 \(A\) 点乘向量 \(\mathbf{B}\) 的值
记翻转多项式 \(B^R(x) = \sum_{i = 0}^n b_{n - i}x^i\)，则点乘的值即为 \([x^n](A \times B^R)\)，即 \(A \times B^R\) 的其中一项。
当 \(A\) 可以分解成两个多项式乘积 \(A = A_1 \times A_2\) 时，若已知 \(A_1 \times B^R\) 以及 \(A_2\)，可以 \(\mathcal{O}(n)\) 求出 \([x^n](A \times B^R)\)。
（例如：给出多项式序列 \(p_1, \cdots, p_n\)，欲求 \(\left(\prod_{i = l}^r p_i\right) \cdot \mathbf{B}\)，即求 \(\mathrm{pre}_r \times \mathrm{ipre}_{l - 1} \times B^R\)。预处理出所有的 \(\mathrm{pre}_i \times B^R\) 以及 \(\mathrm{ipre}_i\) 即可）
Model：分治卷积
分治 NTT - 模型 1：给出 \(n\) 个数 \(a_1, \cdots, a_n\)，求出 \(n\) 个单项式的积 \(\prod_{i = 1}^n(a_ix + 1)\)。
设 \(f(l, r)\) 表示区间 \([l, r]\) 中所有单项式的积，\(f(l, r)\) 是一个次数为 \(r - l + 1\) 的单项式（与区间长度有关）。
取 \(\mathrm{mid} = \lfloor \frac{l + r}{2} \rfloor\)，则有
这是一个分治的过程，不断分治直到 \(l = r\) 即可。
时间复杂度 \(T(n) = T\left( \frac{n}{2} \right) + \mathcal{O}(n \log n) = \mathcal{O}(n \log^2 n)\)。
By the way：求出 \(n\) 个多项式的积，次数之和为 \(m\)。
可以每次取出两个次数最小的多项式相乘后放回，时间复杂度 \(\mathcal{O}(m \log^2 m)\)。
分治 NTT - 模型 2：给出 \(n\) 个数 \(a_1, \cdots, a_n\)，求出 \(n\) 个单项分式的和 \(\sum_{i = 1}^n \frac{1}{1 - a_ix}\)。
对于分式 \(\frac{A_1}{A_2}, \frac{B_1}{B_2}\)，若我们已知分子 \(A_1, B_1\) 与分母 \(A_2, B_2\)，由通分可得 \(\frac{A_1}{A_2} + \frac{B_1}{B_2} = \frac{A_1B_2 + A_2B_1}{A_2B_2}\)，由此可以得到新的分子 \(A_1B_2 + A_2B_1\) 与分母 \(A_2B_2\)。
设 \(A(l, r), B(l, r)\) 分别表示区间 \([l, r]\) 中所有单项分式的和，按照上述合并方法所得的分子与分母。\(A(l, r), B(l, r)\) 的次数同样与区间长度有关。
取 \(\mathrm{mid} = \lfloor \frac{l + r}{2} \rfloor\)，则有
通过分治求出 \(A(1, n), B(1, n)\)，最后求逆一次即可。
时间复杂度 \(\mathcal{O}(n \log^2 n)\)。
分治 NTT - 模型 3：给出 \(n\) 个数 \(a_1, \cdots, a_n\)，求出 \(\sum_{i = 1}^n \prod_{j \neq i} (a_jx + 1)\)。
设 \(f(l, r) = \sum_{i = l}^r \prod_{j \in [l, r], j \neq i} (a_jx + 1)\) 与 \(g(l, r) = \prod_{i = l}^r (a_ix + 1)\)，\(f(l, r), g(l, r)\) 的次数同样与区间长度有关。
取 \(\mathrm{mid} = \lfloor \frac{l + r}{2} \rfloor\)，则有
这是一个分治的过程，不断分治直到 \(l = r\) 即可。
时间复杂度 \(\mathcal{O}(n \log^2 n)\)。
Model：半在线卷积
半在线卷积：给出 \(g_1, \cdots, g_n\)，求出 \(f_0, \cdots, f_n\)，其中初值 \(f_0 = 1\)，递推式 \(f_i = \sum_{j = 1}^i f_{i - j}g_j\)。
使用强制中序遍历转移的 CDQ 分治：设 \(\mathrm{solve}(l, r)\) 表示求出区间 \([l, r]\) 中所有的 \(f\) 值，取 \(\mathrm{mid} = \lfloor \frac{l + r}{2} \rfloor\)。
- 先递归 \(\mathrm{solve}(l, \mathrm{mid})\)，求出 \(f_l, \cdots, f_{\mathrm{mid}}\)。
- 计算 \(f_{[l, \mathrm{mid}]}\) 与 \(g_{[0, r - l]}\) 的卷积结果，求出左半区间对右半区间的贡献。
- 后递归 \(\mathrm{solve}(\mathrm{mid} + 1, r)\)，求出 \(f_{\mathrm{mid} + 1}, \cdots, f_r\)。
- 最后递归到 \(l = r\) 时，请注意根据具体递推式调整 \(f_i\)（e.g. 递推式 \(f_i = g_i - \sum_{j = 1}^{i - 1} \binom{i - 1}{j - 1} f_{i - j}g_j\)）。
时间复杂度 \(T(n) = T\left( \frac{n}{2} \right) + \mathcal{O}(n \log n) = \mathcal{O}(n \log^2 n)\)。
卡常技巧：计算 \(f_{[l, \mathrm{mid}]}\) 与 \(g_{[0, r - l]}\) 的卷积结果时，实际上是 \(\frac{\mathrm{len}}{2}\) 次多项式与 \(\mathrm{len}\) 次多项式相乘。
由于 NTT 是循环卷积（\(a_ib_j \to c_{(i + j) \bmod L}\)，其中 \(L\) 表示补充界），将补充界设定成第一个 \(\geq \mathrm{len}\) 的 \(2\) 的幂，普通卷积的最后 \(\frac{\mathrm{len}}{2}\) 项会破坏掉循环卷积的前 \(\frac{\mathrm{len}}{2}\) 项，但由于我们只需要循环卷积的后 \(\frac{\mathrm{len}}{2}\) 项，所以将补充界调小不会对贡献产生影响。
全在线卷积：求出 \(f_0, \cdots, f_n\)，其中初值 \(f_0 = 1\)，递推式 \(f_i = \sum_{j = 1}^i f_{i - j}g_j\)，其中 \(g_i\) 与 \(f_0, \cdots, f_i\) 强相关。
仍然使用强制中序遍历转移的 CDQ 分治。但 \(g_{[0, r - l]}\) 我们不一定提前知道，进一步分析分治树的结构，发现仅有 \(l = 0\) 时不知道 \(g_{[0, r - l]}\)。所以我们需要修改计算贡献的方式。
- 当 \(l = 0\) 时，计算 \(f_{[0, \mathrm{mid}]} \times g_{[0, \mathrm{mid}]}\) 的贡献（此时还有 \(f_{[0, \mathrm{mid}]} \times g_{[\mathrm{mid} + 1, r]}\) 的贡献没被计算）。
- 当 \(l > 0\) 时，计算 \(f_{[l, \mathrm{mid}]} \times g_{[0, r - l]} + g_{[l, \mathrm{mid}]} \times f_{[0, r - l]}\) 的贡献（后半部分即为 \(l = 0\) 中漏掉的贡献）。
可以证明，点对的贡献不重不漏。
时间复杂度 \(\mathcal{O}(n \log^2 n)\)。
多项式除法（暴力）
多项式除法：给出两个多项式 \(A(x), B(x)\)，你需要求出多项式 \(C(x)\)，满足 \(B(x)C(x) = A(x)\)。
设 \(A(x) = \sum_{i = 0}^n a_ix^i, B(x) = \sum_{i = 0}^m b_ix^i, C(x) = \sum_{i = 0}^{n - m}c_ix^i\)。直接根据定义 \(\sum_{j = 0}^i c_jb_{i - j} = a_i\) 暴力展开（这里需要保证 \(b_0 \neq 0\)，否则需要对 \(b\) 的前导零进行处理）。
初值
递推式
但递推式里的 \(j\) 只需要枚举到 \(i - j \leq m\) 即 \(j \geq i - m\) 就可以了。
时间复杂度 \(\mathcal{O}(nm)\)。
多项式导数/积分
多项式形式导数：
多项式形式积分：
多项式牛顿迭代
问题：已知函数 \(g\)，假设满足 \(g(f(x)) = 0\) 的函数 \(f(x) = h(x)\)，我们需要求出 \(h \bmod x^n\)（在这里，我们将 \(f(x)\) 看成函数 \(g(f(x))\) 的自变量，则我们要求的就是 \(g(f(x)) = 0\) 模 \(x^n\) 意义下的零点）。
当 \(n = 1\) 时，\([x^0]h(x)\) 需要单独求出。
当 \(n > 1\) 时，记 \(h_0 = h \bmod x^{\left\lceil \frac{n}{2} \right\rceil}, h_1 = h \bmod x^n\)，假设我们现在已经求出了 \(h_0\)，我们希望通过 \(h_0\) 推出 \(h_1\)。考虑将函数 \(g(f(x))\) 在 \(f(x) = h_0(x)\) 处泰勒展开。
将 \(h_1(x)\) 代入自变量 \(f(x)\)，观察到 \(h_1(x) - h_0(x)\) 的前 \(\left\lceil \frac{n}{2} \right\rceil\) 项均为 \(0\)（因为模 \(x^{\left\lceil \frac{n}{2} \right\rceil}\) 意义下两者显然相等）。
也就是说 \(h_1(x) - h_0(x)\) 的最低次数为 \(\left\lceil \frac{n}{2} \right\rceil\)，进一步有 \(\left( h_1(x) - h_0(x) \right)^i\) 的最低次数为 \(\left\lceil \frac{n}{2} \right\rceil^i\)。当 \(i \geq 2\) 时，\(\left( h_1(x) - h_0(x) \right)^i \bmod x^n\) 均为 \(0\)，于是可以忽略 \(i \geq 2\) 的项。
进一步可以得到
移项可得
时间复杂度 \(T(n) = T\left(\frac{n}{2}\right) + \mathcal{O}(n \log n) = \mathcal{O}(n \log n)\)。
多项式求逆
多项式逆：给出一个多项式 \(A(x)\)，你需要求出多项式 \(B(x)\)，满足 \(A(x) \times B(x) = 1\)。一般我们只需要求出 \(B \bmod{x^m}\) 的结果。
多项式逆存在条件：多项式 \(A(x)\) 存在逆，当且仅当 \([x^0]A(x) \neq 0\)。
（注意：一个多项式的逆可能是有无穷项的，例如 \(\frac{1}{1 - x} = 1 + x + x^2 + \cdots\)）
少项式求逆（暴力）：
设 \(A(x) = \sum_{i = 0}^n a_ix^i, B(x) = \sum b_ix^i\)。直接根据多项式逆的定义 \(\sum_{j = 0}^{i} b_ja_{i - j} = [i = 0]\) 暴力展开。
初值
递推式
但递推式里的 \(j\) 只需要枚举到 \(i - j \leq n\) 即 \(j \geq i - n\) 就可以了。
时间复杂度 \(\mathcal{O}(nm)\)。当次数 \(n\) 较小时有奇效。
这也验证了 "多项式逆存在条件" 的充分性。
多项式求逆（牛顿迭代）：
构造函数 \(g(f(x)) = \frac{1}{f(x)} - A(x)\)，代入牛顿迭代法得
初始时，\(A^{-1}(x)\) 的常数项即为 \(A(x)\) 常数项的逆元。
时间复杂度 \(\mathcal{O}(m \log m)\)。
多项式开方
多项式开方：给出一个多项式 \(A(x)\)，你需要求出多项式 \(B(x)\)，满足 \(B^2(x) \equiv A(x)\)。一般我们只需要求出 \(B \bmod x^m\) 的结果。
多项式开方存在条件：多项式 \(A(x)\) 存在开方，当且仅当 \([x^0]A(x)\) 存在模意义下的二次剩余。
构造函数 \(g(f(x)) = f^2(x) - A(x)\)，代入牛顿迭代法得
初始时，\(\sqrt{A(x)}\) 的常数项即为 \(A(x)\) 常数项的二次剩余。当 \(A(x)\) 常数项不等于 \(1\) 时，还需要通过二次剩余求解常数项。
时间复杂度 \(\mathcal{O}(m \log m)\)。
By the way：注意到 \(\sqrt{A(x)} = e^{\frac{1}{2}\ln A(x)}\)，可以直接套用多项式 ln 与多项式 exp 的模板。
多项式带余除法
多项式除法 & 取模：给出一个 \(n\) 次多项式 \(A(x)\) 和一个 \(m\) 次多项式 \(B(x)\)，你需要求出多项式 \(Q(x), R(x)\) 满足：
- \(Q(x)\) 的次数为 \(n - m\)，\(R(x)\) 的次数 \(<m\)。
- \(A(x) = Q(x) \times B(x) + R(x)\)。
考虑构造变换 \(f^R(x) = x^{\mathrm{deg} f} f\left(\frac{1}{x}\right)\)，实际上即为翻转 \(f\) 的系数。
将 \(A(x) = Q(x) \times B(x) + R(x)\) 中的 \(x\) 替换为 \(\frac{1}{x}\)，并且等式两边同乘以 \(x^n\)
由于 \(Q(x)\) 的次数小于 \(n - m + 1\)，故对 \(x^{n - m + 1}\) 取模没有影响。
故在模 \(x^{n - m + 1}\) 意义下求出 \(\frac{A^R(x)}{B^R(x)}\) 再翻转一遍即可求出 \(Q(x)\)，再使用 \(A(x) - Q(x) \times B(x)\) 即可求出 \(R(x)\)。
时间复杂度 \(\mathcal{O}(n \log n)\)。
多项式 ln/exp
\(e\) 在模 \(998244353\) 下没有意义。同理，\(\ln 2\) 在模 \(998244353\) 下没有意义。
简单的解释：考察 \(e = \sum_{i = 0}^\infty \frac{1}{i!}\)，在 \(i \geq 998244353\) 时，分式就已经失去了意义。即使有意义，无穷项求和在模 \(998244353\) 下也不会收敛。
多项式 ln 与多项式 exp 的定义：由麦克劳林级数定义
多项式 ln 存在条件：多项式 \(f(x)\) 存在 ln，当且仅当 \([x^0]f(x) = 1\)，并且由此可得 \([x^0] \ln f(x) = 0\)。
多项式 ln：给出一个多项式 \(A(x)\)，你需要求出多项式 \(B(x) = \ln A(x)\)。一般我们只需要求出 \(B \bmod{x^m}\) 的结果。
将 \(B(x) = \ln A(x)\) 等式两边对 \(x\) 求导可得
两边同时积分可得
注意：虽然多项式先求导后积分会丢掉常数项的信息，但是容易验证当 \(A(x)\) 的常数项为 \(1\) 时，\(B(x)\) 的常数项必定为 \(0\)。
时间复杂度 \(\mathcal{O}(m \log m)\)。
多项式 exp 存在条件：多项式 \(f(x)\) 存在 exp，当且仅当 \([x^0]f(x) = 0\)，并且可以推出 \([x^0]\exp f(x) = 1\)。
多项式 exp：给出一个多项式 \(A(x)\)，你需要求出多项式 \(B(x) = \exp A(x)\)。一般我们只需要求出 \(B \bmod{x^m}\) 的结果。
构造函数 \(g(f(x)) = \ln f(x) - A(x)\)，代入牛顿迭代得
初始时，\(\exp A(x)\) 的常数项必为 \(1\)。
时间复杂度 \(\mathcal{O}(m \log m)\)。
多项式快速幂
多项式快速幂：给出一个多项式 \(A(x)\)，你需要求出 \(A^k(x) \bmod{x^m}\)。
多项式快速幂（倍增法）：具体算法流程同快速幂。
时间复杂度 \(\mathcal{O}(m \log m \log k)\)，但可以处理的问题形式更加灵活。
e.g. 多项式快速幂对另一多项式取模。
e.g. 多项式快速幂且保留次数在 \([-m, m]\) 的项。
多项式快速幂（先 ln 后 exp）：注意到 \(A^k(x) = e^{k\ln A(x)}\)，考虑直接套用多项式 ln 与多项式 exp 的模板。
但由于多项式 ln 的一些限制（需要保证 \([x^0]A(x) = 1\)），首先要对 \(A(x)\) 进行一些操作。
找出一个最小的 \(i\)，满足 \([x^i]A(x) \neq 0\)（相当于是找出 \(A(x)\) 的第一个非零系数），记 \(v = [x^i]A(x)\)。则有
注意到 \(\frac{A(x)}{x^iv}\) 满足常数项为 \(1\) 的限制，可以进行多项式 ln。
然后有
先使用多项式 ln 与多项式 exp 求出 \(\left(\frac{A(x)}{x^iv}\right)^k \bmod{x^{m - ik}}\) 的值，最后再乘上 \(x^{ik}v^k\) 即可。
（注意：这里的乘 \(x^{ik}\) 相当于多项式系数向右移动 \(ik\) 位；上面的除以 \(x^i\) 同理，相当于多项式系数向左移动 \(i\) 位，即舍弃掉前 \(i\) 项系数）
时间复杂度 \(\mathcal{O}(m \log m)\)。
设模数为质数 \(P\)。
- 当 \(k \leq P\) 时，可以放心应用该做法。
- 当 \(k > P\) 时，需要进一步讨论 \(k\) 的取模相关问题。
引理 1：设 \(P\) 为质数。则对于任意整数 \(a, b\)
\[(a + b)^P \equiv a^P + b^P \pmod{P} \]证明：因为 \((a + b)^P = \sum_{i = 0}^P \binom{P}{i}a^ib^{P - i}\)，当 \(0 < i < P\) 时，有 \(\binom{P}{i} = \frac{P}{i}\binom{P - 1}{i - 1}\)，由于 \(P\) 是质数，所以该式子是 \(P\) 的倍数。故上式得证。
引理 2：设 \(f(x)\) 为一个 \(n\) 次多项式，\(P\) 为质数。则有
\[f^P(x) \equiv f(x^P) \pmod{P} \]证明：不断地应用引理 1 即可。
考虑计算 \(f^P(x) \bmod x^m\)，由于一般情况下 \(m < p\)，所以 \(f^P(x) \equiv f(x^P) \equiv \left([x^0]f(x)\right)^P \pmod{x^m}\)。由于我们进行了一系列操作使得 \([x^0]f(x) = 1\)，所以 \(f^P(x) \equiv 1 \pmod{x^m}\)。
- 所以在多项式 ln 与多项式 exp 部分，\(k\) 需要对 \(P\) 取模。
- 考虑 \(x^{ik}\)，当 \(i = 0\) 时，可以直接忽略；当 \(i > 0\) 时，\(k\) 显然不能超过 \(P\)，否则右移 \(ik\) 位已经远超需要保留的范围。所以乘以 \(x^{ik}\) 部分，当 \(i > 0\) 时 \(k\) 不能超过 \(P\)。
- 考虑 \(v^k\)，由于是数值计算，所以乘以 \(v^k\) 部分，\(k\) 需要对 \(\varphi(P) = P - 1\) 取模。
少项式快速幂
少项式快速幂：给出一个次数不高的多项式 \(A(x)\)，你需要求出 \(A^k(x) \bmod x^m\)。
设 \(A(x) = \sum_{i = 0}^n a_ix^i, A^k(x) = \sum_{i = 0}^{m - 1}b_ix^i\)。考虑对 \(A^k(x)\) 求导
注意到 \(A^k(x)\) 与 \(A^k(x)'\) 有一定的递推关系。具体地，考虑等式左右两边的 \(i\) 次项
等式左边只需枚举到 \(j \geq i - n\)，等式右边只需枚举到 \(j \geq i - n + 1\)。可以解得
初值 \(b_0 = a_0^k\)，上式可以在已知 \(b_0, b_1, \cdots, b_i\) 的情况下递推出 \(b_{i + 1}\)。
需要线性预处理出 \(1, 2, \cdots, m - 1\) 的逆元，递推式：\(i^{-1} \equiv -\left\lfloor \frac{p}{i} \right\rfloor (p \bmod i)^{-1} \pmod{p}\)。
时间复杂度 \(\mathcal{O}(nm)\)。当次数 \(n\) 较小时有奇效。
By the way：从这个递推式进行归纳可知，当 \(a_0 = 1\) 时，\([x^n]A^k(x)\) 是一个关于 \(k\) 的 \(n\) 次多项式。
*转置原理
待填坑。
多项式多点求值
多项式多点求值：给出一个多项式 \(f(x)\) 和 \(n\) 个点 \(x_1, x_2, \cdots, x_n\)，求 \(f(x_1), f(x_2), \cdots, f(x_n)\)。
多项式取模做法
考虑使用分治将问题规模减半。设 \(\mathrm{mid} = \left\lfloor \frac{1 + n}{2} \right\rfloor\)，将给定的点集分成两部分
构造多项式
由多项式取模，记 \(f(x) = Q_1(x) \times g_1(x) + R_1(x)\)。对于 \(X_1\) 中的元素 \(x_i\)，由于 \(g_1(x_i) = 0\)，则有 \(f(x_i) = R_1(x_i)\)。又由于 \(R_1(x)\) 的次数小于区间长度的一半，于是问题规模成功减半。对于 \(X_2\) 中的元素也是同理。
多项式 \(g\) 可以通过分治 NTT 预处理。
时间复杂度 \(T(n) = 2T(n) + \mathcal{O}(n \log n) = \mathcal{O}(n \log^2 n)\)。
另一种理解角度：\(f(x) \equiv f(x_i) \pmod{x - x_i}\)。证明只需将 \(x^n = \left((x - x_i) + x_i\right)^n\) 二项式定理展开即可。
/* 多项式多点求值（多项式取模做法），作为 poly 的成员函数 */
constexpr std::vector<int> eval(std::vector<int> x) const {
int n = x.size();
if (size() == 0) {
return std::vector<int>(n, 0);
}
std::vector<poly> t(n * 4);
std::vector<int> ans(n);
std::function<void(int, int, int)> build = [&] (int p, int l, int r) {
if (l == r) {
int v = x[l];
t[p] = poly{v ? mod - v : 0, 1};
return;
}
int mid = (l + r) >> 1;
build(p * 2, l, mid), build(p * 2 + 1, mid + 1, r);
t[p] = t[p * 2] * t[p * 2 + 1];
};
build(1, 0, n - 1);
std::function<void(int, int, int, const poly &)> solve = [&] (int p, int l, int r, const poly &f) {
if (l == r) {
ans[l] = f[0];
return;
}
int mid = (l + r) >> 1;
auto [q1, r1] = divmod(f, t[p * 2]);
solve(p * 2, l, mid, r1);
auto [q2, r2] = divmod(f, t[p * 2 + 1]);
solve(p * 2 + 1, mid + 1, r, r2);
};
solve(1, 0, n - 1, *this);
return ans;
}
转置做法
待填坑。
多项式快速插值
多项式快速插值：给出 \(n\) 个点 \((x_1, y_1), (x_2, y_2), \cdots, (x_n, y_n)\)，其中横坐标两两不同。快速地求一个 \(n - 1\) 次多项式，使得这 \(n\) 个点都在 \(f(x)\) 上。
考虑拉格朗日插值
记 \(g(x) = \prod_{j = 1}^n(x - x_j)\)，则
首先使用分治 NTT 求出 \(g'(x)\)，再使用多项式多点求值，求出所有的 \(g'(x_i)\)。
记 \(v_i = \frac{y_i}{g'(x_i)}\)，则原式化为
仍然使用分治 NTT 求出上式。具体地，设 \(\mathrm{mid} = \left\lfloor \frac{1 + n}{2} \right\rfloor\)，设多项式
显然有 \(f(x) = f_1(x)g_2(x) + f_2(x)g_1(x)\)。\(f_1, f_2, g_1, g_2\) 可以进一步分治求解。
时间复杂度 \(\mathcal{O}(n \log^2 n)\)。
多项式平移
多项式平移：给出一个多项式 \(A(x)\) 以及一个常数 \(c\)，你需要求出多项式 \(A(x + c)\)。
多项式平移（二项式定理法）：设 \(A(x) = \sum_{i = 0}^n a_ix^i\)，考虑二项式定理
注意到括号里的求和式，下标之差恒为 \(j\)，于是可以使用 "翻转卷积" 的技巧求出每个关于 \(j\) 的求和式。
时间复杂度 \(\mathcal{O}(n \log n)\)。
连续点值平移：给出 \(n + 1\) 个连续点值 \(f(0), f(1), \cdots, f(n)\)（其中 \(f\) 是一个次数至多为 \(n\) 的多项式）以及一个常数 \(c\)，你需要求出 \(f(c), f(c + 1), \cdots, f(c + n)\)。
由于算法流程会涉及到 \(\frac{1}{c - n + i}\)（其中 \(i = 0 \sim 2n\)），我们不妨设 \(c > n\)。
对于 \(c \leq n\) 的情况，首先 \(f(c), \cdots, f(n)\) 是已知的，然后 \(f(n + 1), \cdots, f(n + c)\) 可以通过 \(f(0), \cdots, f(n)\) 向右平移 \(n + 1\) 然后取前 \(c\) 项获得（注意这里不能只取 \(f(0), \cdots, f(c - 1)\) 向右平移 \(n + 1\)，否则多项式的次数会变小）。
代入拉格朗日插值公式，得
于是
注意到右边的求和式，是一个卷积的形式。
记 \(a_i = \frac{1}{c - n + i}\)（其中 \(i = 0 \sim 2n\)，因为分母的取值范围为 \([c - n, c + n]\)），\(b_i = \frac{y_i(-1)^{n - i}}{i!(n - i)!}\)。则右边的求和式为 \(\sum_{i = 0}^n a_{x - i + n}b_i\)。
将 \(a, b\) 进行卷积得到 \(\mathrm{res}\)，则
时间复杂度 \(\mathcal{O}(n \log n)\)。
给出 \(n + 1\) 个等差点值 \(f(b), f(b + k), \cdots, f(b + nk)\)（其中 \(f\) 是一个次数至多为 \(n\) 的多项式） 以及一个常数 \(c\)，你需要求出 \(f(c + b), f(c + b + k), \cdots, f(c + b + nk)\)。
设 \(g(x) = f(b + xk)\)，可以看成：已知 \(g(0), g(1), \cdots, g(n)\)，向右平移 \(\frac{c}{k}\) 得到 \(g(\frac{c}{k}), g(\frac{c}{k} + 1), \cdots, g(\frac{c}{k} + n)\)。
应用：快速阶乘算法
快速阶乘算法：求 \(n! \bmod p\)，其中 \(p\) 是一个大质数（不保证是一个 NTT 模数的情况下，需要使用 MTT）。
当模数固定时，可以考虑分段打表。
设 \(v = \left\lfloor \sqrt{n} \right\rfloor\)，设多项式 \(g(x) = \prod_{i = 1}^v (x + i)\)，则
首先后半部分的 \(\prod_{i = v^2 + 1}^n i\) 可以 \(\mathcal{O}(\sqrt{v})\) 求出。我们现在只关心前半部分。
设多项式 \(g_d(x) = \prod_{i = 1}^d (x + i)\)，则 \(d + 1\) 个点值 \(g_d(0), g_d(v), \cdots, g_d(dv)\) 可以唯一确定 \(d\) 次多项式 \(g_d(x)\)。
所以对于当前的 \(d\)，我们都需要维护出这 \(d + 1\) 个点值。初始时 \(d = 1\)，只有点值 \(g_1(0) = 1, g_1(v) = v + 1\)。考虑两种操作：
操作 1：\(d \gets 2d\)
- 补齐点值：补齐 \(g_d((d + 1)v), \cdots, g_d(2dv)\)，可以通过 \(g_d(0), g_d(v), \cdots, g_d(dv)\) 向右平移 \(dv\) 得到。
（可以看作连续点值 \(h(x) = g(xv)\) 向右平移 \(d\)） - 修正点值：有公式 \(g_{2d}(x) = g_d(x)g_d(x + d)\)，将所有点值 \(g_d(x)\) 向右平移 \(d\) 即可得到 \(g_d(x + d)\)，相乘即可。
（可以看作连续点值 \(h(x) = g(xv)\) 向右平移 \(\frac{d}{v}\)）
操作 2：\(d \gets d + 1\)
- 修正点值：有公式 \(g_{d + 1}(x) = g_d(x) \cdot (x + d + 1)\)，直接相乘即可。
- 补齐点值：补齐 \(g_{d + 1}((d + 1)v)\)，直接暴力 \(\mathcal{O}(d)\) 求出即可。
于是可以从 \(d = 1\) 通过倍增变换到 \(d = v\)。倍增过程的时间复杂度为 \(T(n) = T(\frac{n}{2}) + \mathcal{O}(n \log n) = \mathcal{O}(n \log n)\)。
而我们只需要 \(\mathcal{O}(\sqrt{n})\) 个点值，故时间复杂度 \(\mathcal{O}(\sqrt{n} \log n)\)。
0x31 快速阶乘算法.cpp：
// 快速阶乘算法
int fastFactorial(int n) {
if (n >= mod) {
return 0;
}
int v = sqrt(n), iv = qpow(v, mod - 2, mod);
std::function< std::vector<int>(int) > solve = [&] (int d) {
if (d == 1) {
return std::vector<int>{1, v + 1};
}
if (d & 1) {
auto y = solve(d - 1);
for (int i = 0; i < y.size(); i ++) {
mul(y[i], i * v + d);
}
int num = 1;
for (int i = 1; i <= d; i ++) {
mul(num, d * v + i);
}
y.push_back(num);
return y;
} else {
auto y = solve(d / 2);
auto y1 = shift_conti(y, d / 2);
y.insert(y.end(), y1.begin() + 1, y1.end());
auto y2 = shift_conti(y, 1ll * (d / 2) * iv % mod);
for (int i = 0; i < y.size(); i ++) {
mul(y[i], y2[i]);
}
return y;
}
};
auto res = solve(v);
int ans = 1;
for (int i = 0; i < v; i ++) {
mul(ans, res[i]);
}
for (int i = v * v + 1; i <= n; i ++) {
mul(ans, i);
}
return ans;
}
⭐多项式全家桶
reference：jiangly 算法模板。
封装理念：
- 构建结构体
poly继承自std::vector<int>，自由度高且易于访问。 - 若函数的功能涉及 "当前多项式的变换或求值"，均作为
poly的内置函数；其余功能作为外置函数。
详细的使用手册位于代码注释之中！
0x31 多项式全家桶.cpp：
/**
* 多项式全家桶 使用手册
* author：Calculatelove
* version：2025.10.14
*
* 版本兼容：
* - C++23 以上：请放心使用
* - C++17 ~ C++20：
* 请删去所有 poly 中的 constexpr（在 poly 上方 #define constexpr 即可）
* - C++14：
* 请将 RootPower 前面的 constexpr 改成 const
* - C++11：
* 请删去所有函数体前面的 constexpr，
* 请将 findRootPower 与 RootPower 前面的 auto 改成 std::array<int, 32>
*
* 模数：
* - 涉及到多项式乘法时：
* 必须为 NTT 模数！否则请将 "多项式乘法" 替换为 "MTT 版本的多项式乘法"
* - 仅涉及到 poly 的四则运算、求导、积分，以及暴力多项式 mul_bf, div_bf, inv_bf，pow_bf 时：
* 只需模数为质数即可
*
* 原根：
* - 默认为 3，请根据具体 NTT 模数确认原根（在 findRootPower() 中修改）
*
* 已有的常数优化：
* - 使用了一些函数辅助取模（0x11 取模.cpp）
* - 使用 findRootPower() 预处理原根的幂（即单位根）
*
* 可能的常数优化：
* - 半在线卷积中，补充界只需 >= r - l + 1
* - 预处理阶乘以及阶乘逆元
* polyinit_factorial() -> shift(), shift_conti()
* - 预处理线性逆元
* polyinit_linearInverse() -> integ(), pow_bf()
*/
// 预处理单位根
constexpr auto findRootPower() {
int g = 3; // 模数对应的原根，需根据实际问题调整
std::array<int, 32> w{};
for (int k = 1, idx = 0; (mod - 1) % (k << 1) == 0; k <<= 1) {
w[idx ++] = qpow(g, (mod - 1) / (k << 1), mod); // 2k 阶单位根
}
return w;
}
constexpr auto RootPower = findRootPower();
std::vector<int> rev;
void dft(std::vector<int> &a) {
int n = a.size();
for (int i = 0; i < n; i ++) {
if (i < rev[i]) {
std::swap(a[i], a[rev[i]]);
}
}
for (int k = 1, idx = 0; k < n; k <<= 1) {
int omega = RootPower[idx ++];
for (int i = 0; i < n; i += (k << 1)) {
int x = 1;
for (int j = 0; j < k; j ++, mul(x, omega)) {
int u = a[i + j], v = 1ll * x * a[i + j + k] % mod;
add(a[i + j] = u, v), dec(a[i + j + k] = u, v);
}
}
}
}
void idft(std::vector<int> &a) {
int n = a.size(), inv = qpow(n, mod - 2, mod);
std::reverse(a.begin() + 1, a.end());
dft(a);
for (int i = 0; i < n; i ++) {
mul(a[i], inv);
}
}
/* 多项式全家桶 */
struct poly : public std::vector<int> {
poly() : std::vector<int>() {}
explicit constexpr poly(int n) : std::vector<int>(n) {}
explicit constexpr poly(const std::vector<int> &a) : std::vector<int>(a) {}
constexpr poly(const std::initializer_list<int> &a) : std::vector<int>(a) {}
template <class InputIt, class = std::_RequireInputIter<InputIt>>
explicit constexpr poly(InputIt st, InputIt ed) : std::vector<int>(st, ed) {}
// 多项式乘法
constexpr friend poly operator * (poly a, poly b) {
int tot = a.size() + b.size() - 1;
if (tot < 128) {
poly c(tot);
for (int i = 0; i < a.size(); i ++) {
for (int j = 0; j < b.size(); j ++) {
c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % mod;
}
}
return c;
}
int L = 1, P = 0;
while (L < tot) L <<= 1, P ++;
rev.resize(L);
for (int i = 0; i < L; i ++) {
rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (P - 1));
}
a.resize(L), b.resize(L);
dft(a), dft(b);
for (int i = 0; i < L; i ++) {
mul(a[i], b[i]);
}
idft(a);
a.resize(tot);
return a;
}
// 多项式乘法（暴力）
constexpr friend poly mul_bf(poly a, poly b) {
poly c(a.size() + b.size() - 1);
for (int i = 0; i < a.size(); i ++) {
for (int j = 0; j < b.size(); j ++) {
c[i + j] = (c[i + j] + 1ll * a[i] * b[j]) % mod;
}
}
return c;
}
// 多项式除法（暴力）
constexpr friend poly div_bf(poly a, poly b) { // 需要保证 b 的常数项不为 0
int m = b.size() - 1, inv = qpow(b[0], mod - 2, mod);
poly c(a.size() - b.size() + 1);
for (int i = 0; i < c.size(); i ++) {
c[i] = a[i];
for (int j = std::max(0, i - m); j < i; j ++) {
dec(c[i], 1ll * c[j] * b[i - j] % mod);
}
mul(c[i], inv);
}
return c;
}
constexpr friend poly operator + (poly a, poly b) {
poly c(std::max(a.size(), b.size()));
for (int i = 0; i < a.size(); i ++) {
add(c[i], a[i]);
}
for (int i = 0; i < b.size(); i ++) {
add(c[i], b[i]);
}
return c;
}
constexpr friend poly operator - (poly a, poly b) {
poly c(std::max(a.size(), b.size()));
for (int i = 0; i < a.size(); i ++) {
add(c[i], a[i]);
}
for (int i = 0; i < b.size(); i ++) {
dec(c[i], b[i]);
}
return c;
}
constexpr friend poly operator - (poly a) {
for (int &x : a) {
neg(x);
}
return a;
}
constexpr friend poly operator * (poly a, int b) {
for (int &x : a) {
mul(x, b);
}
return a;
}
constexpr friend poly operator / (poly a, int b) {
int inv = qpow(b, mod - 2, mod);
for (int &x : a) {
mul(x, inv);
}
return a;
}
constexpr poly &operator += (poly b) {
return (*this) = (*this) + b;
}
constexpr poly &operator -= (poly b) {
return (*this) = (*this) - b;
}
constexpr poly &operator *= (poly b) {
return (*this) = (*this) * b;
}
constexpr poly &operator *= (int b) {
return (*this) = (*this) * b;
}
constexpr poly &operator /= (int b) {
return (*this) = (*this) / b;
}
constexpr poly mulxk(int k) const {
auto res = *this;
res.insert(res.begin(), k, 0);
return res;
}
constexpr poly divxk(int k) const {
if (k >= size()) return poly();
return poly((*this).begin() + k, (*this).end());
}
constexpr poly modxk(int k) const {
if (k >= size()) return *this;
return poly((*this).begin(), (*this).begin() + k);
}
constexpr poly circ(int k) const {
auto res = *this;
for (int i = k; i < size(); i ++) {
add(res[i - k], res[i]);
}
return res.modxk(k);
}
constexpr poly flip() const {
auto res = *this;
std::reverse(res.begin(), res.end());
return res;
}
constexpr poly mulT(poly b) const {
if (b.size() == 0) return poly();
std::reverse(b.begin(), b.end());
return ((*this) * b).divxk(b.size() - 1);
// 等价于“翻转卷积”，a.mulT(b) 相当于是求 a[i] * b[j] -> c[i - j]
}
// 多项式导数
constexpr poly deriv() const {
if (empty()) return poly();
poly res(size() - 1);
for (int i = 1; i < size(); i ++) {
res[i - 1] = 1ll * (*this)[i] * i % mod;
}
return res;
}
// 多项式积分
constexpr poly integ() const {
poly res(size() + 1);
for (int i = 0; i < size(); i ++) {
res[i + 1] = 1ll * (*this)[i] * qpow(i + 1, mod - 2, mod) % mod;
}
return res;
}
// 多项式求逆
constexpr poly inv(int m) const { // 需要保证常数项不为 0
poly x{qpow((*this)[0], mod - 2, mod)};
int k = 1;
while (k < m) {
k <<= 1;
x = (x * (poly{2} - modxk(k) * x)).modxk(k);
}
return x.modxk(m);
}
// 少项式求逆
constexpr poly inv_bf(int m) const { // 需要保证常数项不为 0
auto &&a = *this;
int n = size() - 1, inv = qpow(a[0], mod - 2, mod);
poly b(m);
for (int i = 0; i < m; i ++) {
b[i] = i == 0;
for (int j = std::max(0, i - n); j < i; j ++) {
dec(b[i], 1ll * b[j] * a[i - j] % mod);
}
mul(b[i], inv);
}
return b;
}
// 多项式开方
constexpr poly sqrt(int m) const {
poly x{1}; // 默认常数项为 1
int k = 1;
while (k < m) {
k <<= 1;
x = (x + (modxk(k) * x.inv(k)).modxk(k)) / 2;
}
return x.modxk(m);
}
// 多项式带余除法
constexpr friend std::pair<poly, poly> divmod(poly a, poly b) {
int n = a.size(), m = b.size(), L = n - m + 1;
if (n < m) return {poly{0}, a};
poly q = (a.flip().modxk(L) * b.flip().inv(L)).modxk(L).flip();
poly r = a - b * q; while (r.size() && !r.back()) r.pop_back();
return {q, r};
}
// 多项式 ln
constexpr poly ln(int m) const { // 需要保证常数项为 1
return (deriv() * inv(m)).integ().modxk(m);
}
// 多项式 exp
constexpr poly exp(int m) const { // 需要保证常数项为 0
poly x{1};
int k = 1;
while (k < m) {
k <<= 1;
x = (x * (poly{1} - x.ln(k) + modxk(k))).modxk(k);
}
return x.modxk(m);
}
// 多项式快速幂
constexpr poly pow(int k, int m) const {
int i = 0;
while (i < size() && (*this)[i] == 0) {
i ++;
}
if (i == size() || 1ll * i * k >= m) {
return poly(m);
}
int v = (*this)[i];
poly f = divxk(i) * qpow(v, mod - 2, mod);
return (f.ln(m - i * k) * k).exp(m - i * k).mulxk(i * k) * qpow(v, k, mod);
/*
当实际 k < mod 时，可以放心使用
当实际 k >= mod 时，注意：
- 当 i > 0 时（即常数项为 0），直接返回长度为 m 的全 0 多项式
- 当 i = 0 时
- 在多项式 ln 与多项式 exp 部分，k 需要对 mod 取模
- 在 qpow(v, k, mod) 部分，k 需要对 mod - 1 取模
*/
}
// 少项式快速幂
constexpr poly pow_bf(int k, int m) const {
std::vector<int> v(m);
v[1] = 1;
for (int i = 2; i < m; i ++) {
v[i] = 1ll * v[mod % i] * (mod - mod / i) % mod;
}
auto &&a = *this;
int n = size() - 1, inv = qpow(a[0], mod - 2, mod);
poly b(m);
b[0] = qpow(a[0], k, mod);
for (int i = 0; i + 1 < m; i ++) {
b[i + 1] = 0;
for (int j = std::max(0, i - n + 1); j <= i; j ++) {
add(b[i + 1], 1ll * (i - j + 1) * a[i - j + 1] % mod * b[j] % mod);
}
mul(b[i + 1], k);
for (int j = std::max(0, i - n); j < i; j ++) {
dec(b[i + 1], 1ll * (j + 1) * b[j + 1] % mod * a[i - j] % mod);
}
mul(b[i + 1], v[i + 1]), mul(b[i + 1], inv);
}
return b;
}
// 多项式多点求值
constexpr std::vector<int> eval(std::vector<int> x) const {
if (size() == 0) {
return std::vector<int>(x.size(), 0);
}
int n = std::max(x.size(), this->size());
std::vector<poly> t(n * 4);
std::vector<int> ans(x.size());
x.resize(n);
std::function<void(int, int, int)> build = [&] (int p, int l, int r) {
if (l == r) {
t[p] = poly{1, norm(-x[l])};
return;
}
int mid = (l + r) >> 1;
build(p * 2, l, mid), build(p * 2 + 1, mid + 1, r);
t[p] = t[p * 2] * t[p * 2 + 1];
};
build(1, 0, n - 1);
std::function<void(int, int, int, const poly &)> solve = [&] (int p, int l, int r, const poly &num) {
if (l >= ans.size()) return;
if (l == r) {
ans[l] = num[0];
return;
}
int mid = (l + r) >> 1;
solve(p * 2, l, mid, num.mulT(t[p * 2 + 1]).modxk(mid - l + 1));
solve(p * 2 + 1, mid + 1, r, num.mulT(t[p * 2]).modxk(r - mid));
};
solve(1, 0, n - 1, mulT(t[1].inv(n)));
return ans;
}
// 多项式平移
constexpr poly shift(int c) const {
int n = size() - 1;
std::vector<int> fact(n + 1), facv(n + 1);
fact[0] = 1;
for (int i = 1; i <= n; i ++) {
fact[i] = 1ll * fact[i - 1] * i % mod;
}
facv[n] = qpow(fact[n], mod - 2, mod);
for (int i = n - 1; i >= 0; i --) {
facv[i] = 1ll * facv[i + 1] * (i + 1) % mod;
}
c = norm(c);
poly a(n + 1), b(n + 1);
for (int i = 0; i <= n; i ++) {
a[i] = 1ll * (*this)[i] * fact[i] % mod;
}
for (int i = 0, w = 1; i <= n; i ++, mul(w, c)) {
b[i] = 1ll * w * facv[i] % mod;
}
poly res = a.mulT(b);
for (int i = 0; i <= n; i ++) {
mul(res[i], facv[i]);
}
return res;
}
};
// 多项式快速插值
constexpr poly lagrange(std::vector< std::pair<int, int> > seq) {
int n = seq.size();
std::vector<poly> t(n * 4);
std::vector<int> x(n);
for (int i = 0; i < n; i ++) {
x[i] = seq[i].first;
}
std::function<void(int, int, int)> build = [&] (int p, int l, int r) {
if (l == r) {
t[p] = poly{norm(-x[l]), 1};
return;
}
int mid = (l + r) >> 1;
build(p * 2, l, mid), build(p * 2 + 1, mid + 1, r);
t[p] = t[p * 2] * t[p * 2 + 1];
};
build(1, 0, n - 1);
auto res = (t[1].deriv()).eval(x);
std::function<poly(int, int, int)> solve = [&] (int p, int l, int r) {
if (l == r) {
int v = 1ll * seq[l].second * qpow(res[l], mod - 2, mod) % mod;
return poly{v};
}
int mid = (l + r) >> 1;
return solve(p * 2, l, mid) * t[p * 2 + 1] + solve(p * 2 + 1, mid + 1, r) * t[p * 2];
};
return solve(1, 0, n - 1);
}
// 连续点值平移
constexpr std::vector<int> shift_conti(std::vector<int> y, int c) {
int n = y.size() - 1; // y 的下标从 0 开始
c = norm(c);
if (c == 0) return y;
if (c <= n) {
auto cur = std::vector<int>(y.begin() + c, y.end());
auto net = shift_conti(y, n + 1);
for (int i = 0; i < c; i ++) cur.push_back(net[i]);
return cur;
}
std::vector<int> fact(n + 1), facv(n + 1);
fact[0] = 1;
for (int i = 1; i <= n; i ++) {
fact[i] = 1ll * fact[i - 1] * i % mod;
}
facv[n] = qpow(fact[n], mod - 2, mod);
for (int i = n - 1; i >= 0; i --) {
facv[i] = 1ll * facv[i + 1] * (i + 1) % mod;
}
poly a(2 * n + 1), b(n + 1);
for (int i = 0; i <= 2 * n; i ++) {
a[i] = qpow(c - n + i, mod - 2, mod);
}
for (int i = 0; i <= n; i ++) {
b[i] = 1ll * facv[i] * facv[n - i] % mod * y[i] % mod;
if ((n - i) & 1) neg(b[i]);
}
poly res = a * b;
int num = 1;
for (int i = c - n; i <= c; i ++) mul(num, i);
std::vector<int> ans(n + 1);
for (int i = 0; i <= n; i ++) {
ans[i] = 1ll * num * res[i + n] % mod;
mul(num, c + i + 1), mul(num, qpow(c + i - n, mod - 2, mod));
}
return ans;
}
常系数齐次线性递推
常系数齐次线性递推：已知初值 \(a_0, \cdots, a_{k - 1}\) 与系数 \(c_1, \cdots, c_k\)（常数项为 \(0\)），对于一个 \(k\) 阶齐次线性递推数列 \(a_i = \sum_{j = 1}^k c_ja_{i - j}\)，你需要求出第 \(n\) 项。
Fiduccia 算法
Fiduccia 算法：构造特征多项式 \(\Gamma(x) = x^k - \sum_{i = 1}^{k}c_ix^{k - i}\)，记 \(A(x) = \sum_{i = 0}^{k - 1} a_ix^i\)，则有
（其中 \(\cdot\) 代表点积）
于是应用多项式快速幂（倍增法）求出 \(x^n \bmod \Gamma(x)\)，最后求出点积即可。
时间复杂度 \(\mathcal{O}(k \log k \log n)\)。
Fiduccia 算法的简单解释：我们想要知道 \(a_n\) 具体可以被多少个 \(a_0, \cdots, a_{k - 1}\) 表示。相当于是初始有一个 \(a_n\)，每次选出最高项 \(a_i\)，根据公式 \(a_i = \sum_{j = 1}^kc_ja_{i - j}\) 将 \(a_i\) 消成更低的项。
从生成函数的角度考虑，相当于是将 \(x^i\) 替换成 \(\sum_{j = 1}^kc_jx^{i - j}\)，等价于减去 \(x^{i - k}\left( x^k - \sum_{j = 1}^kc_jx^{k - j} \right)\)。整个过程等价于 \(x^n\) 对 \(x^{k} - \sum_{j = 1}^k c_jx^{k - j}\) 取模。
0x31 常系数齐次线性递推（Fiduccia）.cpp：
// 常系数齐次线性递推（Fiduccia）
int linearRecurrence(std::vector<int> c, std::vector<int> y, int n) {
int k = y.size(); // k 阶
if (n < k) {
return y[n];
}
poly g(k + 1);
g[k] = 1;
for (int i = 1; i <= k; i ++) {
g[k - i] = norm(-c[i]);
}
poly res{1}, a{0, 1};
for (int b = n; b; b >>= 1) {
if (b & 1) res = divmod(res * a, g).second;
a = divmod(a * a, g).second;
}
int ans = 0;
for (int i = 0; i < k; i ++) {
add(ans, 1ll * y[i] * res[i] % mod);
}
return ans;
}
Bostan-Mori 算法
Bostan-Mori 算法（求解分式远处系数）：该算法用于求解分式 \(\frac{P(x)}{Q(x)}\) 的第 \(n\) 项 \([x^n]\frac{P(x)}{Q(x)}\)。
考虑
记 \(V(x) = Q(x)Q(-x)\)，注意到 \(V(x) = V(-x)\)，则 \(V(x)\) 仅有偶次项有值。于是我们有
由于左边只影响奇次项，右边只影响偶次项（因为 \(\frac{1}{U(x^2)} = \frac{1}{Q(x)Q(-x)}\) 也是偶函数，仅有偶次项有值），于是我们只需要根据 \(n\) 的奇偶递归到其中一边即可。
当 \(n\) 为偶数时
当 \(n\) 为奇数时
递归到 \(n = 0\) 时，有 \([x^0]\frac{f(x)}{g(x)} = \frac{f_0}{g_0}\)。
常系数齐次线性递推与分式第 \(n\) 项的关系：对于数列 \(\{a_i\}\)，我们想要构造多项式 \(P(x), Q(x)\) 使得 \([x^n]\frac{P(x)}{Q(x)} = a_n\)。
取 \(Q(x) = \Gamma^R(x) = 1 - \sum_{i = 1}^k c_ix^i\) 以及 \(P(x) = ((\sum_{i = 0}^{k - 1} a_ix^i) \cdot Q(x)) \bmod x^k\)，即可满足 \([x^n]\frac{P(x)}{Q(x)} = a_n\)。
下面是检验，设 \(\frac{P(x)}{Q(x)} = \sum_{i = 0}^{\infty} h_ix^i\)，根据多项式除法的定义：
- 当 \(i = 0\) 时，有 \(h_0 = \frac{p_0}{q_0} = a_0\)。
- 当 \(1 \leq i < k\) 时，根据多项式 \(P, Q\) 的定义，显然有 \(h_i = a_i\)。
- 当 \(i \geq k\) 时，有 \(h_i = -\frac{1}{q_0}\sum_{j = 1}^k q_j h_{i - j}\)，化简得 \(h_i = \sum_{j = 1}^k c_jh_{i - j}\)，发现这是 \(\{a_i\}\) 的递推式，归纳得 \(h_i = a_i\)。
综上 \(h_i = a_i\)，原命题成立。
得出多项式 \(P, Q\) 后，套用上述算法求解 \([x^n]\frac{P(x)}{Q(x)}\) 即可。
时间复杂度 \(\mathcal{O}(k \log k \log n)\)，但常数相较于 Fiduccia 算法更小。
0x31 常系数齐次线性递推（Bostan-Mori）.cpp：
// 常系数齐次线性递推（Bostan-Mori）
int linearRecurrence(std::vector<int> c, std::vector<int> y, int n) {
int k = y.size();
if (n < k) {
return y[n];
}
// 求多项式 p/q 的第 n 项系数
std::function<int(poly, poly, int)> solve = [&] (poly p, poly q, int n) {
for (; n; n >>= 1) {
poly nq = q;
for (int i = 1; i < nq.size(); i += 2) {
neg(nq[i]);
}
p *= nq, q *= nq;
int i = n & 1;
while (i < p.size()) {
p[i / 2] = p[i], i += 2;
}
p.resize(i / 2);
int j = 0;
while (j < q.size()) {
q[j / 2] = q[j], j += 2;
}
q.resize(j / 2);
}
return 1ll * p[0] * qpow(q[0], mod - 2, mod) % mod;
};
poly q = poly{1} - poly(c);
poly p = (poly(y) * q).modxk(k);
return solve(p, q, n);
}
整式递推
待填坑。
拉格朗日反演
待填坑。
0x32 数学知识 集合幂级数
FWT（快速沃尔什变换）
位运算卷积：已知序列 \(\{a_i\}, \{b_i\}\)，求序列 \(\{c_i\}\) 满足 \(c_k = \sum_{i \oplus j = k} a_ib_j\)。其中 \(\oplus\) 表示一种位运算（or / and / xor）。
DWT
FFT 多项式乘法的核心思想：\(A, B\) 的系数转 \(A, B\) 的点值，\(A, B\) 的点值点积转 \(C\) 的点值，\(C\) 的点值转 \(C\) 的系数。
我们想要 DWT 与 IDWT 有着类似的功能。设 \(\mathrm{dwt}(a)\) 表示幂级数 \(a\) 经过 DWT 变换后得到的幂级数（点值），我们希望
由于 DFT 是线性变换，我们希望 DWT 也是线性变换。
记序列长度 \(n = 2^m\)，设 \(\mathrm{dwt}(a)_i = \sum_{j = 0}^{m - 1} c(i, j)a_j\)，其中变换系数 \(c(i, j)\) 表示变换前下标 \(j\) 对变换后下标 \(i\) 的贡献。
由 \(\mathrm{dwt}(a) \cdot \mathrm{dwt}(b) = \mathrm{dwt}(c)\)，得
由 \(a \times b = c\)，得
对比上下两段的最后一个式子，发现只需要满足 \(c(i, j)c(i, k) = c(i, j \oplus k)\) 即可。
另外，由于位运算可以按位考虑的特性，我们可以这样构造 \(c(i, j)\)：先得到单个位的转移系数 \(c([0, 1], [0, 1])\)，对于所有的 \(c(i, j)\)，我们将 \(i, j\) 二进制分解，将对应位的转移系数相乘
其中 \(i_x, j_x\) 分别表示 \(i, j\) 的二进制下第 \(x\) 位。
现在考虑优化计算 DWT。考虑按最高位拆分，记 \(i', j'\) 表示 \(i, j\) 去掉最高位的值，记 \(i_h\) 表示 \(i\) 的最高位
发现括号里的式子是一个子问题。记 \(a_0\) 表示幂级数 \(a\) 最高位为 \(0\) 的部分，\(a_1\) 表示幂级数 \(a\) 最高位为 \(1\) 的部分。
- 对于 \(k < n / 2\)，代入 \(i = k\) 得
- 对于 \(k < n / 2\)，代入 \(i = k + n / 2\) 得
如果我们知道 \(\mathrm{dwt}(a_0)_i, \mathrm{dwt}(a_1)_i\) 在 \(i = 0, 1, \cdots, n / 2 - 1\) 处的值，就可以 \(\mathcal{O}(n)\) 求出 \(\mathrm{dwt}(a)_i\) 在 \(i = 0, 1, \cdots, n - 1\) 处的值。
这是一个分治的过程，不断分治直到仅剩一个项即可。
时间复杂度 \(\mathcal{O}(m2^m)\)。
IDWT
在 DWT 中，我们只需要
这四个数即可完成变换，记该矩阵为位矩阵 \(c\)。
如果我们要进行 IDWT，则需要位矩阵 \(c\) 的逆矩阵 \(c^{-1}\)。类似地，有
套用 DWT 的模板即可。
在构造矩阵 \(c\) 的时候，我们需要小心一点，得确保矩阵 \(c\) 存在逆矩阵 \(c^{-1}\)。
位运算卷积
or 卷积：
or 运算的位矩阵
逆矩阵
or 卷积的正变换等价于子集求和，因为 \(c(i, j) = [i \operatorname{and} j = j]\)。
and 卷积：
and 运算的位矩阵
逆矩阵
and 卷积的正变换等价于超集求和，因为 \(c(i, j) = [i \operatorname{or} j = j]\)。
xor 卷积：
xor 运算的位矩阵
逆矩阵
xor 卷积的正变换，其变换系数 \(c(i, j) = (-1)^{\mathrm{popcount}(i \operatorname{and} j)}\)。
xor 卷积的 IDWT 相较于 DWT，每一位的变换系数需要乘以 \(\frac{1}{2}\)，一共有 \(m\) 位，总共的影响为 \(\frac{1}{2^m} = \frac{1}{n}\)。所以 xor 卷积的 IDWT 只需套用 DWT 的模板，最后乘以 \(\frac{1}{n}\) 即可。
0x32 FWT.cpp：
/* fwt-Or */
namespace bitOr {
void dwt(std::vector<int> &a, int opt) {
int n = a.size();
for (int k = 1; k < n; k <<= 1) {
for (int i = 0; i < n; i += (k << 1)) {
for (int j = 0; j < k; j ++) {
opt > 0 ? add(a[i + j + k], a[i + j]) : dec(a[i + j + k], a[i + j]);
}
}
}
}
std::vector<int> conv(std::vector<int> a, std::vector<int> b) {
dwt(a, +1), dwt(b, +1);
for (int i = 0; i < a.size(); i ++) {
mul(a[i], b[i]);
}
dwt(a, -1);
return a;
}
}
/* fwt-And */
namespace bitAnd {
void dwt(std::vector<int> &a, int opt) {
int n = a.size();
for (int k = 1; k < n; k <<= 1) {
for (int i = 0; i < n; i += (k << 1)) {
for (int j = 0; j < k; j ++) {
opt > 0 ? add(a[i + j], a[i + j + k]) : dec(a[i + j], a[i + j + k]);
}
}
}
}
std::vector<int> conv(std::vector<int> a, std::vector<int> b) {
dwt(a, +1), dwt(b, +1);
for (int i = 0; i < a.size(); i ++) {
mul(a[i], b[i]);
}
dwt(a, -1);
return a;
}
}
/* fwt-Xor */
namespace bitXor {
void dwt(std::vector<int> &a, int opt) {
int n = a.size();
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 = a[i + j], v = a[i + j + k];
add(a[i + j] = u, v), dec(a[i + j + k] = u, v);
}
}
}
if (opt < 0) {
int inv = qpow(n, mod - 2, mod);
for (int i = 0; i < n; i ++) {
mul(a[i], inv);
}
}
}
std::vector<int> conv(std::vector<int> a, std::vector<int> b) {
dwt(a, +1), dwt(b, +1);
for (int i = 0; i < a.size(); i ++) {
mul(a[i], b[i]);
}
dwt(a, -1);
return a;
}
}
子集 dp
子集 dp（Sum Over Subsets DP，SOS DP）：已知序列 \(\{a_i\}\)，求序列 \(\{b_i\}\) 满足
将一个二进制状态视作一个 \(m\) 维的坐标，第 \(i\) 维的取值等于二进制下第 \(i\) 位的值。
则二进制下的子集和，相当于是要求每个维度都小于等于当前坐标的所有坐标的权值和，运用高维前缀和即可。
时间复杂度 \(\mathcal{O}(m2^m)\)。
0x32 子集 dp.cpp：
// 注意：要先枚举维度 i，再枚举状态 S
for (int i = 0; i < m; i ++) {
for (int S = 0; S < (1 << m); S ++) {
if (S >> i & 1) {
f[S] += f[S ^ (1 << i)];
}
}
}
子集卷积
子集卷积：已知序列 \(\{a_i\}, \{b_i\}\)，求序列 \(\{c_i\}\) 满足
注意到 \(i \operatorname{and} j = \varnothing\) 时，必有 \(|i| + |j| = |i \operatorname{or} j|\)。于是引入序列 \(f\) 的占位多项式 \(F\)
将占位多项式代入到子集卷积的式子中
从 or 卷积的角度来看，有
于是将 \(A_p\) 和 \(B_q\) 进行 DWT，求出所有的 \(C_s\) 后再进行 IDWT 即可。
时空复杂度 \(\mathcal{O}(m^22^m)\)。
0x32 子集卷积.cpp：
/* 子集卷积 */
std::vector<int> subset_conv(std::vector<int> a, std::vector<int> b) {
int n = a.size(), m = std::__lg(n);
std::vector< std::vector<int> >
A(m + 1, std::vector<int>(n, 0)),
B(m + 1, std::vector<int>(n, 0)),
C(m + 1, std::vector<int>(n, 0));
for (int i = 0; i < n; i ++) {
A[__builtin_popcount(i)][i] = a[i];
B[__builtin_popcount(i)][i] = b[i];
}
for (int i = 0; i <= m; i ++) {
bitOr::dwt(A[i], +1), bitOr::dwt(B[i], +1);
}
for (int i = 0; i <= m; i ++) {
for (int j = 0; j <= i; j ++) {
for (int k = 0; k < n; k ++) {
add(C[i][k], 1ll * A[j][k] * B[i - j][k] % mod);
}
}
}
for (int i = 0; i <= m; i ++) {
bitOr::dwt(C[i], -1);
}
std::vector<int> c(n);
for (int i = 0; i < n; i ++) {
c[i] = C[__builtin_popcount(i)][i];
}
return c;
}
0x33 数学知识 组合数学
排列组合
排列数：\(A_{n}^m\)，从 \(n\) 个有标号元素中，取出 \(m\) 个元素排成一列的方案数（需要考虑取数顺序）。
排列数：\(C_n^m\) 或 \(\binom{n}{m}\)，从 \(n\) 个有标号元素中，取出 \(m\) 个元素组成集合的方案数（不需要考虑取数顺序）。
组合数常见恒等式
阶乘展开式：\(n, m \in \Z\)，\(n \geq m \geq 0\)
归纳恒等式：\(m \in \Z\)
对称恒等式：\(n, m \in \Z\)，\(n \geq 0\)
吸收恒等式：\(m \in \Z\)，\(m \neq 0\)
上指标反转：\(m \in \Z\)
三项式版恒等式：\(m, k \in \Z\)
二项式定理：\(n \in \Z\)
一些常见特例：
- \(x = 1\) 且 \(y = a\)：
\[\sum_{i = 0}^{n} \binom{n}{i}a^i = (1 + a)^n \]
- \(x = 1\) 且 \(y = -1\)：
\[\sum_{i = 0}^n (-1)^i \binom{n}{i} = [n = 0] \]
广义二项式定理：\(r \in \R\)
平行求和法：\(m \in \Z\)
上指标求和法：\(m \in \Z\)，\(n \geq 0\)
范德蒙德卷积：\(m \in \Z\)
斜对角线求和法：\(n \in \Z\)
等价于长度为 \(n - 1\) 的没有相邻 \(1\) 的二进制串个数。
组合数带权和：\(n \in \Z\)
组合数带权平方和：\(n \in \Z\)
插板法
问题 1：将 \(n\) 个完全相同的元素分成 \(k\) 组，保证每组至少有一个元素。等价于 \(\sum_{i = 1}^k x_i = n\) 且 \(\forall x_i \geq 1\) 的整数解个数。
相当于将 \(k - 1\) 块板，插入 \(n - 1\) 个间隙中。
问题 2：将 \(n\) 个完全相同的元素分成 \(k\) 组，允许组内元素为空。等价于 \(\sum_{i = 1}^k x_i = n\) 且 \(\forall x_i \geq 0\) 的整数解个数。
给 \(k\) 组分别补一个元素，转化为问题 1。
问题 3：给 \(n\) 个完全相同的元素分成 \(k\) 组，保证第 \(i\) 组至少分到 \(a_i\) 个元素。
给第 \(i\) 组提前安排 \(a_i\) 个元素，转化为问题 2。
Lucas 定理
Lucas 定理：对于质数 \(p\)
相当于是将 \(n, m\) 转化为 \(p\) 进制，将每一位上对应的组合数 \(\binom{n_i}{m_i}\) 相乘后对 \(p\) 取模。
组合数 mod 2：当 \(p = 2\) 时，\(\binom{n}{m} \equiv \binom{\lfloor n / 2 \rfloor}{\lfloor m / 2 \rfloor}\binom{n \bmod 2}{m \bmod 2} \pmod 2\)。
故 \(\binom{n}{m} \bmod 2 = 1\)，当且仅当在二进制下 \(m\) 是 \(n\) 的子集。
组合数常见求法
求法 1：递推
适用情况：对模数不做要求，可以承受平方时间的预处理（\(n \leq 5000\)）。
时间复杂度：预处理 \(\mathcal{O}(n^2)\)，查询 \(\mathcal{O}(1)\)。
初值 \(\binom{n}{0} = \binom{n}{n} = 1\)，递推式 \(\binom{n}{m} = \binom{n - 1}{m} + \binom{n - 1}{m - 1}\)。
for (int i = 0; i <= n; i ++) {
C[i][0] = C[i][i] = 1;
for (int j = 1; j < i; j ++) {
C[i][j] = (C[i - 1][j] + C[i - 1][j - 1]) % mod;
}
}
求法 2：预处理阶乘与阶乘逆元
适用情况：要求模数为质数，可以承受线性时间的预处理（\(n \leq 10^6\)）。
时间复杂度：预处理 \(\mathcal{O}(n)\)，查询 \(\mathcal{O}(1)\)。
将组合数使用阶乘形式表达
考虑预处理出阶乘与阶乘逆元。
设 \(\mathrm{fact}_i = i!\)，则有初值 \(\mathrm{fact}_0 = 1\)，递推式 \(\mathrm{fact}_i =\mathrm{fact}_{i - 1} \times i\)。
设 \(\mathrm{facv} = \frac{1}{i!}\)，则有初值 \(\mathrm{facv}_n = \mathrm{fact}_n^{-1}\)，递推式 \(\mathrm{facv}_i = \mathrm{facv}_{i + 1} \times(i + 1)\)。
于是先正着递推一遍得到 \(\mathrm{fact}_{0, \cdots, n}\)，然后通过一次求逆元得到 \(\mathrm{facv}_n\)，最后再反着递推一遍得到 \(\mathrm{facv}_{n, \cdots, 0}\)。
0x33 组合数（预处理阶乘与阶乘逆元）.cpp：
/* 组合数（预处理阶乘与阶乘逆元） */
struct BinomCoef {
std::vector<int> fact, facv;
BinomCoef() {}
BinomCoef(int n) {
init(n);
}
void init(const int &n) {
fact.resize(n + 1), facv.resize(n + 1);
fact[0] = 1;
for (int i = 1; i <= n; i ++) {
fact[i] = 1ll * fact[i - 1] * i % mod;
}
facv[n] = qpow(fact[n], mod - 2, mod);
for (int i = n - 1; i >= 0; i --) {
facv[i] = 1ll * facv[i + 1] * (i + 1) % mod;
}
}
int binom(int n, int m) {
if (n < m || m < 0) {
return 0;
}
return 1ll * facv[m] * facv[n - m] % mod * fact[n] % mod;
}
} bc;
求法 3：暴力统计
适用范围：要求模数为质数，\(m\) 的范围可以承受线性枚举。
时间复杂度：无预处理，查询 \(\mathcal{O}(m)\)。
将组合数表达为
注意到分子分母均只有 \(m\) 项，当 \(n\) 很大但 \(m\) 很小时，可以暴力枚举。
具体实现时，可以先分别求出分子与分母，然后再对分母求一次逆元。
求法 4：Lucas 定理
适用范围：要求模数为质数，模数以内的组合数可以预处理。
时间复杂度：预处理 \(\mathcal{O}(p)\)，查询 \(\mathcal{O}(\log_{p} n)\)。
由 Lucas 定理，得
先预处理出模数以内的组合数。然后不断地使用 Lucas 定理，直到 \(n, m\) 递降至模数以内。
求法 5：exLucas 定理
待填坑。
求法 6：分段打表
适用范围：要求模数为固定质数，代码长度限制不小，可以承受段长 \(\mathcal{O}(W)\) 时间的查询。
（正常数据范围 \(n \leq 10^9\)，配合 Lucas 定理可以做到 \(\mathrm{mod} \leq 10^9\) 且 \(n \leq 10^{18}\)）
时间复杂度：打表预处理，查询 \(\mathcal{O}(W)\)。
仍然将组合数使用阶乘形式表达，每次查询相当于要得到三个阶乘的值。
考虑分段打表。记段长为 \(W\)（取 \(10^5 \sim 10^6\) 相对合适），先通过打表预处理出所有 \(W\) 的倍数的阶乘。
每次查询 \(n!\) 的时候，就找到最近的 \(W\) 的倍数，根据该答案继续向上乘。
邪修：快速阶乘算法
适用范围：要求模数为不固定质数（\(n \leq 10^9\)）。
时间复杂度：无预处理，查询 \(\mathcal{O}(\sqrt{n} \log n)\)。
直接使用快速阶乘算法！得到三个阶乘的值。
求法 7：阶乘分解质因数
适用范围 1：对模数不做要求，可以承受线性时间的查询（\(n \leq 10^6\)）。
时间复杂度：预处理 \(\mathcal{O}(n)\)，查询 \(\mathcal{O}(n)\)。
适用范围 2：高精度组合数（\(n \leq 5000\)）。
仍然将组合数适用阶乘形式表达，我们考虑 \(1 \sim n\) 中的每个质数出现的次数。
首先先将 \(1 \sim n\) 中所有的质数提前筛出来，大约有 \(\frac{n}{\ln n}\) 个质数。对于质数 \(p\)，在阶乘 \(n!\) 中 \(p\) 的出现次数为
于是将三个阶乘都进行质因数分解，得出每个质数的出现次数，然后将所有质数的贡献相乘即可。
时间复杂度 \(\mathcal{O}(n)\)（高精度组合数的复杂度另算）。
斯特林数
第一类斯特林数
第二类斯特林数
卡特兰数
卡特兰数：\(\mathrm{Cat}_n\)，从原点 \((0, 0)\) 出发，每次只能向 \(x\) 轴或 \(y\) 轴正方向走，且不能越过直线 \(y = x\)，到达 \((n, n)\) 的方案数。
下标从 \(0\) 开始的若干项：
卡特兰数通项公式：
简单推导：越过直线 \(y = x\)，相当于是触碰直线 \(y = x + 1\)。
将触碰 \(y = x + 1\) 的第一个点以后的路径，关于 \(y = x + 1\) 对称，就可以得到一条从 \((0, 0)\) 到 \((n - 1, n + 1)\) 的路径。
于是「\((0, 0) \to (n, n)\) 且越过 \(y = x\) 的路径」与「\((0, 0) \to (n - 1, n + 1)\) 的路径」是一一对应的，方案数均为 \(\binom{2n}{n - 1}\)。
正难则反，不越过直线 \(y = x\) 的方案数，即为总方案数 \(\binom{2n}{n}\) 减去越过直线的方案数 \(\binom{2n}{n - 1}\)。
卡特兰数递推公式：
简单推导：枚举第一次触碰直线 \(y = x\) 的位置 \((i, i)\)，其中 \(i \geq 1\)。
由于在 \((i, i)\) 之前不能触碰直线 \(y = x\)，相当于是 \((1, 0) \to (i, i - 1)\) 不越过直线 \(y = x - 1\)，于是这部分的方案数为 \(\mathrm{Cat}_{i - 1}\mathrm{Cat}_{n - i}\)。
卡特兰数递推公式 2：
常见卡特兰数问题：
- 从原点 \((0, 0)\) 出发，每次只能向 \(x\) 轴或 \(y\) 轴正方向走，且不能越过直线 \(y = x\)，到达 \((n, n)\) 的方案数 \(\mathrm{Cat}_n\)。
- 以 \(1, 2, \cdots, n\) 作为 \(n\) 个元素的进栈序列时，出栈序列的方案数 \(\mathrm{Cat}_n\)。
- 以 \(n\) 对括号组成的合法括号序列的方案数 \(\mathrm{Cat}_n\)。
- 以 \(n\) 个节点构造无标号、左右儿子区分的二叉树个数 \(\mathrm{Cat}_n\)。
- 以 \(n\) 个节点构造无标号、儿子区分的有根树个数 \(\mathrm{Cat}_{n - 1}\)。
- 以圆周上的 \(2n\) 个点为端点，连互不相交的 \(n\) 条弦的方案数 \(\mathrm{Cat}_n\)。
- 将凸 \(n\) 边形用其中 \(n - 3\) 条对角线分割为互不重叠的 \(n - 2\) 个三角形的方案数 \(\mathrm{Cat}_{n - 2}\)。
- 将 \(1 \sim 2n\) 中的每个数，不重不漏地填入一个 \(2 \times n\) 的矩阵，每个元素大于其左方、上方元素的方案数 \(\mathrm{Cat}_n\)。
贝尔数
贝尔数：\(B_n\)，表示将 \(n\) 个两两不同的元素，划分成若干个互不区分的非空子集的方案数。
下标从 \(0\) 开始的若干项：
贝尔数递推公式：
简单推导：枚举与元素 \(n\) 在同一集合的元素个数 \(i\)，共有 \(\binom{n - 1}{i}\) 种选法，剩下的元素继续划分，于是这部分的方案数为 \(\binom{n - 1}{i}B_{n - 1 - i}\)。
贝尔数级别的暴搜：当我们需要将 \(n\) 个两两不同的元素，划分成若干个互不区分的非空子集时。若 \(n\) 较小可以直接暴搜，复杂度为 \(B_n\) 级别。
设 \(s_1, \cdots, s_{t}\) 表示当前时刻已经存在的集合，每次加入一个新元素 \(x\) 时有两种转移
- 选择一个已经存在的集合 \(s_i\)，将 \(x\) 加入 \(s_i\)。
- 新开一个集合 \(s_{t + 1}\)，将 \(x\) 加入 \(s_{t + 1}\)。
贝尔数与第二类斯特林数的关系：
因为第二类斯特林数 \({n \brace i}\) 表示将 \(n\) 个两两不同的元素，划分成 \(i\) 个互不区分的非空子集的方案数。
斐波那契数列
斐波那契数列：\(\mathrm{Fib}_n\)。
斐波那契数列递推公式：
斐波那契数列递推公式（矩阵形式）：
即
该转移矩阵可以求逆，从而实现逆向转移。
斐波那契数列通项公式：
当 \(5\) 是模 \(p\) 意义下的二次剩余时（一个常见模数是 \(10^9 + 9\)），可以直接暴力预先找出 \(\sqrt{5}\)，然后使用通项公式计算。
否则，考虑使用扩展数域进行计算，定义扩展数域 \(\mathbb{F}_p(\sqrt{5})\) 为
可以验证，\(\mathbb{F}_p(\sqrt{5})\) 关于加、减、乘、除封闭。记 \((a, b)\) 表示 \(a + b\sqrt{5}\)，则有
广义斐波那契数列：将递推公式反过来写，有 \(\mathrm{Fib}_n = \mathrm{Fib}_{n + 2} - \mathrm{Fib}_{n + 1}\)。
考虑将斐波那契数列扩展到负数域，经过推导可得
类斐波那契数列：已知数列 \(G_0 = a, G_1 = b\)，\(G_n = G_{n - 1} + G_{n - 2}\)，则有
类斐波那契数列相加，仍然是类斐波那契数列。只需将相应的 \(a, b\) 相加即可。
斐波那契数列常见恒等式
斐波那契数列求和：
斐波那契数列平方和：
卡西尼性质：
简单推导：\(\mathrm{Fib}_{n - 1}\mathrm{Fib}_{n + 1} - \mathrm{Fib}_n^2 = \det\begin{pmatrix} \mathrm{Fib}_{n + 1} & \mathrm{Fib}_n \\\mathrm{Fib}_n & \mathrm{Fib}_{n - 1} \end{pmatrix} = \det\begin{pmatrix} 1 & 1 \\1 & 0 \end{pmatrix}^n = (-1)^n\)。
附加性质：
简单推导：可以使用 "类斐波那契数列" 的知识解释。
记 \(G_0 = \mathrm{Fib}_n, G_1 = \mathrm{Fib}_{n + 1}\)，\(G_n = G_{n - 1} + G_{n - 2}\)，则有 \(\mathrm{Fib}_{n + m} = G_m\)。
斐波那契数列相邻项互质性质：
斐波那契数列 gcd 性质：
简单推导：不妨设 \(n < m\)，则有
由上文提到的 \(\mathrm{Fib}_{n + 1}\) 与 \(\mathrm{Fib}_n\) 互质。则有 \(\gcd(\mathrm{Fib}_n, \mathrm{Fib}_m) = \gcd(\mathrm{Fib}_n, \mathrm{Fib}_{m - n})\)。该式子与辗转相除法相似，得证。
斐波那契数列整除性质：
斐波那契数列模意义下的循环节
在模意义下，二元组 \((\mathrm{Fib}_n \bmod p, \mathrm{Fib}_{n + 1} \bmod p)\) 仅有 \(p^2\) 种取值。因此，当 \(n\) 足够大时，必定会出现循环节。
根据递推式，有 \(\mathrm{Fib}_n = \mathrm{Fib}_{n + 2} - \mathrm{Fib}_{n + 1}\)，也就是说 \(\mathrm{Fib}_{n + 1}, \mathrm{Fib}_{n + 2}\) 的前面必为 \(\mathrm{Fib}_n\)。
根据这两个事实我们可以得知，斐波那契数列在模意义下必有循环节，并且一定是纯循环的（循环节的起点为 \(\mathrm{Fib}_0, \mathrm{Fib}_1\)）。
最小循环节（皮萨诺周期）：定义斐波那契数列模 \(p\) 意义下的最小循环节 \(\pi(p)\) 为最小的正整数 \(m\)，使得 \(\mathrm{Fib}_k \equiv 0 \pmod p\) 且 \(\mathrm{Fib}_{k + 1} \equiv 1 \pmod p\)。
注意：周期一定是最小循环节的倍数，即
很多时候我们只需求出周期即可。
最小循环节的上界：斐波那契数列 \(\operatorname{mod} m\) 的最小正周期不超过 \(6m\)，即 \(\pi(m) \leq 6m\)。
求斐波那契数列循环节 - 结论法：
- 质数情景：
- \(\pi(2) = 3\)，\(\pi(5) = 20\)。
- 当 \(p\) 为奇质数且 \(p \equiv \pm 1 \pmod 5\) 时，\(\pi(p) \mid p - 1\)。
- 当 \(p\) 为奇质数且 \(p \equiv \pm 2 \pmod 5\) 时，\(\pi(p) \mid 2p + 2\)。
- 一般情景：
- 对于质数 \(p\) 以及任意正整数 \(k\)，有 \(\pi(p^k) \mid p^{k - 1}\pi(p)\)。
- 对于 \(m = \prod p_i^{c_i}\)，有 \(\pi(m) = \mathrm{lcm}\{ \pi(p_i^{c_i}) \}\)，同时有 \(\pi(m) \mid \mathrm{lcm}\{ p_i^{c_i - 1}\pi(p_i) \}\)。
在一般应用中，通常只需要求出一个循环节即可。
对于质数 \(p\)，记
然后将模数 \(m\) 质因数分解成 \(\prod p_i^{c_i}\)，使用 \(\mathrm{lcm}\{ p_i^{c_i - 1}g(p_i) \}\) 作为循环节即可（注意这并不是最小循环节），这个值在多数情况下已经够用了。
求斐波那契数列循环节 - 矩阵 BSGS：
适用场景：最小循环节的上界确定，转移矩阵可逆。
考虑斐波那契数列的矩阵表示
记转移矩阵 \(A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}\)，那么循环节即为满足 \(A^m \equiv I \pmod p\) 的正整数 \(m\)。
在转移矩阵可逆的情况下，套用 BSGS 的流程：取 \(S\) 表示上界的根号，设 \(m = iS - j\)（其中 \(i \geq 1\) 且 \(0 \leq j < S\)），则有
预处理 \(0 \leq j < S\) 的 \(A^j\) 存入哈希表，每次在哈希表中查询 \(A^{iS}\) 即可。
时间复杂度 \(\mathcal{O}(\sqrt{p})\)。
求斐波那契数列循环节 - 随机化（生日悖论）：
若存在两个数 \(0 \leq i, j \leq 6p\)，满足 \((\mathrm{Fib}_i, \mathrm{Fib}_{i + 1}) \equiv (\mathrm{Fib}_j, \mathrm{Fib}_{j + 1}) \pmod p\)，则 \(|i - j|\) 为一个循环节（注意这并不是最小循环节）。
每次新随机一个数 \(x\)，将二元组 \((\mathrm{Fib}_x, \mathrm{Fib}_{x + 1})\) 存入哈希表中，并查询哈希表中有没有与之相同的 \((\mathrm{Fib}_y, \mathrm{Fib}_{y + 1})\)。
发现这与生日悖论的形式很像，于是时间复杂度正确。
生日悖论：在 \([1, n]\) 中随机选择正整数，选出两个相同数的期望轮数为 \(\mathcal{O}(\sqrt{n})\) 级别（近似 \(\sqrt{\frac{\pi n}{2}}\)）。
实际取数过程可能需要在 \([0, 12p]\) 中随机以保证期望次数。求 \(\mathrm{Fib}_{i}\) 需要用到 "光速幂" 的技巧。
0x34 数学知识 初等数论
质数
质数：若一个正整数 \(n\) 无法被除了 \(1\) 和它自身以外的任何正整数整除，则称该数为质数（或素数）。否则称该数为合数。
质数分布函数：\(\pi(n)\) 表示 \(1 \sim n\) 中的质数个数。当 \(n\) 足够大时，\(\pi(n) \sim \frac{n}{\ln n}\)。
唯一分解定理：任意一个大于 \(1\) 的正整数 \(n\)，都可以唯一分解成有限个质数的乘积，记作
其中 \(p_i\) 均为质数且满足 \(p_1 < p_2 < \cdots < p_t\)，\(c_i\) 均为正整数。
质因数的根号分治：一个正整数 \(n\)，大于 \(\sqrt{n}\) 的质因子最多只有 \(1\) 个。
质数判定
若一个正整数 \(n\) 为合数，则至少存在一个正整数 \(d\)，其中 \(2 \le d \leq \sqrt{n}\)。
于是我们只需枚举 \(2 \sim \sqrt{n}\) 之间的所有整数，依次检查是否能够整除 \(n\) 即可。
时间复杂度 \(\mathcal{O}(\sqrt{n})\)。
0x34 质数判定.cpp：
/* 质数判定 */
bool isPrime(int n) {
// assert(n > 1);
for (int i = 2; 1ll * i * i <= n; i ++) {
if (n % i == 0) {
return 0;
}
}
return 1;
}
质因数分解
尝试枚举 \(2 \sim \sqrt{n}\) 之间的所有整数 \(d\)。若 \(d\) 能够整除 \(n\)，则从 \(n\) 之中删除 \(d\) 的所有因子，并且累计 \(d\) 的个数。
由于 \(n\) 的合数因子，必定在枚举该因子之前就从 \(n\) 中删除掉了。所以上述过程中枚举到的因数 \(d\) 必定是质因数。
若最后的 \(n\) 没有分解完毕（仍有 \(n > 1\)），则剩下的 \(n\) 也为质因数。
时间复杂度 \(\mathcal{O}(\sqrt{n})\)。
0x34 质因数分解.cpp：
/* 质因数分解 */
auto PrimeFactor(int n) {
std::vector< std::pair<int, int> > seq;
for (int i = 2; 1ll * i * i <= n; i ++) {
if (n % i == 0) {
int c = 0;
while (n % i == 0) {
n /= i, c ++;
}
seq.push_back({i, c});
}
}
if (n > 1) {
seq.push_back({n, 1});
}
return seq;
}
阶乘的质因数分解
对于阶乘 \(n! = 1 \times 2 \times \cdots \times n\)，由于每一项均为 \([1, n]\) 中的正整数，故分解出的质因数也均为 \([1, n]\) 中的质数。
筛出 \(1 \sim n\) 中的所有质数，对于每个质数 \(p\)，考虑 \(n!\) 中一共包含了多少个质因子 \(p\)。
共有大约 \(\frac{n}{\ln n}\) 个质数，每个质数 \(p\) 的求解需要 \(\log_p n\) 的开销，故时间复杂度 \(\mathcal{O}(n)\)。
阶乘 \(n!\) 中包含质数 \(p\) 的次数：
简单推导：设 \(c_i\) 表示 \(i\) 包含的质因子 \(p\) 个数，则阶乘 \(n!\) 包含 \(p\) 的个数为 \(\sum_{j = 1}^{\infty} \sum_{i = 1}^n j[c_i = j]\)，进一步得 \(\sum_{j = 1}^{\infty}\sum_{i = 1}^n [c_i \geq j]\)，不难发现其中 \(\sum_{i = 1}^n [c_i \geq j] = \left\lfloor \frac{n}{p^j} \right\rfloor\)。
阶乘 \(n!\) 中包含合数 \(m\) 的次数：将 \(m\) 唯一分解成 \(m = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\) 的形式。
记 \(v_p(n) = \sum_{i = 1}^{\lfloor \log_p n \rfloor} \left\lfloor \frac{n}{p^i} \right\rfloor\)，答案为
埃氏筛
注意到，任意整数 \(x(x > 1)\) 的倍数 \(2x, 3x, \cdots\) 都不是质数。
从 \(2\) 到 \(n\) 依次扫描每一个数 \(x\)，若 \(x\) 尚未被标记，则说明 \(x\) 为质数，并且将其倍数 \(2x, 3x, \cdots, \left\lfloor \frac{n}{x} \right\rfloor x\) 均标记为合数。
时间复杂度 \(\mathcal{O}(n \log \log n)\)。
在开 O2 优化的情况下，std::bitset 连续读写效率较高，可以加以优化埃氏筛。
埃氏筛的时间复杂度证明：相当于是证明 \(1 \sim n\) 中所有质数的倒数之和为 \(\mathcal{O}(\log \log n)\) 级别。
由质数分布函数，可以推测第 \(n\) 个质数的大小为 \(\Theta(n \log n)\)，于是有
\[\begin{aligned} \sum_{k = 1}^{\pi(n)} \frac{1}{p_k} & = \mathcal{O}\left( \sum_{k = 2}^{\pi(n)} \frac{1}{k \log k} \right) \\ & = \mathcal{O}\left(\int_2^{\pi(n)} \frac{1}{x \log x} \mathrm{dx} \right) \\ & = \mathcal{O}(\log \log \pi(n)) \\ & = \mathcal{O}(\log \log n) \end{aligned} \]
0x34 埃氏筛.cpp：
std::bitset<MaxV> vis;
/* 埃氏筛 */
void sieve(const int &n) {
for (int i = 2; i <= n; i ++) {
if (vis[i]) {
continue;
}
for (int j = 2; j <= n / i; j ++) {
vis[i * j] = 1;
}
}
}
线性筛
为了避免埃氏筛重复标记合数。在线性筛中，我们希望确定一个合数的方式能够唯一。
从 \(2\) 到 \(n\) 依次扫描每一个数 \(x\)，在标记合数时，我们只在 \(x\) 的基础上乘上一个质因子 \(p\)，并且保证 \(p\) 是 \(xp\) 的最小质因子。这样做，每个合数只会被自己的最小质因子标记一次。并且顺带还能筛出每个数的最小质因子。
时间复杂度 \(\mathcal{O}(n)\)。
由于线性筛每次只扩展一个最小质因子的结构，可以额外处理很多信息（例如筛出积性函数）。
0x34 线性筛.cpp：
int primeCount, prime[MaxV], fac[MaxV];
/* 线性筛 */
void sieve(const int &n) {
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
}
}
}
Miller-Rabin 素性测试
回顾一些内容：
费马小定理：设 \(p\) 是素数，对于任意整数 \(a\)（\(a\) 不能是 \(p\) 的倍数），满足 \(a^{p - 1} \equiv 1 \pmod p\)。
Fermat 伪素数：若满足 \(a^n \equiv 1 \pmod n\)，但 \(n\) 不是素数，则称 \(n\) 为以 \(a\) 为底的 Fermat 伪素数。
对于任何固定的基底 \(a\)，都存在无穷多个的 Fermat 伪素数。
Fermat 素性测试：不断地选取在 \([2, n - 1]\) 中的基底 \(a\)，并检验是否每次都满足 \(a^{n - 1} \equiv 1\pmod n\)。
即使检查了所有与 \(n\) 互质的基底 \(a\)，依然无法保证 \(n\) 是素数。这样的数 \(n\) 称为 Carmichael 数，有无穷多个。
二次探测定理：设 \(p\) 是奇素数，则同余方程 \(x^2 \equiv 1 \pmod p\) 的解为 \(x \equiv 1 \pmod p\) 或 \(x \equiv p - 1\pmod p\)。
Miller-Rabin 素性测试：
Miller-Rabin 素性测试，结合了费马小定理和二次探测定理，是一种概率性素性测试。
实际上，没有已知的数字通过了 Miller-Rabin 素性测试但实际上是合数。可以放心使用。
先将 \(a^{n - 1} \equiv 1 \pmod n\) 中的指数 \(n - 1\)，分解成 \(n - 1 = u \times 2^t\) 的形式。
每次随机一个基底 \(a\)，进行如下的判断：
- 先使用快速幂求出 \(v = a^u \bmod n\)。如果 \(v = 1\) 或 \(v = p - 1\)，则该轮测试通过。
- 然后对这个值进行 \(t - 1\) 次平方操作，途中如果该值等于 \(p - 1\)，则该轮测试通过。
- 否则该轮测试不通过。
在不考虑乘法的开销时，对 \(n\) 进行 \(k\) 轮测试的时间复杂度为 \(\mathcal{O}(k \log n)\)。
可以证明，奇合数 \(n(n > 9)\) 通过随机选取的一个基底的概率不超过 \(\frac{1}{4}\)。因此随机选取 \(k\) 个基底后，仍然误判的概率不超过 \(\frac{1}{4^k}\)。
对于 \([2, 2^{64})\) 范围内的数，工业上通常选用 \(\{ 2, 325, 9375, 28178, 450775, 9780504, 1795265022 \}\) 作为检验的基底。为了方便也可以选择前 \(12\) 个质数作为基底。
Pollard-Rho 质因数分解
Pollard-Rho：用于找出合数 \(n\) 的某一个非平凡因数 \(d\)（非平凡因数：非 \(1\) 且非 \(n\) 的因数）。
若一个正整数 \(n\) 为合数，那么 \(n\) 的最小质因子 \(p\) 满足 \(p \leq \sqrt{n}\)。
如果随机生成 \(\mathcal{O}(\sqrt{p})\) 个正整数，由生日悖论，可能存在两个正整数 \(a, b\) 满足 \(a \equiv b \pmod p\)，此时 \(|a - b|\) 是 \(p\) 的倍数，于是查询 \(\gcd(|a - b|, n)\) 即可得到一个 \(n\) 的一个因子。
生日悖论：在 \([1, n]\) 中随机选择正整数，选出两个相同数的期望轮数为 \(\mathcal{O}(\sqrt{n})\) 级别（近似 \(\sqrt{\frac{\pi n}{2}}\)）。
但这样做，枚举数对的时间复杂度还是太高。
记映射 \(f(x) = (x^2 + c) \bmod n\)，其中 \(c\) 是一个随机选取的常数（\(c \neq 0, -2\)）。Pollard 构造了这样的一个伪随机序列 \(\{x_i\}\)，满足 \(x_i = f(x_{i - 1})\)。
考虑序列 \(\{ x_i \bmod p \}\)，由于 \(x_i\) 仅仅与 \(x_{i - 1}\) 有关，且 \(x_i \bmod p \in [0, p)\)，所以 \(\{ x_i \bmod p \}\) 是一个混循环序列。
建图后形似 \(\rho\) 形（rho），由生日悖论可得尾长与环长均为 \(\mathcal{O}(\sqrt{p})\) 级别。
我们发现，若 \(x \equiv y \pmod p\)，必有 \(f(x) \equiv f(y) \pmod p\)。由此可得在模 \(p\) 的环上，只要某个相距 \(d\) 步的点对满足条件，则所有相距 \(d\) 步的点对均满足条件。所以我们每次最好要检查不同的距离。
Brent 倍增迭代：初始有两个点 \(a, b\)，在第 \(k\) 轮，\(a\) 原地不动，\(b\) 不断地向前共计 \(2^k\) 步，最后再让 \(a\) 直接跳到 \(b\) 的位置。\(b\) 的每一步移动，我们都计算 \(d = \gcd(|a - b|, n)\)，当第一次出现 \(d > 1\) 时
- 若 \(1 < d < n\)，则我们已经找到了一个非平凡因数，直接返回 \(d\)。
- 若 \(d = n\)，则说明我们已经完整地跑完了一个环，继续在环上跑也没有意义。
此时我们应该换个映射，或者换一个初始值继续做。
进一步，我们想要降低求解 \(\gcd\) 的次数。注意到，如果 \(\gcd(a, n) > 1\)，则对于任意正整数 \(b\) 均有 \(\gcd(ab \bmod n, n) = \gcd(ab, n) > 1\)。
也就是说，如果 \(\gcd\left(\prod |x_i - x_j| \bmod n, n\right) > 1\)，那么其中必有一对 \(x_i, x_j\) 满足 \(\gcd(|x_i - x_j|, n) > 1\)。
于是我们每 \(k\) 对 \(a, b\) 就计算一次 \(\gcd\)。时间复杂度 \(\mathcal{O}(\sqrt{p} + k^{-1}\sqrt{p} \log n)\)，当 \(k\) 与 \(\log n\) 大致相等时，可以得到 \(\mathcal{O}(\sqrt{p})\) 的理论复杂度。具体实现中，一般取 \(k = 128\)。
Pollard-Rho 质因数分解：如果要对 \(n\) 进行质因数分解，先使用 Miller-Rabin 判断 \(n\) 是否是质数，若 \(n\) 为合数，则使用 Pollard-Rho 找出 \(n\) 的一个非平凡因数 \(d\)，得到两个规模分别为 \(d\) 和 \(\frac{n}{d}\) 的两个子问题，递归进行求解。
0x34 素性测试与质因数分解（Miller-Rabin + Pollard-Rho）.cpp：
/* Miller-Rabin + Pollard-Rho */
template <class T, auto qmul>
struct MRPR {
constexpr T qpow(T a, T b, T p) {
T ans = 1;
for (; b; b >>= 1) {
if (b & 1) ans = qmul(ans, a, p);
a = qmul(a, a, p);
}
return ans;
}
std::mt19937_64 mtrand{std::random_device{}()};
constexpr T rand(T l, T r) {
std::uniform_int_distribution<T> range(l, r);
return range(mtrand);
}
/* Miller-Rabin 素性测试 */
bool isPrime(T n) {
// assert(n > 1);
if (n % 2 == 0) return n == 2;
if (n % 3 == 0) return n == 3;
T t = __builtin_ctzll(n - 1);
T u = (n - 1) >> t;
int test = 6;
while (test --) {
T a = rand(2, n - 1);
T v = qpow(a, u, n);
if (v == 1 || v == n - 1) {
continue;
}
for (int i = 1; i < t; i ++) {
v = qmul(v, v, n);
if (v == n - 1) {
break;
}
}
if (v != n - 1) {
return 0;
}
}
return 1;
}
T walk(T n) {
T c = rand(1, n - 1);
T x = rand(1, n - 1), y, v, d;
for (int k = 1;; k <<= 1) {
y = x, v = 1;
for (int i = 1; i <= k; i ++) {
x = qmul(x, x, n) + c; if (x >= n) x -= n;
v = qmul(v, x - y, n);
if (i % 127 == 0 || i == k) {
if (v < 0) v += n;
d = std::gcd(v, n);
if (d > 1) return d;
}
}
}
}
/* Pollard-Rho 质因数分解 */
auto PrimeFactor(T n) {
std::map<T, int> buc;
std::function<void(T)> solve = [&] (T n) {
if (n <= 10000) {
for (int i = 2; i * i <= n; i ++) {
while (n % i == 0) {
n /= i;
buc[i] ++;
}
}
if (n > 1) {
buc[n] ++;
}
return;
}
if (isPrime(n)) {
buc[n] ++;
return;
}
T d = 0;
while (d = walk(n), d == n);
solve(d), solve(n / d);
};
solve(n);
std::vector< std::pair<T, int> > seq;
for (auto dat : buc) {
seq.push_back(dat);
}
return seq;
}
};
// MRPR<int, [] (int a, int b, int p) { return 1ll * a * b % p; }> mrpr;
// MRPR<s64, [] (s64 a, s64 b, s64 p) { return static_cast<__int128>(a) * b % p; }> mrpr;
约数
整除：设 \(d, n \in \Z\) 且 \(d \neq 0\)，如果 \(\exist q \in \Z\) 使得 \(qd = n\)，那么称 \(d\) 整除 \(n\)（记作 \(d \mid n\)）。
整除的简单性质：
- \(a \mid b \land b \mid c \Longrightarrow a\mid c\)。
- \(a \mid b \land a \mid c \Longrightarrow \forall x, y \in \Z, a \mid (xb + yc)\)。
- 设 \(m \neq 0\)，则 \(a \mid b \Longrightarrow ma \mid mb\)。
- 设 \(\gcd(a, m) = 1\)，则 \(a \mid mb \Longrightarrow a\mid b\)。
单个数的约数集合
若 \(d \geq \sqrt{n}\) 是 \(n\) 的约数，则 \(\frac{n}{d} \leq \sqrt{n}\) 也一定是 \(n\) 的约数。
也就是说，约数 \(d\) 和 \(\frac{n}{d}\) 一定是成对出现的（除了完全平方数，\(\sqrt{n}\) 单独出现），并且其中一个 \(\leq \sqrt{n}\)，另一个 \(> \sqrt{n}\)。
于是我们只需枚举 \(2 \sim \sqrt{n}\) 之间的所有整数 \(d\)。若 \(d \mid n\)，则将 \(d\) 和 \(\frac{n}{d}\) 一起加入 \(n\) 的约数集合。
时间复杂度 \(\mathcal{O}(\sqrt{n})\)。
由上述分析，可知 \(n\) 的约数个数必定不超过 \(2\sqrt{n}\)。实际上很难达到该上界，在算法竞赛中可以简单认为 \(10^6 \sim 10^{18}\) 范围内约数个数为 \(\mathcal{O}(\sqrt[3]{n})\) 级别。
如果我们想要得到排序后的结果。注意到枚举的 \(d\) 是升序的，\(\frac{n}{d}\) 是降序的，并且 \(\frac{n}{d}\) 的部分一定比 \(d\) 的部分数值更大。
于是我们只需将 \(\frac{n}{d}\) 的部分反转后插入 \(d\) 的部分的末尾即可。
0x34 单个数的约数集合.cpp：
/* 单个数的约数集合 */
auto Factor(int n) {
std::vector<int> seq;
for (int i = 1; 1ll * i * i <= n; i ++) {
if (n % i == 0) {
seq.push_back(i);
if (i * i < n) {
seq.push_back(n / i);
}
}
}
return seq;
}
0x34 单个数的约数集合（排序）.cpp：
/* 单个数的约数集合（排序） */
auto Factor(int n) {
std::vector<int> lhs, rhs;
for (int i = 1; 1ll * i * i <= n; i ++) {
if (n % i == 0) {
lhs.push_back(i);
if (i * i < n) {
rhs.push_back(n / i);
}
}
}
for (int i = (int)rhs.size() - 1; i >= 0; i --) {
lhs.push_back(rhs[i]);
}
return lhs;
}
1~n 的约数集合
反过来考虑 \(i\) 可以作为哪些数的约数，显然 \(i\) 可以作为 \(i, 2i, \cdots, \left\lfloor \frac{n}{i} \right\rfloor i\) 的约数。
时间复杂度 \(\sum_{i = 1}^n \frac{n}{i} = \mathcal{O}(n \log n)\)。
由上述分析，可知 \(1 \sim n\) 的约数个数总和为 \(\mathcal{O}(n \log n)\) 级别。
0x34 1~n 的约数集合.cpp：
std::vector<int> factor[MaxV];
/* 1~n 的约数集合 */
void FactorInit(const int &n) {
for (int i = 1; i <= n; i ++) {
for (int j = i; j <= n; j += i) {
factor[j].push_back(i);
}
}
}
最大公约数与最小公倍数
最大公约数：最大公共约数（Greatest Common Divisor，简写 GCD）。
最小公倍数：最小公共倍数（Least Common Multiple，简写 LCM）。
唯一分解角度下的 gcd 与 lcm：设 \(a = p_1^{c_1}p_2^{c_2}\cdots\)，\(b = p_1^{d_1}p_2^{d_2} \cdots\)，则有
在唯一分解角度下，gcd 相当于是对每个指数取 min，lcm 相当于是对每个指数取 max。
gcd 的简单性质：
- \(\gcd(a_1, \cdots, a_n) = \gcd(|a_1|, \cdots, |a_n|)\)。
- \(\gcd(a, 0) = \gcd(a, a) = |a|\)。
- \(\gcd(bq + r, b) = \gcd(r, b)\)。
- 设 \(\gcd(a, m) = 1\)，则 \(\gcd(a, mb) = \gcd(a, b)\)。
- \(\gcd(ma_1, \cdots, ma_n) = |m| \gcd(a_1, \cdots, a_n)\)。
lcm 的简单性质：
- \(\mathrm{lcm}(a_1, \cdots, a_n) = \mathrm{lcm}(|a_1|, \cdots, |a_n|)\)。
- \(\mathrm{lcm}(a, 1) = \mathrm{lcm}(a, a) = |a|\)。
- \(\mathrm{lcm}(ma_1, \cdots, ma_n) = |m| \operatorname{lcm}(a_1, \cdots, a_n)\)。
gcd 和 lcm 的常见恒等式
gcd-lcm 恒等式：
请注意！该式子扩展到三元以上几乎不成立！
事实上当 \(n \geq 3\) 时，\(\gcd(a_1, \cdots, a_n) \operatorname{lcm}(a_1, \cdots, a_n) = a_1\cdots a_n\) 当且仅当 \(a_1, \cdots, a_n\) 两两互质。
gcd 和 lcm 的结合律：
该式子扩展到任意多元仍然成立。
于是可以使用恒等式 \(\mathrm{lcm}(a, b) = \frac{ab}{\gcd(a, b)}\)，用来求解 \(\mathrm{lcm}(a, b, c)\)，以及 \(\mathrm{lcm}(a_1, \cdots, a_n)\)。
gcd 对乘法的分配律：
根据对乘法的分配律计算 \(\gcd(a, b_1\cdots b_n)\)。每次处理 \(b_i\) 时，计算 \(t = \gcd(a, b_i)\)，然后将 \(a\) 除以 \(t\)（相当于是将已计算的质因数贡献去掉），答案乘以 \(t\)。
gcd 的差分律：
通过差分律，我们可以维护差分数组的 gcd，以实现 "区间加" 与 "区间 gcd"。
斐波那契数列 gcd 恒等式：
幂差 gcd 恒等式：不妨设 \(a \geq b\)，当 \(\gcd(a, b) = 1\) 时，有
简单推导：
右边整除左边：可以证明 \(x - y \mid x^k - y^k\)。
\(k\) 次方差公式：\(x^k - y^k = (x - y)(x^{k - 1} + x^{k - 2}y + \cdots + xy^{k - 2} + y^{k - 1})\)。
对于 \(n, m\) 的公约数 \(d\) 来说，可以证得 \(a^d - b^d \mid a^n - b^n\) 与 \(a^d - b^d \mid a^m - b^m\)。
于是 \(a^d - b^d \mid \gcd(a^n - b^n, a^m - b^m)\)。所以
左边整除右边：
"降幂引理"：不妨设 \(r \geq s\)，当 \(\gcd(a, b) = 1\) 时，有
\[\gcd(a^r - b^r, a^s - b^s) \mid a^{r - s} - b^{r - s} \]证明：设 \(d\) 为 \(a^r - b^r, a^s - b^s\) 的公约数。
因为 \((a^r - b^r) - a^{r - s}(a^s - b^s) = b^s(a^{r - s} - b^{r - s})\)，可以证得 \(d \mid b^s(a^{r - s} - b^{r - s})\)；同理可得 \(d \mid a^s(a^{r - s} - b^{r - s})\)。
又由于 \(\gcd(a, b) = 1\)，于是 \(d \mid a^{r - s} - b^{r - s}\)，"降幂引理" 得证。
对 \(\gcd(a^n - b^n, a^m - b^m)\) 反复应用 "降幂引理"，过程等价于辗转相除法。所以
最大公约数常见求法
已有的轮子：std::gcd(a, b)。
求法 1：辗转相除法
（也称之为欧几里得算法）
对于 \(\forall a, b \in \N\) 且 \(b \neq 0\)，有
- 当 \(a < b\) 时，一次迭代后就会转化为 \(a \geq b\)。
- 当 \(a \geq b\) 时，由于 \(a \bmod b \leq \frac{a}{2}\)，这样的迭代不超过 \(\log_2 a\) 次。
时间复杂度 \(\mathcal{O}(\log \max(a, b))\)。
使用辗转相除法求解斐波那契数列相邻两项的 gcd，会让该算法达到最坏复杂度。
0x34 gcd（辗转相除法）.cpp：
/* 辗转相除法 */
int gcd(int a, int b) {
return b == 0 ? a : gcd(b, a % b);
}
求法 2：更相减损术
（也称之为 Stein 算法）
对于 \(\forall a, b \in \N\)，不妨设 \(a\geq b\)，有
- 当 \(a = b\) 时，\(\gcd(a, b) = a\)。
- "两个偶数"：\(\gcd(a, b) = 2\gcd(\frac{a}{2}, \frac{b}{2})\)。
- "一奇一偶"：不妨设 \(a\) 为奇数，\(\gcd(a, b) = \gcd(\frac{a}{2}, b)\)。
- "两个奇数"：不妨设 \(a \geq b\)，\(\gcd(a, b) = \gcd(a - b, b)\)。
当至少有一个偶数的时候，每次迭代至少会将 \(a, b\) 之一减半。而两个奇数相减会得到偶数。
时间复杂度 \(\mathcal{O}(\log \max(a, b))\)。
更相减损术常用于计算高精度 gcd（因为 "高精度取模" 复杂度较高，"高精度减法" 以及 "高精除低精" 复杂度较低）。
0x34 gcd（更相减损术）.cpp：
/* 更相减损术 */
int gcd(int a, int b) {
if (a == b) {
return a;
}
if (a == 0 || b == 0) {
return a ^ b;
}
int az = __builtin_ctz(a), bz = __builtin_ctz(b), z = az < bz ? az : bz, v = 0;
a >>= az, b >>= bz;
while (b) {
v = b - a;
a = a < b ? a : b;
b = v < 0 ? -v : v;
if (b) b >>= __builtin_ctz(b);
}
return a << z;
}
求法 3：递推
初值 \(\gcd(i, 0) = \gcd(0, i) = i\)，递推式 \(\gcd(i, j) = \gcd(i - j, j)\)。
时间复杂度：预处理 \(\mathcal{O}(n^2)\)，查询 \(\mathcal{O}(1)\)。
for (int i = 0; i <= n; i ++) {
g[i][0] = g[0][i] = i;
for (int j = 1; j <= i; j ++) {
g[i][j] = g[j][i] = g[i - j][j];
}
}
求法 4：基于值域预处理的快速 gcd
引理：对于任意正整数 \(n\)，一定可以被拆分成三个正整数 \(a, b, c(a \leq b \leq c)\) 的乘积，其中 \(a, b\leq \sqrt{n}\)，且 \(c \leq \sqrt{n}\) 与 \(c \in \mathbb{P}\) 至少有一个成立。
证明 1：假设 \(c > \sqrt{n}\) 且 \(c \notin \mathbb{P}\)。因为 \(c > \sqrt{n}\) 且 \(abc = n\)，所以 \(ab \leq \sqrt{n}\)。由于 \(c\) 不是质数，设 \(c = xy\)，显然 \(x, y\) 不可能同时 \(>\sqrt{n}\)，不妨设 \(x \leq y\)，则必有 \(x \leq \sqrt{n}\)。此时 \(n\) 可以被拆分成 \(x, ab, y\) 的乘积。不断对 \(y\) 进行相同的讨论。
由于数值不断递减，迭代到最后 \(c \leq \sqrt{n}\) 与 \(c \in \mathbb{P}\) 至少有一个成立。
证明 2：考虑归纳。当 \(n = 1\) 时显然成立。当 \(n > 1\) 时，设 \(p\) 为 \(n\) 的最小质因数，设 \(a', b', c'(a' \leq b' \leq c')\) 为 \(\frac{n}{p}\) 的一组合法拆分。
- 当 \(p \leq \sqrt[4]{n}\) 时。由于 \(a' \leq \sqrt[3]{\frac{n}{p}}\)，所以 \(pa' \leq \sqrt[3]{np^2} \leq \sqrt{n}\)。此时 \(pa', b', c'\) 为一组合法拆分。
- 当 \(p > \sqrt[4]{n}\) 时。假设 \(a' \geq p\)，则 \(n = pa'b'c' \geq p^4 > n\) 矛盾，故 \(a' < p\)。由于 \(p\) 为 \(n\) 的最小质因数，所以 \(a' = 1\)。此时 \(p, b', c'\) 为一组合法拆分。
预处理：通过线性筛得出 \(1 \sim n\) 所有数的最小质因数。然后根据证明 2，递推得出 \(1 \sim n\) 所有数的合法拆分，对于 \(i > 1\)，设 \(p\) 为 \(i\) 的最小质因数，将 \(\frac{i}{p}\) 的合法拆分中的最小数乘以 \(p\) 即可得到 \(i\) 的合法拆分。最后预处理出所有 \(0 \leq x, y \leq \sqrt{n}\) 的 \(\gcd(x, y)\)。
查询：当询问 \(\gcd(p, q)\) 时，根据 gcd 对乘法的分配律，依次考虑 \(q\) 的合法拆分 \(a, b, c\)，设当前考虑到的数为 \(x\)。
- 计算 \(t = \gcd(p, x)\)，具体地：
- 若 \(x\) 为合数（此时 \(x \leq \sqrt{n}\)），则 \(\gcd(p, x) = \gcd(x, p \bmod x)\)。使用预处理的 gcd 值回答即可。
- 若 \(x\) 为质数，则当且仅当 \(p\) 为 \(x\) 的倍数时 \(\gcd(p, x) = x\)，否则 \(\gcd(p, x) = 1\)。
- 然后将 \(p\) 除以 \(t\)，答案乘以 \(t\)。
时间复杂度：预处理 \(\mathcal{O}(n)\)，查询 \(\mathcal{O}(1)\)。
0x34 gcd（基于值域预处理的快速 gcd）.cpp：
/* O(n)-O(1) gcd */
struct flashGCD {
int b;
std::vector<int> prime, fac;
std::vector< std::array<int, 3> > sp;
std::vector< std::vector<int> > g;
flashGCD() {}
flashGCD(int n) {
init(n);
}
void init(int n) {
b = sqrt(n);
prime.clear(), fac.assign(n + 1, 0);
sp.assign(n + 1, {0, 0, 0});
g.assign(b + 1, std::vector<int>(b + 1, 0));
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime.push_back(i), fac[i] = i;
}
for (int p : prime) {
if (p > fac[i] || p > n / i) break;
fac[i * p] = p;
}
}
sp[1] = {1, 1, 1};
for (int i = 2; i <= n; i ++) {
sp[i] = sp[i / fac[i]];
sp[i][0] *= fac[i];
if (sp[i][0] > sp[i][1]) {
std::swap(sp[i][0], sp[i][1]);
}
if (sp[i][1] > sp[i][2]) {
std::swap(sp[i][1], sp[i][2]);
}
}
for (int i = 0; i <= b; i ++) {
g[i][0] = g[0][i] = i;
for (int j = 1; j <= i; j ++) {
g[i][j] = g[j][i] = g[i - j][j];
}
}
}
int gcd(int p, int q) {
if (p == 0 || q == 0) {
return p ^ q;
}
int ans = 1, t;
for (int x : sp[q]) {
if (x <= b) {
t = g[x][p % x];
} else {
t = p % x ? 1 : x;
}
p /= t, ans *= t;
}
return ans;
}
} fg;
类欧几里得算法
类欧几里得算法：用来解决一类形如 \(\left\lfloor \frac{ai + b}{c} \right\rfloor\) 结构的求和问题。
基本模型为
首先，可以转化为 \(0 \leq a, b < c\) 的情况
然后，令 \(m = \left\lfloor \frac{an + b}{c} \right\rfloor\)
结合这两部操作，二元组 \((a, c) \to (c, a \bmod c)\)，相当于对 \((a, c)\) 进行辗转相除法，故称之为类欧几里德算法。
时间复杂度 \(\mathcal{O}(\log \min(a, c))\)。
0x34 类欧几里德算法.cpp：
/* 类欧几里德算法（要求 a, b >= 0 且 c > 0） */
s64 euclid(s64 a, s64 b, s64 c, s64 n) {
if (a == 0) {
return (n + 1) * (b / c);
}
if (a >= c || b >= c) {
return euclid(a % c, b % c, c, n) +
(n * (n + 1) / 2) * (a / c) + (n + 1) * (b / c);
}
s64 m = (a * n + b) / c;
if (m == 0) {
return 0;
}
return n * m - euclid(c, c - b - 1, a, m - 1);
}
/* 类欧几里德算法 */
s64 euclid_full(s64 a, s64 b, s64 c, s64 n) {
if (c < 0) {
a = -a, b = -b, c = -c;
}
s64 res = 0;
if (a < 0) {
s64 t = (a - c + 1) / c;
a -= t * c, res += (n * (n + 1) / 2) * t;
}
if (b < 0) {
s64 t = (b - c + 1) / c;
b -= t * c, res += (n + 1) * t;
}
return res + euclid(a, b, c, n);
}
取整
向下取整：\(\lfloor x \rfloor\)，表示不超过 \(x\) 的最大整数（floor）。
向上取整：\(\lceil x \rceil\)，表示不小于 \(x\) 的最小整数（ceiling）。
注意：C++ 中的整数除法结果向零取整，详见 "0x11 整数除法上下取整"。
取整的简单性质：
- \(x - 1 < \lfloor x \rfloor \leq x\)，\(x \leq \lceil x \rceil < x + 1\)。
- \(\lfloor -x \rfloor = - \lceil x \rceil\)，\(\lceil -x \rceil = - \lfloor x \rfloor\)。
取整的常见恒等式
上下取整的相互转化：\(n, m \in \Z\)，\(m > 0\)
上下取整的嵌套恒等式：\(a, b, c \in \N_+\)
上下取整对加法的恒等式：\(A, B, C, D \in \N_+\)
数论分块
也称之为整除分块。
下取整数论分块
给出一个正整数 \(n\)，记 \(D_n\) 表示 \(\left\lfloor \frac{n}{d} \right\rfloor\) 所有可能的取值构成的集合，即 \(D_n = \{ \left\lfloor \frac{n}{d} \right\rfloor \mid 1 \leq d \leq n \}\)。
性质 1：\(\left\lfloor \frac{n}{d} \right\rfloor\) 的取值不超过 \(2\sqrt{n}\) 种，即 \(|D_n| \leq 2\sqrt{n}\)。
简单推导：
- 当 \(d \leq \sqrt{n}\) 时，因为 \(d\) 的取值不超过 \(\sqrt{n}\) 种，所以此时 \(\left\lfloor \frac{n}{d} \right\rfloor\) 的取值不超过 \(\sqrt{n}\) 种。
- 当 \(d > \sqrt{n}\) 时，因为 \(\left\lfloor \frac{n}{d} \right\rfloor \leq \sqrt{n}\)，所以此时 \(\left\lfloor \frac{n}{d} \right\rfloor\) 的取值不超过 \(\sqrt{n}\) 种。
性质 2：若 \(m \in D_n\)，则 \(D_m \subseteq D_n\)。
简单推导：设 \(m = \left\lfloor \frac{n}{k} \right\rfloor\)，此时对于任意 \(d \in \N_+\)
这意味着，如果我们想要递归地应用数论分块（类似杜教筛），在计算过程中涉及到的取值集合均为 \(D_n\)。
下取整数论分块：快速地找出 \(\left\lfloor \frac{n}{d} \right\rfloor\) 取值相同的所有区间。
由于函数 \(\frac{n}{x}\) 在区间 \((0, +\infty)\) 上单调递减，所以 \(\left\lfloor \frac{n}{x} \right\rfloor\) 取值相同的所有自变量 \(x\) 构成一段连续的区间。并且由上述分析，这样的区间不超过 \(2\sqrt{n}\) 个。
对于任意 \(i(1 \leq i \leq n)\)，需要找到一个最大的 \(j(i \leq j \leq n)\) 满足 \(\left\lfloor \frac{n}{i} \right\rfloor = \left\lfloor \frac{n}{j} \right\rfloor\)。由不等式 \(\left\lfloor \frac{n}{i} \right\rfloor \leq \frac{n}{j}\) 可得 \(\boldsymbol{j = \left\lfloor \frac{n}{\left\lfloor \frac{n}{i} \right\rfloor} \right\rfloor}\)。
根据该公式，我们就可以根据区间的左端点计算出右端点。于是依次遍历所有区间即可。
for (int x = 1, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
// 区间 [x, nx] 内，下取整的值均为 n / x
}
By the way：数论分块套数论分块的时间复杂度为 \(\mathcal{O}(n^{\frac{3}{4}})\)。
上取整数论分块
上取整数论分块：快速地找出 \(\lceil \frac{n}{d} \rceil\) 取值相同的所有区间。
由恒等式 \(\left\lceil \frac{n}{m} \right\rceil = \left\lfloor \frac{n - 1}{m} \right\rfloor + 1\) 可知，\(n\) 的上取整数论分块与 \(n - 1\) 的下取整数论分块，分出的区间是一样的。
于是套用 \(n - 1\) 的下取整数论分块，注意统计 \(d = n\) 的特殊情况。
for (int x = 1, nx; x <= n - 1; x = nx + 1) {
nx = (n - 1) / ((n - 1) / x);
// 区间 [x, nx] 内，上取整的值均为 (n - 1) / x + 1
}
// 注意统计 d = n 时的特殊情况
多维数论分块
二维数论分块：快速地找出 \(\left\lfloor \frac{n}{d} \right\rfloor, \left\lfloor \frac{m}{d} \right\rfloor\) 取值均相同的所有区间。
显然，这样的区间不超过 \(2\sqrt{n} + 2\sqrt{m}\) 个。
对于任意 \(i\)，需要找到一个最大的 \(j\) 满足 \(\left\lfloor \frac{n}{i} \right\rfloor = \left\lfloor \frac{n}{j} \right\rfloor, \left\lfloor \frac{m}{i} \right\rfloor = \left\lfloor \frac{m}{j} \right\rfloor\)。则有 \(\boldsymbol{j = \min\left\{ \left\lfloor \frac{n}{ \left\lfloor \frac{n}{i} \right\rfloor } \right\rfloor, \left\lfloor \frac{m}{ \left\lfloor \frac{m}{i} \right\rfloor } \right\rfloor \right\}}\)。
扩展到多维也是一样。
积性函数与迪利克雷卷积
积性函数：若数论函数 \(f(n)\)，对于任意互质正整数 \(x, y\)，都有 \(f(xy) = f(x)f(y)\)，则称 \(f(n)\) 为积性函数。
完全积性函数：若数论函数 \(f(n)\)，对于任意正整数 \(x, y\)，都有 \(f(xy) = f(x)f(y)\)，则称 \(f(n)\) 为完全积性函数。
积性函数的简单性质：
- 对于任意积性函数 \(f(x)\)（除了恒为 \(0\) 的函数），有 \(f(1) = 1\)。
- 对于任意积性函数 \(f(x), g(x)\)，以下的函数仍为积性函数
- 对于任意积性函数 \(f(x)\)，将 \(x\) 唯一分解成 \(x = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\)，有 \(f(x) = \prod_{i = 1}^t f(p_i^{c_i})\)。
常见积性函数（及性质）
单位函数（完全积性）
恒等函数（完全积性）
当 \(k = 1\) 时，亦可简写为 \(\mathrm{ID}(n) = n\)。
常数函数（完全积性）
约数函数（积性）
约数函数：
当 \(k = 0\) 时，\(\sigma_0(n)\) 表示 \(n\) 的约数个数函数。
约数个数函数的前缀和：
乘积的约数个数函数：
简单推导（仅证二维）：从唯一分解的角度，对每个质因子 \(p\) 单独考虑。\(p\) 在 \((x, y)\) 中的出现情况数，为 \(i, j\) 中质因子 \(p\) 的个数之和再加 \(1\)。将所有质因子 \(p\) 的贡献相乘，这正是 \(\sigma(ij)\) 的表达式。
可以扩展至更多维。
欧拉函数（积性）
欧拉函数：\(\varphi(n)\)，表示 \(1 \sim n\) 中与 \(n\) 互质的数的个数。将 \(n\) 唯一分解成 \(n = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\)，则欧拉函数为
乘积的欧拉函数：
简单推导：设 \(P(n)\) 表示 \(n\) 的质因数集合，则有 \(P(ij) = P(i) \cup P(j)\) 与 \(P(\gcd(i, j)) = P(i) \cap P(j)\)。
由容斥原理得
莫比乌斯函数（积性）
莫比乌斯函数：\(\mu(n)\)，将 \(n\) 唯一分解成 \(n = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\)，则莫比乌斯函数有如下的定义
特别地，\(\mu(1) = 1\)。
线性筛筛积性函数
对于积性函数 \(f(x)\)，若质数幂的函数值 \(f(p^c)\) 便于求出。则 \(f(x)\) 的前 \(n\) 项函数值可以通过线性筛得出。
具体地，记 \(\mathrm{low}_i\) 表示：设 \(i\) 的最小质因子为 \(p_1\)，其次数为 \(c_1\)，则 \(\mathrm{low}_i = p_1^{c_1}\)。
在线性筛的过程中，额外进行如下的维护
- 维护 \(\mathrm{low}_i\) 数组。
- 每当筛出一个质数 \(p\)，就将 \(f(p), f(p^2), \cdots, f(p^c)\) 全都求出来。
- 其余满足 \(i \neq\mathrm{low}_i\) 的 \(i\)，使用递推式 \(f(i) = f(\mathrm{low}_i) \times f\left(\frac{i}{\mathrm{low}_i}\right)\) 计算函数值。
0x34 线性筛筛积性函数.cpp：
int primeCount, prime[MaxV], low[MaxV];
int f[MaxV];
/* 线性筛筛积性函数 */
void sieve(const int &n) {
f[1] = 1;
for (int i = 2; i <= n; i ++) {
if (!low[i]) {
prime[++ primeCount] = i, low[i] = i;
for (s64 v = i, c = 1; v <= n; v *= i, c ++) {
// 计算 f[v]
}
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > n / i) break;
low[i * prime[j]] = i % prime[j] ? prime[j] : low[i] * prime[j];
if (i % prime[j] == 0) {
break;
}
}
if (low[i] < i) {
f[i] = f[low[i]] * f[i / low[i]];
}
}
}
对于一些简单的积性函数（例如 \(\varphi, \mu\)）就有更简单的筛法，不一定需要按照上述流程维护。
0x34 线性筛筛欧拉函数.cpp：
int primeCount, prime[MaxV], fac[MaxV];
int phi[MaxV];
/* 线性筛筛欧拉函数 */
void sieve(const int &n) {
phi[1] = 1;
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
phi[i] = i - 1;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
phi[i * prime[j]] = phi[i] * (i % prime[j] ? prime[j] - 1 : prime[j]);
}
}
}
0x34 线性筛筛莫比乌斯函数.cpp：
int primeCount, prime[MaxV], fac[MaxV];
int mu[MaxV];
/* 线性筛筛莫比乌斯函数 */
void sieve(const int &n) {
mu[1] = 1;
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
mu[i] = -1;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
mu[i * prime[j]] = mu[i] * (i % prime[j] ? -1 : 0);
}
}
}
迪利克雷卷积
迪利克雷卷积：对于两个数论函数 \(f, g\)，定义 \(f, g\) 的迪利克雷卷积为
迪利克雷卷积的简单性质：
- 交换律：\(f \ast g = g \ast f\)。
- 结合律：\((f \ast g) \ast h = f \ast (g \ast h)\)。
- 分配律：\(f \ast (g + h) = f \ast g + f \ast h\)。
- 单位元：\(f \ast \epsilon = \epsilon \ast f = f\)，因此 \(\epsilon\) 又被称为迪利克雷卷积的单位元。
常见迪利克雷卷积式
式 1
简单推导：由积性函数的简单性质可知，关于 \(n\) 的函数 \(f(n) = \sum_{d \mid n} \mu(d)\) 为积性函数。
- 当 \(n = 1\) 时，等式成立。
- 当 \(n > 1\) 时，将 \(n\) 唯一分解成 \(n = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\)，则有
综上，\(\sum_{d \mid n} \mu(d) = [n = 1]\)。
式 2
简单推导：由积性函数的简单性质可知，关于 \(n\) 的函数 \(f(n) = \sum_{d \mid n} \varphi(d)\) 为积性函数。
- 当 \(n = 1\) 时，等式成立。
- 当 \(n > 1\) 时，将 \(n\) 唯一分解成 \(n = p_1^{c_1}p_2^{c_2}\cdots p_t^{c_t}\)，则有
综上，\(\sum_{d \mid n} \varphi(d) = n\)。
式 3
简单推导：
杜教筛
杜教筛：用于 \(\mathcal{O}(n^{\frac{2}{3}})\) 时间内，求解一类数论函数 \(f(x)\) 的前缀和。
- \(f(x)\) 为积性函数。
- 存在辅助函数 \(g(x)\)，满足函数 \(f \ast g\) 与 \(g\) 的前缀和便于求出。
记 \(S(n) = \sum_{i = 1}^n f(i)\)。取杜教筛中的辅助函数 \(g(x)\)，则有
移项得
此时得到了一个 \(S(n)\) 关于 \(S\left( \left\lfloor \frac{n}{i} \right\rfloor \right)\) 的递推式。若函数 \(f \ast g\) 与 \(g\) 的前缀和便于求出，就可以使用数论分块处理该式。因此找到一个合适的辅助函数 \(g\) 是杜教筛的关键。
杜教筛的时间复杂度分析：
由数论分块的性质，对于 \(m \in D_n\)，有 \(D_m \subseteq D_n\)。于是我们需要使用记忆化保证杜教筛的复杂度，使用记忆化以后，只需要对所有的 \(k \in D_n\) 计算一次 \(S(k)\) 即可。
假设函数 \(f \ast g\) 与 \(g\) 的前缀和均可 \(\mathcal{O}(1)\) 求出，设计算 \(S(n)\) 的时间复杂度为 \(T(n)\)，则有
进一步优化，我们预处理出数论函数 \(f(x)\) 的前 \(m\) 项（其中 \(m \geq \sqrt{n}\)），并预处理出其前缀和。设预处理的时间复杂度为 \(T_0(m)\)，则有
通常我们使用线性筛预处理 \(f(x)\) 的前 \(m\) 项，此时 \(T_0(m) = \mathcal{O}(m)\)。由均值不等式，取 \(m = \mathcal{O}(n^{\frac{2}{3}})\) 可得杜教筛的最优时间复杂度 \(\mathcal{O}(n^{\frac{2}{3}})\)。
总结：预处理出 \(f(x)\) 的前 \(n^{\frac{2}{3}}\) 项及其前缀和，使用关于 \(S\left( \left\lfloor \frac{n}{i} \right\rfloor \right)\) 的递推式求解 \(S(n)\)，并使用 std::map 记忆化。
0x34 杜教筛.cpp：
/* 杜教筛 */
namespace xudyh {
const int V = 1e6, MaxV = V + 10; // 预处理到的范围（需根据实际问题调整）
int primeCount, prime[MaxV], low[MaxV];
int f[MaxV];
s64 fArr[MaxV];
std::map<int, s64> fBuc;
void sieve(const int &n = V) {
f[1] = 1;
for (int i = 2; i <= n; i ++) {
if (!low[i]) {
prime[++ primeCount] = i, low[i] = i;
for (s64 v = i, c = 1; v <= n; v *= i, c ++) {
// 计算 f[v]
}
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > n / i) break;
low[i * prime[j]] = i % prime[j] ? prime[j] : low[i] * prime[j];
if (i % prime[j] == 0) {
break;
}
}
if (low[i] < i) {
f[i] = f[low[i]] * f[i / low[i]];
}
}
for (int i = 1; i <= n; i ++) {
fArr[i] = fArr[i - 1] + f[i];
}
}
s64 fSum(int n) {
if (n <= V) {
return fArr[n];
}
if (fBuc.contains(n)) {
return fBuc[n];
}
s64 ans = 0; // 计算函数 f*g 第 1~n 项的和
for (int x = 2, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
s64 sum = 0; // 计算函数 g 第 x~nx 项的和
ans -= sum * fSum(n / x);
}
return fBuc[n] = ans;
}
}
杜教筛求莫比乌斯函数与欧拉函数前缀和：
给莫比乌斯函数选择的辅助函数为 \(1(n) = 1\)，由于 \(\mu \ast 1 = \epsilon\)，函数 \(\epsilon\) 与 \(1\) 的前缀和均便于求出。代入杜教筛公式得
给欧拉函数选择的辅助函数为 \(1(n) = 1\)，由于 \(\varphi \ast 1 = \mathrm{ID}\)，函数 \(\mathrm{ID}\) 与 \(1\) 的前缀和均便于求出。代入杜教筛公式得
求欧拉函数前缀和，更简便的方法是莫比乌斯反演，注意到 \(\sum_{i = 1}^n \varphi(i) = \sum_{1 \leq i \leq j \leq n} [\gcd(i, j) = 1]\)。
然而 \(\sum_{i = 1}^n \sum_{j = 1}^n [\gcd(i, j) = 1] = \sum_{d = 1}^n \mu(d)\left\lfloor \frac{n}{d} \right\rfloor^2\)，使用莫比乌斯函数前缀和配合数论分块处理该式即可。
0x34 杜教筛求莫比乌斯函数与欧拉函数前缀和.cpp：
/* 杜教筛求莫比乌斯函数与欧拉函数前缀和 */
namespace xudyh {
const int V = 1e6, MaxV = V + 10;
int primeCount, prime[MaxV], fac[MaxV];
int mu[MaxV];
int muArr[MaxV];
std::map<int, int> muBuc;
void sieve(const int &n = V) {
mu[1] = 1;
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
mu[i] = -1;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
mu[i * prime[j]] = mu[i] * (i % prime[j] ? -1 : 0);
}
}
for (int i = 1; i <= n; i ++) {
muArr[i] = muArr[i - 1] + mu[i];
}
}
int muSum(int n) {
if (n <= V) {
return muArr[n];
}
if (muBuc.contains(n)) {
return muBuc[n];
}
int ans = 1;
for (int x = 2, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
ans -= (nx - x + 1) * muSum(n / x);
}
return muBuc[n] = ans;
}
s64 phiSum(int n) {
s64 ans = 0;
for (int x = 1, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
ans += 1ll * (n / x) * (n / x) * (muSum(nx) - muSum(x - 1));
}
return (ans + 1) / 2;
}
}
Min_25 筛
Min_25 筛：用于 \(\mathcal{O}(n^{1 - \epsilon})\) 时间内，求解一类数论函数 \(f(x)\) 的前缀和。
- \(f(x)\) 为积性函数。
- 对于 \(p \in \mathbb{P}\)，\(f(p)\) 是关于 \(p\) 的低次多项式。
- 对于 \(p \in \mathbb{P}\)，\(f(p^c)\) 便于求出。
一些约定：
- 记 \(p_k\) 表示第 \(k\) 小的质数（特别地，\(p_0 = 1\)）。
- 记 \(\mathrm{lpf}(i)\) 表示 \(i\) 的最小质因子。
记 \(F(n, k) = \sum_{i = 2}^n \left[p_k < \mathrm{lpf}(i)\right] f(i)\)，记 \(F_{\mathrm{prime}}(n) = \sum_{2 \leq p \leq n, p \in \mathbb{P}} f(p)\)。
现在的目标是求出 \(\sum_{i = 1}^n f(i) = F(n, 0) + f(1)\)。
考虑求出 \(F(n, k)\)，枚举 \(2 \sim n\) 每个数的最小质因子及其次数，以计算贡献
边界情况：当 \(p_k > n\) 时，\(F(n, k) = 0\)。
考虑求出 \(F_{\mathrm{prime}}(n)\)，根据 \(F(n, k)\) 的递推式，我们发现只有 \(x \in D_n\) 处的点值 \(F_{\mathrm{prime}}(x)\) 是需要的。
若对于 \(p \in \mathbb{P}\)，\(f(p)\) 是关于 \(p\) 的低次多项式。可以表示成 \(f(p) = \sum a_ip^i\)，那么可以分开考虑每个 \(p^i\) 的贡献，其对 \(F_{\mathrm{prime}}(n)\) 的贡献为 \(a_i\left( \sum_{2 \leq p \leq n, p \in \mathbb{P}} p^i \right)\)。
对于一个固定的次数 \(c\)，记 \(G(n, k) = \sum_{i = 2}^n \left[p_k < \mathrm{lpf}(i) \lor i \in \mathbb{P}\right] i^c\)。
现在的目标是，对于所有 \(x \in D_n\) 求出 \(G(x, \lambda)\)，其中 \(p_{\lambda}\) 是最后一个不超过 \(\sqrt{n}\) 的质数。
考虑求出 \(G(n, k)\)，考虑如何从 \(G(n, k - 1)\) 转移到 \(G(n, k)\)，相当于是要去掉 \(2 \sim n\) 中满足 \(\mathrm{lpf}(i) = p_k\) 的合数 \(i\) 的贡献。先提取出一个 \(p_k\)，剩下数的范围不超过 \(\left\lfloor \frac{n}{p_k} \right\rfloor\)。于是剩下数的贡献即为 \(G\left( \left\lfloor \frac{n}{p_k} \right\rfloor, k - 1 \right)\)，但此时我们会把 \([1, p_{k - 1}]\) 范围内的质数贡献给扣除掉，所以还要补回这些质数的贡献 \(G(p_{k - 1}, k - 1)\)
注意到式子中有一项 \(G(p_{k - 1}, k - 1)\)，即质数处的点值。由于此时 \(p_{k - 1} < \lfloor \sqrt{n} \rfloor\)，可以证明 \(p_{k - 1} \in D_n\)。于是质数处的点值也会被处理出来。
引理：对于任意 \(1 \leq k < \lfloor \sqrt{n} \rfloor\)，均有 \(k \in D_n\)。
边界情况：\(G(n, 0) = \sum_{i = 2}^n i^c\)，一般可以根据公式快速算出。
0x34 Min_25 筛.cpp：
// 以求解函数 f(p^c) = p^c(p^c - 1) 为例，注意 n >= 1e10 情况下的取模与 long long 问题
/* Min_25 筛 */
namespace Min_25 {
const int V = 1e5, MaxV = V + 10;
int primeCount, prime[MaxV], fac[MaxV];
void sieve(const int &n = V) {
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
}
}
}
constexpr int inv6 = qpow(6, mod - 2, mod);
int S1(int n) {
return (1ll * n * (n + 1) / 2) % mod;
}
int S2(int n) {
return 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
int fSum(int n) {
int b = sqrt(n);
int t = std::upper_bound(prime + 1, prime + 1 + primeCount, b) - prime - 1;
std::vector<int> id1(b + 1), id2(b + 1);
auto id = [&] (int x) {
return x <= b ? id1[x] : id2[n / x];
};
std::vector<int> seq(1); // 集合 D_n，数值按照从大到小排序
std::vector<int> g1(1), g2(1);
for (int x = 1, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
int v = n / x;
if (v <= b) {
id1[v] = seq.size();
} else {
id2[n / v] = seq.size();
}
seq.push_back(v);
g1.push_back(norm(S1(v) - 1));
g2.push_back(norm(S2(v) - 1));
}
for (int j = 1; j <= t; j ++) {
int p = prime[j];
for (int i = 1; i < seq.size(); i ++) {
if (1ll * p * p > seq[i]) {
break;
}
int x1 = id(seq[i] / p), x2 = id(prime[j - 1]);
dec(g1[i], 1ll * p * (g1[x1] - g1[x2] + mod) % mod);
dec(g2[i], 1ll * p * p % mod * (g2[x1] - g2[x2] + mod) % mod);
}
}
std::function<int(int, int)> solve = [&] (int n, int j) {
if (prime[j] > n) {
return 0;
}
int x1 = id(n), x2 = id(prime[j]);
int ans = norm((g2[x1] - g1[x1]) - (g2[x2] - g1[x2]));
for (int i = j + 1; i <= t; i ++) {
int p = prime[i];
if (1ll * p * p > n) {
break;
}
for (s64 v = p; v <= n; v *= p) {
add(ans, 1ll * v * (v - 1) % mod * (solve(n / v, i) + (v > p)) % mod);
}
}
return ans;
};
return (solve(n, 0) + 1) % mod;
}
}
Min_25 筛求前缀质数个数：相当于是求出 \(2 \sim n\) 中所有质数的 \(0\) 次方，仅使用 Min_25 筛前半部分的求质数贡献即可。
0x34 Min_25 筛求前缀质数个数.cpp：
/* Min_25 筛求前缀质数个数 */
namespace Min_25 {
const int V = 1e5, MaxV = V + 10;
int primeCount, prime[MaxV], fac[MaxV];
void sieve(const int &n = V) {
for (int i = 2; i <= n; i ++) {
if (!fac[i]) {
prime[++ primeCount] = i, fac[i] = i;
}
for (int j = 1; j <= primeCount; j ++) {
if (prime[j] > fac[i] || prime[j] > n / i) break;
fac[i * prime[j]] = prime[j];
}
}
}
int fSum(int n) {
int b = sqrt(n);
int t = std::upper_bound(prime + 1, prime + 1 + primeCount, b) - prime - 1;
std::vector<int> id1(b + 1), id2(b + 1);
auto id = [&] (int x) {
return x <= b ? id1[x] : id2[n / x];
};
std::vector<int> seq(1); // 集合 D_n，数值按照从大到小排序
std::vector<int> g(1);
for (int x = 1, nx; x <= n; x = nx + 1) {
nx = n / (n / x);
int v = n / x;
if (v <= b) {
id1[v] = seq.size();
} else {
id2[n / v] = seq.size();
}
seq.push_back(v);
g.push_back(v - 1);
}
for (int j = 1; j <= t; j ++) {
int p = prime[j];
for (int i = 1; i < seq.size(); i ++) {
if (1ll * p * p > seq[i]) {
break;
}
int x1 = id(seq[i] / p), x2 = id(prime[j - 1]);
g[i] -= g[x1] - g[x2];
}
}
return g[id(n)];
}
}
同余
带余除法：对于任意 \(n, m \in \Z\) 且 \(m \neq 0\)，存在唯一的 \(q, r \in \Z\)，使得 \(n = qm + r\)（\(0 \leq r < |m|\)）。
这里记 \(q = \left\lfloor \frac{n}{m} \right\rfloor\)，\(r = n \bmod m\)。
同余：若 \(a \bmod m = b \bmod m\)，则称 \(a, b\) 模 \(m\) 同余，记作 \(a \equiv b \pmod m\)。
同余的简单性质：
- 若 \(a \equiv b \pmod m, c \equiv d \pmod m\)，则
- \(a + c \equiv b + d \pmod m\)。
- \(ac \equiv bd \pmod m\)。
- \(ac \equiv bc \pmod m \Longrightarrow a \equiv b \pmod{\frac{m}{\gcd(c, m)}}\)。
同余类：集合 \(\overline{a} = \{a + km \mid k \in \Z\}\)（其中 \(a\in [0, m - 1]\)）中的所有数模 \(m\) 同余，余数均为 \(a\)，该集合称为一个模 \(m\) 的同余类。
（完全）剩余系：\(m\) 的完全剩余系包含 \(m\) 个元素 \(\overline{0}, \overline{1}, \cdots, \overline{m - 1}\)。
（简化）剩余系：\(m\) 的简化剩余系包含 \(\varphi(m)\) 个元素，\(1 \sim m\) 中与 \(m\) 互质的数代表的同余类构成该剩余系。
剩余系的简单性质：
- 对于 \(i = 0, 1, \cdots, m - 1\)，\(ia \bmod m\) 两两不同当且仅当 \(\gcd(a, m) = 1\)。
- 对于与 \(m\) 互质的整数 \(i, j\)，\(ij \bmod m\) 仍然与 \(m\) 互质。
费马小定理与欧拉定理
费马小定理：设 \(p\) 是质数，对于任意整数 \(a\)（\(a\) 不能是 \(p\) 的倍数），满足 \(a^{p - 1} \equiv 1 \pmod p\)。
另一个形式：设 \(p\) 是质数。对于任意整数 \(a\)，满足 \(a^p \equiv a \pmod p\)。
简单推导：由于 \(p\) 是质数且 \(p \nmid a\)，对于 \(i = 0, 1, \cdots, p - 1\)，\(ia \bmod p\) 两两不同。
因此 \(\{ ia \bmod p \}\) 也构成 \(0 \sim p - 1\) 的排列。所以
又由于
所以
由于 \(1, 2, \cdots, p - 1\) 均不是 \(p\) 的倍数，所以 \(a^{p - 1} - 1\) 是 \(p\) 的倍数。即 \(a^{p - 1}\equiv 1 \pmod p\)。
By the way：费马小定理的逆命题并不成立。
即使对于所有与 \(p\) 互质的整数 \(a\) 均满足 \(a^{p - 1} \equiv 1 \pmod p\)，\(p\) 也不一定是质数。
欧拉定理：设整数 \(m > 0\)，对于任意满足 \(\gcd(a, m) = 1\) 的整数 \(a\)，满足 \(a^{\varphi(m)} \equiv 1 \pmod m\)。
简单推导：设 \(R\) 表示 \(1 \sim m\) 中与 \(m\) 互质的数构成的集合。已知 \(a\in R\)，对于所有 \(R\) 中的元素 \(r\)，\(ra \bmod m\) 仍然属于 \(R\) 且两两不同。所以
又由于
所以
由于 \(r \in R\) 均不是 \(m\) 的倍数，所以 \(a^{\varphi(m)} - 1\) 是 \(m\) 的倍数。即 \(a^{\varphi(m)} \equiv 1 \pmod m\)。
扩展欧拉定理：设整数 \(m > 0\)，对于任意整数 \(a, b(b \geq 0)\)，有
直观理解：考虑 \(a^b \bmod m\) 的循环情况。当 \(\gcd(a, m) = 1\) 时，该循环为纯循环；当 \(\gcd(a, m) \neq 1\) 时，该循环为混循环（类似 \(\rho\) 形），需要对尾巴和环形分别讨论。
通常使用扩展欧拉定理，以达到降幂的效果。
模意义下逆元
模意义下逆元：对于非零整数 \(a, m\)，如果存在整数 \(b\) 使得 \(ab \equiv 1 \pmod m\)，则称 \(b\) 为 \(a\) 在模 \(m\) 意义下的逆元。记作 \(a^{-1} \bmod p\)。
模质数意义下逆元
当模数 \(p\) 为质数时。对于任意整数 \(a\)（\(a\) 不能是 \(p\) 的倍数），有
此时，\(a^{p - 2} \bmod p\) 即为 \(a\) 在模 \(p\) 意义下的逆元。
单个逆元的求解
若要求解 \(a\) 在模 \(m\) 意义下的逆元，相当于是要求解线性同余方程 \(ax \equiv 1 \pmod m\)，使用 exgcd 即可。
时间复杂度 \(\mathcal{O}(\log V)\)。
0x34 单个逆元的求解.cpp：
/* 单个逆元的求解 */
s64 inverse(s64 a, s64 m) {
auto [d, x, y] = exgcd(a, m);
// assert(d == 1);
return (x % m + m) % m;
}
多个逆元的求解
若要求解多个整数 \(a_1, \cdots, a_n\) 在模 \(m\) 意义下的逆元，处理出前缀积与前缀积逆元，通过 \(a^{-1}_i = S_{i - 1}S^{-1}_i \bmod m\) 求出所有数的逆元即可。
时间复杂度 \(\mathcal{O}(n + \log V)\)。
0x34 多个逆元的求解.cpp：
/* 多个逆元的求解 */
auto inverse(std::vector<int> a, int mod) {
int n = a.size(), x;
std::vector<int> res(n);
x = 1;
for (int i = 0; i < n; i ++) {
res[i] = x;
x = 1ll * x * a[i] % mod;
}
x = qpow(x, mod - 2, mod); // 模数为质数时可用
for (int i = n - 1; i >= 0; i --) {
res[i] = 1ll * res[i] * x % mod;
x = 1ll * x * a[i] % mod;
}
return res;
}
预处理线性逆元
若要求解 \(1 \sim n\) 所有整数模质数 \(p\) 意义下的逆元。对于 \(1 < i < p\)，有
将等式两边同乘以 \(i^{-1}(p \bmod i)^{-1}\)，得
根据上式递推即可。
时间复杂度 \(\mathcal{O}(n)\)。
0x34 预处理线性逆元.cpp：
std::vector<int> v;
/* 预处理线性逆元 */
void linearInverse(const int &n) {
v.resize(n + 1);
v[1] = 1;
for (int i = 2; i <= n; i ++) {
v[i] = 1ll * v[mod % i] * (mod - mod / i) % mod;
}
}
线性同余方程
裴蜀定理：设 \(a, b\) 是不均为 \(0\) 的整数。
- 存在整数 \(x, y\)，使得 \(ax + by = \gcd(a, b)\)。
- 对于任意整数 \(x, y\)，均有 \(\gcd(a, b) \mid ax + by\)。
裴蜀定理可以扩展到更多维：设 \(a_1, \cdots, a_n\) 是不均为 \(0\) 的整数。
- 存在整数 \(x_1, \cdots, x_n\)，使得 \(\sum_{i = 1}^n a_ix_i = \gcd(a_1, \cdots, a_n)\)。
- 对于任意整数 \(x_1, \cdots, x_n\)，均有 \(\gcd(a_1, \cdots, a_n) \mid \sum_{i = 1}^n a_ix_i\)。
可以使用 gcd 的结合律进行归纳证明。
exgcd（扩展欧几里得算法）：用于求出二元一次不定方程 \(ax + by = \gcd(a, b)\) 的一组可行解 \(x, y\)。
在欧几里得算法的基础上归纳求解：
- 当 \(b = 0\) 时。此时有 \(x = 1, y = 0\)，使得 \(a \cdot 1 + 0 \cdot 0 = \gcd(a, 0)\)。
- 当 \(b > 0\) 时。由于 \(\gcd(a, b) = \gcd(b, a \bmod b)\)，假设存在一对 \(x, y\) 使得
利用 \(a \bmod b = a - \left\lfloor \frac{a}{b} \right\rfloor b\)，并移项得
此时取 \(x' = y, y' = x - \left\lfloor \frac{a}{b} \right\rfloor y\)，使得 \(ax' + by' = \gcd(a, b)\)。
时间复杂度 \(\mathcal{O}(\log \max(a, b))\)。
可以证明，扩展欧几里得算法求出的可行解 \(x, y\) 满足 \(|x| \leq b, |y| \leq a\)。
0x34 exgcd（扩展欧几里得算法）.cpp：
/* 扩展欧几里得算法 */
std::array<s64, 3> exgcd(s64 a, s64 b) {
if (b == 0) return {a, 1, 0};
auto [d, x, y] = exgcd(b, a % b);
return {d, y, x - (a / b) * y};
}
二元一次不定方程 \(ax + by = c\)：
- 有解性：当且仅当 \(\gcd(a, b) \mid c\)。
- 通解：使用 exgcd 找出方程 \(ax + by = \gcd(a, b)\) 的一组特解 \(x_0, y_0\)，则有通解
线性同余方程 \(ax \equiv b \pmod m\)：
该方程等价于 \(ax - b\) 为 \(m\) 的倍数，不妨设为 \(-y\) 倍，则原方程转化为 \(ax + my = b\)。
- 有解性：当且仅当 \(\gcd(a, m) \mid b\)。
- 通解：使用 exgcd 找出方程 \(ax + my = \gcd(a, m)\) 的一组特解 \(x_0, y_0\)，则有通解
将线性同余方程中 \(x\) 的系数消成 \(1\)：
对于线性同余方程 \(ax \equiv b \pmod m\)，有解的充要条件是 \(\gcd(a, m) \mid b\)。当满足该条件时，方程可以转化为
\[\frac{a}{\gcd(a, m)}x \equiv \frac{b}{\gcd(a, m)} \pmod{\frac{m}{\gcd(a, m)}} \]此时 \(\gcd\left( \frac{a}{\gcd(a, m)}, \frac{m}{\gcd(a, m)} \right) = 1\)，\(x\) 在模 \(\frac{m}{\gcd(a, m)}\) 意义下有唯一解 \(x'\)。
于是线性同余方程 \(ax \equiv b \pmod m\) 转化为了 \(x \equiv x' \pmod{m'}\) 的形式。该形式对线性同余方程组的求解大有用处。
线性同余方程组
模数两两互质的一元线性同余方程组（中国剩余定理，CRT）：
- 记所有模数的积为 \(M = \prod_{i = 1}^n m_i\)。
- 对于第 \(i\) 个方程：
- 记 \(M_i = \frac{M}{m_i}\)。
- 记 \(V_i\) 表示 \(M_i\) 在模 \(m_i\) 意义下的逆元。
- 特解：
- 通解：
0x34 CRT（中国剩余定理）.cpp：
// {a, m} 表示方程 x = a (mod m)
using Equ = std::pair<s64, s64>;
/* 中国剩余定理 */
Equ CRT(std::vector<Equ> seq) {
auto inverse = [&] (s64 a, s64 m) -> s64 {
auto [d, x, y] = exgcd(a, m);
return (x % m + m) % m;
};
s64 M = 1;
for (auto [a, m] : seq) {
M *= m;
}
s64 x = 0;
for (auto [a, m] : seq) {
s64 Mi = M / m;
s64 Vi = inverse(Mi, m);
x = (x + a * Mi * Vi) % M; // 必要时，使用“快速乘”防止爆 long long
}
x = (x % M + M) % M;
return {x, M};
}
不保证模数两两互质的一元线性同余方程组（扩展中国剩余定理，exCRT）：
仅考虑两个方程的合并（多个方程的情况，可以依次合并）。设这两个方程为
则
使用 exgcd 找出方程 \(m_1p - m_2q = a_2 - a_1\) 的一组特解 \(p_0, q_0\)。
此时取 \(a' = p_0m_1 + a_1\) 与 \(m' = \operatorname{lcm}(m_1, m_2)\)，两个方程合并为
0x34 exCRT（扩展中国剩余定理）.cpp：
// {a, m} 表示方程 x = a (mod m)
using Equ = std::pair<s64, s64>;
/* 扩展中国剩余定理 */
Equ exCRT(std::vector<Equ> seq) {
auto [a1, m1] = seq[0];
for (int i = 1; i < seq.size(); i ++) {
auto [a2, m2] = seq[i];
auto [d, x, y] = exgcd(m1, m2);
if ((a2 - a1) % d) {
return {-1, -1};
}
x = x * ((a2 - a1) / d) % (m2 / d); // 必要时，使用“快速乘”防止爆 long long
s64 m = m1 * (m2 / d);
s64 a = (x * m1 + a1) % m; // 必要时，使用“快速乘”防止爆 long long
a1 = a, m1 = m;
}
a1 = (a1 % m1 + m1) % m1;
return {a1, m1};
}
阶与原根
阶：设整数 \(m > 0\)，满足 \(a^x \equiv 1 \pmod m\) 的最小正整数 \(x\)，称作 \(a\) 模 \(m\) 的阶。记作 \(\mathrm{ord}_m(a)\)。
阶的存在条件：当且仅当 \(\gcd(a, m) = 1\) 时，\(a\) 模 \(m\) 的阶存在。
必要性：当 \(\gcd(a, m) \neq 1\) 时，\(a\) 在模 \(m\) 意义下的逆元不存在。故必定有 \(\gcd(a, m) = 1\)。
充分性：当 \(\gcd(a, m) = 1\) 时，由欧拉定理 \(a^{\varphi(m)} \equiv 1 \pmod m\)，故此时阶必定存在，并且 \(\mathrm{ord}_m(a) \mid \varphi(m)\)。
阶的简单性质：
- 若 \(a\perp m\)，则 \(a, a^2, \cdots, a^{\mathrm{ord}_m(a)}\) 模 \(m\) 两两不同。
- 若 \(a \perp m\)，则 \(\mathrm{ord}_m(a) \mid \varphi(m)\)。
- 若 \(a \perp m\)，则
简单推导：\(\mathrm{ord}_m(a^k) = \frac{\mathrm{lcm}(\mathrm{ord}_m(a), k)}{k}\)。
- 若 \(a\perp m, b\perp m\)，则
- 若 \(a \perp m, b \perp m\)，则
- 若 \(a \perp m, b\perp m\)，则总是存在满足 \(c\perp m\) 的 \(c\)，使得
简单构造：将 \(\mathrm{ord}_m(a), \mathrm{ord}_m(b)\) 唯一分解
\[\mathrm{ord}_m(a) = \prod p_i^{s_i}, \quad \mathrm{ord}_m(b) = \prod p_i^{t_i} \]根据 \(s_i\) 与 \(t_i\) 的大小关系，可以将所有质因子分成两类
\[A = \{ s_i \geq t_i \mid i \}, \quad B = \{ s_i < t_i \mid i \} \]记
\[x_A = \prod_{i \in A} p_i^{s_i}, \quad x_B = \prod_{i \in B} p_i^{s_i}, \quad y_A = \prod_{i \in A} p_i^{t_i}, \quad y_B = \prod_{i \in B} p_i^{t_i} \]此时 \(\mathrm{ord}_m(a) = x_Ax_B, \mathrm{ord}_m(b) = y_Ay_B\)，可得
\[\mathrm{ord}_m(a^{x_B}) = x_A, \quad \mathrm{ord}_m(b^{y_A}) = y_B \]又由于 \(x_A \perp y_B\)，所以
\[\mathrm{ord}_m(a^{x_B}b^{y_A}) = x_Ay_B = \prod p_i^{\max(s_i, t_i)} = \mathrm{lcm}(\mathrm{ord}_m(a), \mathrm{ord}_m(b)) \]取 \(c = a^{x_B}b^{y_A}\) 即可。
原根：设整数 \(m > 0\)，若 \(\mathrm{ord}_m(g) = \varphi(m)\)，则称 \(g\) 为模 \(m\) 的原根。
此时 \(g, g^2, \cdots, g^{\varphi(m)}\) 构成了 \(m\) 的简化剩余系。
原根判定定理：设整数 \(m \geq 3\)，对于整数 \(g \perp m\)，若 \(g\) 是模 \(m\) 的原根，当且仅当对于 \(\varphi(m)\) 的每个质因子 \(p\)，都有
原根存在定理：\(m\) 存在原根，当且仅当 \(m = 1, 2, 4, p^c, 2p^c\)（其中 \(p\) 为奇质数且 \(c \geq 1\)）。
阶的数量：若 \(m\) 存在原根，则对于 \(d \mid \varphi(m)\)，阶等于 \(d\) 的元素恰好有 \(\varphi(d)\) 个。
简单推导：设 \(m\) 的原根为 \(g\)，则 \(m\) 的简化剩余系中的元素可以写成 \(g^k \bmod m\) 的形式，其中 \(1 \leq k \leq \varphi(m)\)。则
\[\mathrm{ord}_m(g^k) = \frac{\varphi(m)}{\gcd(\varphi(m), k)} \]令 \(\mathrm{ord}_m(g^k) = d\)，则有 \(\gcd(\varphi(m), k) = \frac{\varphi(m)}{d}\)。整个式子除以 \(\frac{\varphi(m)}{d}\) 可得 \(\gcd(d, k') = 1\)，此时 \(k'\) 的数量为 \(\varphi(d)\)。
原根数量：若 \(m\) 存在原根，则 \(m\) 的原根个数为 \(\varphi(\varphi(m))\)。
相当于是阶等于 \(\varphi(m)\) 的元素个数，由上述分析，恰好有 \(\varphi(\varphi(m))\) 个。
最小原根的上界：若 \(m\) 存在原根，则 \(m\) 的最小原根的上界为 \(m^{0.25}\) 级别。
原根求解：先求出 \(\varphi(m)\)，再求出 \(\varphi(m)\) 的唯一分解形式。此时就可以根据 "原根判定定理" 配合快速幂判断一个数是否为 \(m\) 的原根。
从小到大暴力寻找原根，即可找到最小原根 \(g\)。
当我们已知其中一个原根 \(g\) 时，可以根据 \(g\) 生成其他所有的原根。对于所有满足 \(k \perp \varphi(m)\) 的 \(k\)，这些 \(g^k \bmod m\) 不重不漏地覆盖了所有原根。
离散对数
离散对数：设整数 \(m > 0\) 存在原根 \(g\)，对于 \(m\) 的简化剩余系中的元素 \(a\)，必存在唯一的整数 \(0 \leq k < \varphi(m)\) 使得 \(g^k \equiv a \pmod m\)。
此时称 \(k\) 为以 \(g\) 为底，模 \(m\) 的离散对数。记作 \(k = \mathrm{ind}_g(a)\)。
离散对数的简单性质：
- \(\mathrm{ind}_g(ab) \equiv \mathrm{ind}_g(a) + \mathrm{ind}_g(b) \pmod{\varphi(m)}\)。
说明离散对数可以将 "模 \(m\) 意义下乘法" 转换成 "模 \(\varphi(m)\) 意义下加法"。 - \(\mathrm{ind}_g(a^k) \equiv k\operatorname{ind}_g(a) \pmod{\varphi(m)}\)。
- \(\mathrm{ind}_{g_1}(a) \equiv \operatorname{ind}_{g_2}(a) \operatorname{ind}_{g_1}(g_2) \pmod{\varphi(m)}\)。其中 \(g_1, g_2\) 均为 \(m\) 的原根。
高次同余方程
固定底数且与模数互质的高次同余方程 \(a^x \equiv b \pmod{p}\)（大步小步算法，BSGS）：
取 \(S = \left\lceil \sqrt{p} \right\rceil\)，设 \(m = iS - j\)（其中 \(i \geq 1\) 且 \(0 \leq j < S\)），则有
预处理 \(0 \leq j < S\) 的 \(b\times a^j \bmod p\) 存入哈希表，枚举所有 \(1 \leq i \leq S\)，在哈希表中查询 \((a^S)^i \bmod p\) 及其对应的 \(j\)。如果查询成功则返回 \(iS - j\)。如果一直查询失败则方程无解。
时间复杂度 \(\mathcal{O}(\sqrt{p})\)。
0x34 BSGS（大步小步算法）.cpp：
/* 大步小步算法 */
int BSGS(int a, int b, int p) { // a^x = b (mod p)
if (b == 1 || p == 1) return 0;
if (a == 0) return b == 0 ? 1 : -1;
int S = sqrt(p) + 1;
std::unordered_map<int, int> buc;
int w = 1;
for (int j = 0; j < S; j ++, w = 1ll * w * a % p) {
int v = 1ll * w * b % p;
buc[v] = j;
}
int x = w;
for (int i = 1; i <= S; i ++, x = 1ll * x * w % p) {
if (buc.contains(x)) {
return i * S - buc[x];
}
}
return -1;
}
固定底数且不保证与模数互质的高次同余方程 \(a^x \equiv b \pmod{p}\)（扩展大步小步算法，exBSGS）：
该方程等价于 \(a^x - b\) 为 \(p\) 的倍数，不妨设为 \(-k\) 倍，则原方程转化为 \(a\cdot a^{x - 1} + kp = b\)。
若 \(\gcd(a, p) \nmid b\) 则无解，否则
每次将 \(\frac{a}{\gcd(a, p)}\) 提取出来，并将 \(b, p\) 除以 \(\gcd(a, p)\)，就变成了一个规模更小的子问题。
记已经提取出的系数之积为 \(r = \prod \frac{a}{\gcd(a, p)}\)。在递降的过程中，若 \(r \equiv b' \pmod{p'}\)，则直接返回答案。
否则一直递降直到 \(\gcd(a, p') = 1\)，此时得到一个固定底数且与模数互质的高次同余方程 \(ra^x \equiv b' \pmod{p'}\)，套用 BSGS。
0x34 exBSGS（扩展大步小步算法）.cpp：
/* 扩展大步小步算法 */
int exBSGS(int a, int b, int p) { // a^x = b (mod p)
if (b == 1 || p == 1) return 0;
if (a == 0) return b == 0 ? 1 : -1;
int d, r = 1, step = 0;
while (d = std::gcd(a, p), d > 1) {
if (b % d) {
return -1;
}
b /= d, p /= d, r = 1ll * r * (a / d) % p;
step ++;
if (r == b) {
return step;
}
}
if (p == 1) {
return step;
}
int S = sqrt(p) + 1;
std::unordered_map<int, int> buc;
int w = 1;
for (int j = 0; j < S; j ++, w = 1ll * w * a % p) {
int v = 1ll * w * b % p;
buc[v] = j;
}
int x = 1ll * r * w % p;
for (int i = 1; i <= S; i ++, x = 1ll * x * w % p) {
if (buc.contains(x)) {
return step + i * S - buc[x];
}
}
return -1;
}
固定次数的高次同余方程 \(x^a \equiv b \pmod p\)：当 \(p\) 存在原根 \(g\) 时，记 \(x = g^c\)（这里 \(c\) 是未知数），问题转化为
记 \(v = g^a\)，此时得到一个固定底数的高次同余方程 \(v^c \equiv b \pmod p\)。
- 特解：\(x \equiv g^{ac} \pmod p\)，使用 (ex)BSGS 得到。
- 通解：\(x \equiv g^{ac + k\frac{\varphi(p)}{\gcd(\varphi(p), a)}} \pmod{p}\)，其中 \(\frac{\varphi(p)}{\gcd(\varphi(p), a)}\) 表示 \(\mathrm{ord}_m(g^a)\)。
Misc
平方差
平方差表示：当且仅当 \(1, 2, 4\)，以及所有形如 \(4k + 2\) 的正整数，不能被表示成两个不同正整数的平方差。
平方差表示证明：对于正整数 \(n\)，相当于判断是否存在两个不同的正整数 \(a, b(a > b)\) 使得 \(n = a^2 - b^2 = (a + b)(a - b)\)。
设 \(x = a + b, y = a - b\)，由于 \(x, y\) 之间的差为 \(2b\)，所以 \(x, y\) 奇偶性相同。对于同奇偶的两个 \(x, y(x > y)\)，也都可以解出 \(a, b\) 的值。
故对于正整数 \(n\)，相当于是要找出 \(n\) 的两个不同的约数 \(x, y\)，且 \(x, y\) 奇偶性相同。
- 当 \(n\) 为奇数时，取 \(x = n, y = 1\) 即可。
- 当 \(n\) 为偶数时：
- 当 \(n = 4k\) 时，取 \(x = 2k, y = 2\) 即可。
- 当 \(n = 4k + 2\) 时，分解出的两个约数，必为一奇一偶。
- \(n = 1, 2, 4\) 为特殊情况，均不能被表示。
0x35 数学知识 离散数学
集合论基础
集合并：\(A \cup B = \{ x \mid x \in A \ \operatorname{or} \ x \in B \}\)。
集合交：\(A \cap B = \{ x \mid x \in A \ \operatorname{and} \ x \in B \}\)。
集合补：\(\complement_U A = \{ x \mid x \notin A \ \operatorname{and} \ x \in U \}\)。
集合相对补：\(B \setminus A = \{ x \mid x \notin A \ \operatorname{and} \ x \in B \}\)。
集合对称差：\(A \operatorname{\triangle} B = \{ x \mid x \in A \ \operatorname{xor} \ x \in B \}\)。同时
集合笛卡尔积：\(A \times B = \{(x, y) \mid x \in A, y \in B \}\)。
集合反演律（德 · 摩根定律）：
置换
置换：有限集合到自身的双射（即一一对应）称为置换。不可重集合 \(S = \{a_1, a_2, \cdots, a_n\}\) 上的置换可以表示为
表示将所有 \(a_i\) 映射成 \(a_{p_i}\)，即 \(f(a_i) = a_{p_i}\)。其中 \(p_1, p_2, \cdots, p_n\) 为 \(1 \sim n\) 的一个排列。
置换的乘法（复合）：给出两个置换
\(f\) 与 \(g\) 的乘法（复合）记作 \(f \circ g\)
即 \((f \circ g)(x) = f(g(x))\)，即先经过了 \(g\) 的映射再经过了 \(f\) 的映射（从右到左）。
置换的逆：给出置换
\(f\) 的逆置换记作 \(f^{-1}\)
一个关于 \(n\) 的排列，与其每个元素排名的序列，互为逆排列。
对称群：集合 \(S\) 上的所有置换构成一个群：满足封闭性、结合律、有单位元（恒等置换）、存在逆元（逆置换）。
通常记集合 \(\{1, 2, \cdots, n\}\) 上的所有置换构成的群为 \(n\) 元对称群，记作 \(S_n\)。
置换群：对称群 \(S_n\) 的任意一个子群，即为置换群。
轮换：一类特殊的置换，可表示为
置换的轮换表示：任意一个置换，都可以分解成若干个不相交的轮换的乘积。轮换可以看作是构成置换的基本单位。
轮换的几何意义：将元素视为图中的点，将映射关系视为图中的有向边，由于每个节点的入度与出度均为 \(1\)，因此形成的图必定由若干个不相交的环构成。
1-轮换（不动点）：对于集合 \(X\) 上的置换 \(\sigma\)，常用 \(X^\sigma\) 表示 \(\sigma\) 的不动点集合。即 \(X^{\sigma} = \{x \in X \mid \sigma(x) = x \}\)。
2-轮换（对换）：
- 任何置换都可以写作一系列对换的乘积。
- 进一步，任何置换都可以写作一系列相邻对换的乘积（冒泡排序）。
置换的奇偶性
排列的奇偶性：排列的奇偶性，为排列逆序对数的奇偶性。
- 交换排列的任意两个相邻元素 \(p_i, p_{i + 1}\)，会使得整个排列的奇偶性取反。
- 交换排列的任意两个不同元素 \(p_i, p_j\)，可以使用 \(2|i - j| - 1\) 次邻项交换实现。故交换排列的任意两个不同的元素，也会使得整个排列的奇偶性取反。
置换的奇偶性：置换的奇偶性，为将该置换分解成若干个对换，所需对换个数的奇偶性。
群/环/域 基础
群
| 群类型 | 性质 | 举例 |
|---|---|---|
| 群 | 封闭性、结合律、一一一、单位元、逆元 | 对称群 \(S_n\)（集合 \(\{1, 2, \cdots, n\}\) 上的所有置换） |
| 半群 | 封闭性、结合律、一一一、一一一、一一 | \((\N_+, +)\) |
| 幺半群 | 封闭性、结合律、一一一、单位元、一一 | \((\Z, \times)\) |
| 交换群 | 封闭性、结合律、交换律、单位元、逆元 | \((\Z, +)\) |
群：设 \(G\) 为非空集合，其上有 \(G \times G \to G\)（封闭性）的二元运算 \(\cdot\)，如果满足以下性质，则称 \((G, \cdot)\) 是一个群：
- 结合律：对于任意 \(a, b, c \in G\)，成立 \(a \cdot (b \cdot c) = (a \cdot b) \cdot c\)。
- 有单位元：存在 \(e \in G\)，使得对于任意 \(a \in G\)，都成立 \(a \cdot e = e \cdot a = a\)。这里称 \(e\) 为 \(G\) 的单位元。
- 存在逆元：对于任意 \(a \in G\)，都存在相应的 \(b \in G\) 使得 \(a \cdot b = b \cdot a = e\)。这里称 \(b\) 为 \(a\) 的逆元。
群的基本性质：
- 结合律推论：对于任意有限长的列 \(g_1, g_2, \cdots, g_k\)，乘积 \(g_1 \cdot g_2 \cdot \ldots \cdot g_k\) 的结果与加括号的方式无关。
- 单位元唯一：单位元 \(e\) 总是唯一的。
- 逆元唯一：对于任意 \(x \in G\)，其逆元 \(x^{-1}\) 也是唯一的。
- 消去律：对于任意 \(a, b, c \in G\)，若 \(a \cdot c = b \cdot c\) 或 \(c \cdot a = c \cdot b\)，则 \(a = b\)。
半群
半群：设 \(G\) 为非空集合，其上有 \(G \times G \to G\)（封闭性）的二元运算 \(\cdot\)，若其满足结合律，则称 \((G, \cdot)\) 是一个半群。
幺半群
幺半群：设半群 \((G, \cdot)\)，若其有单位元，则称 \((G, \cdot)\) 是一个幺半群。
数据结构中，线段树等结构可以维护幺半群的信息。
交换群（Abel 群）
交换群：设群 \((G, \cdot)\)，若其满足交换律，则称 \((G, \cdot)\) 是一个交换群。
环
| 环类型 | 性质 | 举例 |
|---|---|---|
| 环 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、一一一、一一一、一一。 分配律 |
|
| 幺环 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、一一一、单位元、一一。 分配律 |
多项式环 \(R[x]\) |
| 除环 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、一一一、单位元、逆元。 分配律 |
四元数环 \(\mathbb{H}\) |
| 交换环 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、交换律、一一一、一一。 分配律 |
\(2\Z\) |
| 整环 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、交换律、单位元、一一。 分配律、无零因子 |
\((\Z, +, \times)\) |
环：设 \(R\) 为非空集合，其上有两个 \(R\times R \to R\)（封闭性）的二元运算 \(+, \cdot\)，如果满足以下性质，则称 \((R, +, \cdot)\) 是一个环：
- \((R, +)\) 构成交换群：即 \(+\) 满足结合律，有单位元（记作 \(0\)），存在逆元（记作 \(-a\)）。
- \((R, \cdot)\) 构成半群：即 \(\cdot\) 满足结合律。
- 分配律：对于任意 \(a, b, c \in R\)，有 \(a \cdot (b + c) = a \cdot b + a \cdot c\) 与 \((a + b) \cdot c = a \cdot c + b \cdot c\)。
为方便表述，记 $+, \cdot $ 分别表示该环的加法和乘法。加法单位元称作零元，乘法单位元（如果存在）称作幺元。
幺环
幺环：设环 \((R, +, \cdot)\)，若其存在乘法单位元（记作 \(1\)），则称 \((R, +, \cdot)\) 是一个幺环。
除环
除环：设幺环 \((R, +, \cdot)\)，若其所有非 \(0\) 元素 \(a \in R\) 都存在乘法逆元 \(a^{-1}\)，则称 \((R, +, \cdot)\) 是一个除环。
交换环
交换环：设环 \((R, +, \cdot)\)，若其满足乘法交换律，则称 \((R, +, \cdot)\) 是一个交换环。
整环
零元乘以任何元素都得到零元（即 \(0 \cdot a = a \cdot 0 = 0\)）。理解一般环的乘法结构时，要去除加法单位元的影响，考察 \(R \setminus \{0\}\)。基于这一想法，有如下定义：
- 零因子：对于环 \((R, +, \cdot)\) 中的非零元素 \(a \in R\)，若存在非零 \(b \in R\) 满足 \(a \cdot b = 0\) 或 \(b \cdot a = 0\)，则称 \(a\) 为一个零因子。
- 可逆元：对于环 \((R, +, \cdot)\) 中的非零元素 \(a \in R\)，若存在 \(b \in R\) 满足 \(a \cdot b = b \cdot a = 1\)，则称 \(a\) 为一个可逆元。
零因子不可能是可逆元，可逆元不可能是零因子。即零因子与可逆元互斥。但是一个非零元素可以既不是零因子，也不是可逆元。
整环：设交换幺环 \((R, +, \cdot)\)，且无零因子，则称 \((R, +, \cdot)\) 为整环。
整环的性质：
- 整环的消去律：对于任意 \(a, b, c\in R\) 且 \(a \neq 0\)，若 \(a \cdot b = a \cdot c\)，则 \(b = c\)。
- 非零元素对乘法封闭，但元素未必可逆，故仅构成半群。
乘法群
乘法群（单位群）：设幺环 \((R, +, \cdot)\)，设 \(R^\times\) 为 \(R\) 中全体可逆元的集合，则 \((R^\times, \cdot)\) 构成群，成为幺环 \(R\) 的乘法群。
域
| 域类型 | 性质 | 举例 |
|---|---|---|
| 域 | 加法：封闭性、结合律、交换律、单位元、逆元。 乘法：封闭性、结合律、交换律、单位元、逆元。 分配律 |
数域 \(\Q, \R, \C\) |
域（交换除环）：设 \(F\) 为非空集合，其上有两个 \(F \times F \to F\)（封闭性）的二元运算 \(+, \cdot\)，如果满足以下性质，则称 \((F, +, \cdot)\) 是一个域：
- \((F, +)\) 构成一个交换群。其单位元记作 \(0\)，元素 \(a \in F\) 在 \(+\) 下的逆元记作 \(-a\)。
- \((F \setminus \{0\}, \cdot)\) 构成一个交换群。其单位元记作 \(1\)，元素 \(a \in F\) 在 \(\cdot\) 下的逆元记作 \(a^{-1}\)。
域是对加、减、乘、除四则运算都封闭的代数结构。非零元素不仅封闭，且构成乘法群。
浙公网安备 33010602011771号