第十四章 数论与数学¶
数论与数学是 ACM 竞赛中出题频率最高的专题之一。从质数筛法、GCD/LCM 到组合数学、博弈论，掌握这些基础知识和算法是解决大量竞赛题目的前提。本章系统梳理核心知识点，配合推导与代码实现，帮助选手建立完整的数论知识体系。
14.1 质数与筛法¶
14.1.1 质数判定¶
定义：大于 1 的正整数，除了 1 和它本身外没有其他因数的数称为质数（素数）。
试除法判断质数：对于 \(n\)，只需检查 \(2 \sim \lfloor\sqrt{n}\rfloor\) 范围内是否有 \(n\) 的因子。
#include <iostream>
using namespace std;
// 判断 n 是否为质数
// 时间复杂度: O(sqrt(n))
bool isPrime(long long n) {
if (n < 2) return false;
// 只需检查到 sqrt(n)
for (long long i = 2; i * i <= n; i++) {
if (n % i == 0) return false;
}
return true;
}
int main() {
long long n;
cin >> n;
cout << (isPrime(n) ? "质数" : "非质数") << endl;
return 0;
}
14.1.2 埃拉托斯特尼筛法（埃氏筛）¶
基本思想：从 2 开始，将每个质数的所有倍数标记为合数。
graph TD
subgraph 埃氏筛过程 - 求 2~20 的质数
A["初始化: 2~20 全部标记为质数"] --> B["p=2: 标记 4,6,8,10,12,14,16,18,20 为合数"]
B --> C["p=3: 标记 9,15 为合数（6,12,18已被标记）"]
C --> D["p=5: 无新标记（25>20）"]
D --> E["剩余质数: 2, 3, 5, 7, 11, 13, 17, 19"]
end#include <iostream>
#include <cstring>
using namespace std;
const int MAXN = 10000010;
bool isComp[MAXN]; // isComp[i] = true 表示 i 是合数
int primes[MAXN]; // 存储找到的质数
int pcnt; // 质数计数
// 埃拉托斯特尼筛法
// 求 [2, n] 范围内的所有质数
// 时间复杂度: O(n log log n)
void sieveEratosthenes(int n) {
memset(isComp, 0, sizeof(isComp));
pcnt = 0;
for (int i = 2; i <= n; i++) {
if (!isComp[i]) {
primes[pcnt++] = i; // i 是质数
// 从 i*i 开始标记，因为 i*2, i*3, ..., i*(i-1)
// 已经被更小的质数标记过了
for (long long j = (long long)i * i; j <= n; j += i) {
isComp[j] = true;
}
}
}
}
int main() {
int n;
cin >> n;
sieveEratosthenes(n);
cout << n << " 以内共有 " << pcnt << " 个质数" << endl;
return 0;
}
埃氏筛的优化
- 内层循环从 \(i^2\) 而非 \(2i\) 开始，因为 \(i \cdot 2, i \cdot 3, \ldots, i \cdot (i-1)\) 已经被更小的质数筛过
- 可以只筛奇数（除了 2），将空间和时间减半
- 实现简单，常数小，在 \(n \leq 10^7\) 时效率很好
14.1.3 欧拉筛（线性筛）¶
埃氏筛中，合数可能被其多个质因子重复标记。欧拉筛确保每个合数只被其最小质因子筛掉一次，达到线性时间复杂度 \(O(n)\)。
核心思想：遍历每个数 \(i\)，用已知的质数 \(primes[j]\) 去筛 \(i \cdot primes[j]\)。当 \(i \bmod primes[j] == 0\) 时停止，保证每个合数只被最小质因子筛掉。
#include <iostream>
#include <cstring>
using namespace std;
const int MAXN = 10000010;
bool isComp[MAXN];
int primes[MAXN];
int pcnt;
// 欧拉筛（线性筛）
// 时间复杂度: O(n)
void sieveEuler(int n) {
memset(isComp, 0, sizeof(isComp));
pcnt = 0;
for (int i = 2; i <= n; i++) {
if (!isComp[i]) {
primes[pcnt++] = i; // i 是质数
}
// 用已知质数去筛 i * primes[j]
for (int j = 0; j < pcnt && (long long)i * primes[j] <= n; j++) {
isComp[i * primes[j]] = true;
// 关键: i % primes[j] == 0 时停止
// 此时 primes[j] 是 i 的最小质因子
// 继续筛下去会导致重复标记
if (i % primes[j] == 0) break;
}
}
}
int main() {
int n;
cin >> n;
sieveEuler(n);
cout << n << " 以内共有 " << pcnt << " 个质数" << endl;
return 0;
}
为什么 \(i \bmod primes[j] == 0\) 时要停止？
设 \(i = k \cdot primes[j]\)。若继续用 \(primes[j+1]\) 筛 \(i \cdot primes[j+1] = k \cdot primes[j] \cdot primes[j+1]\)，这个合数会在后续 \(i' = k \cdot primes[j+1]\) 时被 \(primes[j]\) 再筛一次。为了保证每个数只被筛一次，在此处停止。
| 筛法 | 时间复杂度 | 空间复杂度 | 特点 |
|---|---|---|---|
| 试除法（单个判断） | \(O(\sqrt{n})\) | \(O(1)\) | 适合判断单个大数 |
| 埃氏筛 | \(O(n \log \log n)\) | \(O(n)\) | 实现简单，常数小 |
| 欧拉筛 | \(O(n)\) | \(O(n)\) | 线性，不重复标记 |
14.2 质因数分解¶
14.2.1 基本分解¶
算术基本定理：每个大于 1 的正整数都能唯一地表示为质数的乘积（不考虑顺序）：
试除法分解：从小到大尝试每个可能的因子。
#include <iostream>
#include <vector>
using namespace std;
// 质因数分解
// 返回: {质因子, 指数} 的列表
// 时间复杂度: O(sqrt(n))
vector<pair<long long, int>> factorize(long long n) {
vector<pair<long long, int>> factors;
for (long long i = 2; i * i <= n; i++) {
if (n % i == 0) {
int cnt = 0;
while (n % i == 0) {
n /= i;
cnt++;
}
factors.push_back({i, cnt});
}
}
// 如果 n > 1，说明 n 本身是一个质因子
if (n > 1) {
factors.push_back({n, 1});
}
return factors;
}
int main() {
long long n;
cin >> n;
auto factors = factorize(n);
cout << n << " = ";
for (int i = 0; i < (int)factors.size(); i++) {
if (i > 0) cout << " * ";
cout << factors[i].first;
if (factors[i].second > 1) cout << "^" << factors[i].second;
}
cout << endl;
return 0;
}
14.2.2 利用筛法预处理最小质因子¶
配合欧拉筛，可以预处理 \([1, n]\) 中每个数的最小质因子（LPF），之后任意数的分解只需 \(O(\log n)\)。
#include <iostream>
#include <vector>
using namespace std;
const int MAXN = 10000010;
int lpf[MAXN]; // lpf[i] = i 的最小质因子
int primes[MAXN], pcnt;
// 预处理最小质因子
void preprocess(int n) {
lpf[1] = 1;
for (int i = 2; i <= n; i++) {
if (lpf[i] == 0) {
lpf[i] = i;
primes[pcnt++] = i;
}
for (int j = 0; j < pcnt && primes[j] <= lpf[i] && (long long)i * primes[j] <= n; j++) {
lpf[i * primes[j]] = primes[j];
}
}
}
// 利用 lpf 数组快速分解质因数
// 时间复杂度: O(log n)
vector<pair<int, int>> fastFactorize(int n) {
vector<pair<int, int>> factors;
while (n > 1) {
int p = lpf[n], cnt = 0;
while (n % p == 0) {
n /= p;
cnt++;
}
factors.push_back({p, cnt});
}
return factors;
}
int main() {
preprocess(10000000);
int n;
cin >> n;
auto factors = fastFactorize(n);
for (auto& [p, e] : factors) {
cout << p << "^" << e << " ";
}
cout << endl;
return 0;
}
质因数分解的竞赛应用
- 因数个数：\(n = p_1^{a_1} \cdots p_k^{a_k}\) 的因数个数为 \((a_1+1)(a_2+1)\cdots(a_k+1)\)
- 因数之和：\(\sigma(n) = \prod_{i=1}^{k} \frac{p_i^{a_i+1} - 1}{p_i - 1}\)
- 欧拉函数：\(\varphi(n) = n \prod_{i=1}^{k} \left(1 - \frac{1}{p_i}\right)\)
14.3 GCD / LCM 与扩展欧几里得¶
14.3.1 辗转相除法（GCD）¶
欧几里得算法：\(\gcd(a, b) = \gcd(b, a \bmod b)\)，直到 \(b = 0\) 时，\(\gcd(a, 0) = a\)。
// 辗转相除法求最大公约数
// 时间复杂度: O(log min(a, b))
long long gcd(long long a, long long b) {
return b == 0 ? a : gcd(b, a % b);
}
// 最小公倍数
long long lcm(long long a, long long b) {
return a / gcd(a, b) * b; // 先除后乘防溢出
}
14.3.2 扩展欧几里得算法¶
问题：求整数 \(x, y\) 使得 \(ax + by = \gcd(a, b)\)。
原理推导：
当 \(b = 0\) 时，\(\gcd(a, 0) = a\)，此时 \(x = 1, y = 0\)。
当 \(b \neq 0\) 时，设已求出 \(bx' + (a \bmod b)y' = \gcd(b, a \bmod b)\)。
由于 \(a \bmod b = a - \lfloor a/b \rfloor \cdot b\)，代入得：
因此：\(x = y'\)，\(y = x' - \lfloor a/b \rfloor \cdot y'\)。
graph TD
subgraph 扩展欧几里得递推
A["ax + by = gcd(a, b)"] --> B{"b == 0 ?"}
B -->|是| C["gcd = a, x = 1, y = 0"]
B -->|否| D["递归求 bx' + (a%b)y' = gcd"]
D --> E["x = y', y = x' - (a/b)*y'"]
end#include <iostream>
using namespace std;
// 扩展欧几里得算法
// 求 ax + by = gcd(a, b) 的一组整数解
// 返回值: gcd(a, b)
// 时间复杂度: O(log min(a, b))
long long exgcd(long long a, long long b, long long& x, long long& y) {
if (b == 0) {
x = 1;
y = 0;
return a;
}
long long g = exgcd(b, a % b, y, x);
// 注意: 这里 y 和 x 交换了
// 实际上 x = y', y = x' - (a/b) * y'
y -= a / b * x;
return g;
}
int main() {
long long a, b, x, y;
cin >> a >> b;
long long g = exgcd(a, b, x, y);
cout << "gcd = " << g << endl;
cout << "x = " << x << ", y = " << y << endl;
cout << a << " * " << x << " + " << b << " * " << y << " = " << g << endl;
return 0;
}
14.3.3 解线性同余方程¶
问题：求 \(ax \equiv b \pmod{m}\) 的解。
转化为：\(ax + my = b\)，即求解不定方程。当且仅当 \(\gcd(a, m) \mid b\) 时有解。
// 解 ax ≡ b (mod m)
// 返回: 最小非负整数解，无解返回 -1
long long solveModEq(long long a, long long b, long long m) {
long long x, y;
long long g = exgcd(a, m, x, y);
if (b % g != 0) return -1; // 无解
// 调整解到 [0, m/g) 范围内
x = x * (b / g) % (m / g);
if (x < 0) x += m / g;
return x;
}
扩展欧几里得的重要性
扩展欧几里得算法是数论中的核心工具。它不仅用于求 GCD，还是求解模逆元、中国剩余定理等高级算法的基础。竞赛中经常需要灵活运用它来解决各种模方程问题。
14.4 快速幂与矩阵快速幂¶
14.4.1 快速幂¶
问题：计算 \(a^b \bmod m\)。
朴素方法需要 \(O(b)\) 次乘法。快速幂利用二进制分解将复杂度降至 \(O(\log b)\)。
原理：将指数 \(b\) 写成二进制形式 \(b = (c_k c_{k-1} \ldots c_1 c_0)_2\)，则：
graph LR
subgraph 快速幂示例: 计算 3^13
A["13 = 1101 二进制"] --> B["3^13 = 3^8 * 3^4 * 3^1"]
B --> C["3^1 = 3"]
C --> D["3^2 = 9"]
D --> E["3^4 = 81"]
E --> F["3^8 = 6561"]
F --> G["结果 = 6561 * 81 * 3 = 1594323"]
end#include <iostream>
using namespace std;
// 快速幂: 计算 a^b mod m
// 时间复杂度: O(log b)
long long qpow(long long a, long long b, long long m) {
long long res = 1 % m;
a %= m;
while (b > 0) {
// 如果 b 的当前最低位为 1，则乘上当前的 a
if (b & 1) {
res = res * a % m;
}
// a 自乘，对应指数的下一位
a = a * a % m;
// b 右移一位
b >>= 1;
}
return res;
}
int main() {
long long a, b, m;
cin >> a >> b >> m;
cout << a << "^" << b << " mod " << m << " = " << qpow(a, b, m) << endl;
return 0;
}
14.4.2 矩阵快速幂¶
问题：计算矩阵 \(A^b\)，其中 \(A\) 是 \(n \times n\) 矩阵。
将快速幂中的标量乘法替换为矩阵乘法即可。常用于加速线性递推。
经典应用：计算 Fibonacci 数列第 \(n\) 项。
#include <iostream>
#include <cstring>
#include <vector>
using namespace std;
typedef long long ll;
typedef vector<vector<ll>> Matrix;
const int MOD = 1e9 + 7;
// 矩阵乘法: C = A * B
// A: n*m 矩阵, B: m*p 矩阵
// 时间复杂度: O(n*m*p)
Matrix mul(const Matrix& A, const Matrix& B) {
int n = A.size(), m = B.size(), p = B[0].size();
Matrix C(n, vector<ll>(p, 0));
for (int i = 0; i < n; i++) {
for (int k = 0; k < m; k++) {
if (A[i][k] == 0) continue; // 优化常数
for (int j = 0; j < p; j++) {
C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % MOD;
}
}
}
return C;
}
// 矩阵快速幂: A^b
// 时间复杂度: O(n^3 log b)
Matrix qpow(Matrix A, ll b) {
int n = A.size();
Matrix res(n, vector<ll>(n, 0));
for (int i = 0; i < n; i++) res[i][i] = 1; // 单位矩阵
while (b > 0) {
if (b & 1) res = mul(res, A);
A = mul(A, A);
b >>= 1;
}
return res;
}
// 利用矩阵快速幂求 Fibonacci 数列第 n 项
ll fib(ll n) {
if (n == 0) return 0;
if (n == 1) return 1;
Matrix A = {{1, 1}, {1, 0}};
Matrix res = qpow(A, n - 1);
return res[0][0]; // F(n)
}
int main() {
ll n;
cin >> n;
cout << "F(" << n << ") = " << fib(n) << endl;
return 0;
}
矩阵快速幂的使用场景
- 线性递推关系（Fibonacci、Tribonacci 等）
- 图的路径计数（\(A^k\) 的 \((i,j)\) 元素表示 \(i\) 到 \(j\) 恰好 \(k\) 步的路径数）
- 带状态转移的 DP 加速
- 关键：递推必须是线性的，才能用矩阵表示
14.5 模逆元与费马小定理¶
14.5.1 模逆元的定义¶
定义：对于 \(a\) 和模数 \(m\)，若 \(\gcd(a, m) = 1\)，则存在整数 \(x\) 使得 \(ax \equiv 1 \pmod{m}\)。\(x\) 称为 \(a\) 模 \(m\) 的逆元，记为 \(a^{-1}\)。
用途：在模运算中，除法不能直接做。\(a / b \bmod m\) 需要转化为 \(a \cdot b^{-1} \bmod m\)。
14.5.2 费马小定理¶
定理：若 \(p\) 是质数，且 \(\gcd(a, p) = 1\)，则：
由此可得：
利用快速幂计算 \(a^{p-2} \bmod p\)，时间复杂度 \(O(\log p)\)。
#include <iostream>
using namespace std;
const int MOD = 1e9 + 7;
// 快速幂
long long qpow(long long a, long long b, long long m) {
long long res = 1 % m;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
// 利用费马小定理求模逆元（要求 MOD 为质数）
// 时间复杂度: O(log MOD)
long long invFermat(long long a, long long p = MOD) {
return qpow(a, p - 2, p);
}
int main() {
long long a;
cin >> a;
cout << a << "^(-1) mod " << MOD << " = " << invFermat(a) << endl;
// 验证: a * a^(-1) % MOD 应该等于 1
cout << "验证: " << a << " * " << invFermat(a) << " mod " << MOD
<< " = " << a % MOD * invFermat(a) % MOD << endl;
return 0;
}
14.5.3 线性递推求逆元¶
当需要求 \(1, 2, \ldots, n\) 所有数的逆元时，逐一用快速幂效率为 \(O(n \log p)\)。使用递推公式可优化到 \(O(n)\)。
递推公式（设 \(p\) 为质数）：
推导：设 \(p = k \cdot i + r\)，其中 \(k = \lfloor p/i \rfloor\)，\(r = p \bmod i\)。则 \(k \cdot i + r \equiv 0 \pmod{p}\)，两边同乘 \(i^{-1} \cdot r^{-1}\) 得 \(k \cdot r^{-1} + i^{-1} \equiv 0\)，即 \(i^{-1} \equiv -k \cdot r^{-1} \pmod{p}\)。
#include <iostream>
using namespace std;
const int MAXN = 3000010;
const int MOD = 1e9 + 7;
long long inv[MAXN];
// 线性递推求 [1, n] 所有数的模逆元
// 要求 MOD 为质数
// 时间复杂度: O(n)
void preprocessInv(int n) {
inv[1] = 1;
for (int i = 2; i <= n; i++) {
inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
}
}
int main() {
int n;
cin >> n;
preprocessInv(n);
for (int i = 1; i <= n; i++) {
cout << i << "^(-1) mod " << MOD << " = " << inv[i] << endl;
}
return 0;
}
| 方法 | 时间复杂度 | 适用场景 |
|---|---|---|
| 费马小定理 | \(O(\log p)\) 每个 | 模数为质数，求单个逆元 |
| 扩展欧几里得 | \(O(\log p)\) 每个 | 模数不必为质数 |
| 线性递推 | \(O(n)\) 总计 | 批量求 \(1 \sim n\) 的逆元 |
14.6 组合数学¶
14.6.1 排列与组合¶
排列数：从 \(n\) 个不同元素中取 \(r\) 个排列：
组合数：从 \(n\) 个不同元素中取 \(r\) 个（不考虑顺序）：
递推关系（杨辉三角）：
边界条件：\(\binom{n}{0} = \binom{n}{n} = 1\)。
14.6.2 预处理组合数（模意义下）¶
利用阶乘和逆元，可以在 \(O(n)\) 预处理后 \(O(1)\) 查询组合数。
#include <iostream>
using namespace std;
const int MAXN = 2000010;
const int MOD = 1e9 + 7;
long long fac[MAXN]; // 阶乘: fac[i] = i! mod MOD
long long invFac[MAXN]; // 阶乘的逆元: invFac[i] = (i!)^(-1) mod MOD
// 快速幂
long long qpow(long long a, long long b, long long m) {
long long res = 1 % m;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
// 预处理阶乘和逆元
// 时间复杂度: O(n)
void preprocessComb(int n) {
fac[0] = 1;
for (int i = 1; i <= n; i++) {
fac[i] = fac[i - 1] * i % MOD;
}
invFac[n] = qpow(fac[n], MOD - 2, MOD); // 费马小定理求逆元
for (int i = n - 1; i >= 0; i--) {
invFac[i] = invFac[i + 1] * (i + 1) % MOD; // 递推求逆元
}
}
// 组合数 C(n, r) mod MOD
// 时间复杂度: O(1) 查询
long long C(int n, int r) {
if (r < 0 || r > n) return 0;
return fac[n] % MOD * invFac[r] % MOD * invFac[n - r] % MOD;
}
int main() {
preprocessComb(2000000);
int n, r;
cin >> n >> r;
cout << "C(" << n << ", " << r << ") = " << C(n, r) << endl;
return 0;
}
14.6.3 Lucas 定理¶
当 \(n, r\) 很大但模数 \(p\) 较小（且为质数）时，使用 Lucas 定理。
定理：设 \(p\) 为质数，将 \(n, r\) 写成 \(p\) 进制：
则：
其中当 \(r_i > n_i\) 时，\(\binom{n_i}{r_i} = 0\)。
适用条件：\(p\) 为质数，\(p\) 不太大（通常 \(p \leq 10^5\)），\(n, r\) 可以很大。
#include <iostream>
using namespace std;
long long MOD;
// 快速幂
long long qpow(long long a, long long b, long long m) {
long long res = 1 % m;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
// 预处理阶乘（模 p 意义下）
const int MAXP = 100010;
long long fac[MAXP];
void preprocessFac() {
fac[0] = 1;
for (int i = 1; i < MOD; i++) {
fac[i] = fac[i - 1] * i % MOD;
}
}
// 小范围组合数 C(n, r) mod p，用费马小定理
long long combSmall(long long n, long long r) {
if (r < 0 || r > n) return 0;
return fac[n] % MOD * qpow(fac[r], MOD - 2, MOD) % MOD
* qpow(fac[n - r], MOD - 2, MOD) % MOD;
}
// Lucas 定理
// 求 C(n, r) mod p，p 为质数
// 时间复杂度: O(p + log_p(n) * log p)
long long lucas(long long n, long long r) {
if (r == 0) return 1;
return lucas(n / MOD, r / MOD) * combSmall(n % MOD, r % MOD) % MOD;
}
int main() {
long long n, r;
cin >> n >> r >> MOD;
preprocessFac();
cout << "C(" << n << ", " << r << ") mod " << MOD << " = " << lucas(n, r) << endl;
return 0;
}
Lucas 定理的适用范围
- 模数 \(p\) 必须是质数
- 当模数不是质数但可以分解为质数幂 \(p^k\) 时，需要使用扩展 Lucas 定理（exLucas），实现更复杂
- 当 \(n, r < p\) 时，直接用预处理的阶乘 + 逆元即可，无需 Lucas
14.7 中国剩余定理（CRT）与扩展 CRT¶
14.7.1 中国剩余定理（CRT）¶
问题：求解同余方程组：
CRT 定理：当 \(m_1, m_2, \ldots, m_k\) 两两互质时，方程组在 \(\bmod\ M\)（\(M = m_1 m_2 \cdots m_k\)）意义下有唯一解。
解法：
- 计算 \(M = m_1 \cdot m_2 \cdots m_k\)
- 对每个 \(i\)，计算 \(M_i = M / m_i\)
- 求 \(M_i\) 模 \(m_i\) 的逆元 \(t_i\)（因为 \(m_i\) 两两互质，\(\gcd(M_i, m_i) = 1\)，逆元存在）
- 解为 \(x = \sum_{i=1}^{k} a_i \cdot M_i \cdot t_i \bmod M\)
#include <iostream>
#include <vector>
using namespace std;
// 扩展欧几里得（用于求逆元）
long long exgcd(long long a, long long b, long long& x, long long& y) {
if (b == 0) { x = 1; y = 0; return a; }
long long g = exgcd(b, a % b, y, x);
y -= a / b * x;
return g;
}
// 中国剩余定理
// a[i]: 余数, m[i]: 模数（两两互质）
// 返回: x mod M，其中 M = m[0]*m[1]*...*m[k-1]
long long crt(const vector<long long>& a, const vector<long long>& m) {
int k = a.size();
long long M = 1;
for (int i = 0; i < k; i++) M *= m[i];
long long x = 0;
for (int i = 0; i < k; i++) {
long long Mi = M / m[i];
long long ti, y;
exgcd(Mi, m[i], ti, y); // Mi * ti ≡ 1 (mod m[i])
ti = (ti % m[i] + m[i]) % m[i];
x = (x + a[i] * Mi % M * ti % M) % M;
}
return x;
}
int main() {
int k;
cin >> k;
vector<long long> a(k), m(k);
for (int i = 0; i < k; i++) cin >> a[i] >> m[i];
cout << "x = " << crt(a, m) << endl;
return 0;
}
CRT 中的乘法溢出风险
上面代码中 x = (x + a[i] * Mi % M * ti % M) % M 这一步，a[i] * Mi 与后续的 * ti 都是在模 \(M\) 之前先做乘法。当 \(M\) 接近 long long 上限（约 \(9.2 \times 10^{18}\)）时，两个 \([0, M)\) 内的数相乘会直接溢出，得到错误结果。只要 \(M > 3 \times 10^9\) 左右就有溢出可能，竞赛中模数乘积很容易超过这个范围。
解决方法是用 __int128 快速乘（mulmod）替换代码中的乘法：
// __int128 快速乘: 计算 a * b % m, 中间结果用 128 位整数存放, 不会溢出
long long mulmod(long long a, long long b, long long m) {
return (__int128)a * b % m;
}
然后把累加一行改写为：
__int128 在主流评测机（GCC/Clang，64 位环境）下均可用。若比赛环境不支持，可改用「龟速乘」（类似快速幂的加法版本，\(O(\log b)\)）或 long double 技巧。14.7.2 的 EXCRT 中 cur_m * t 一步同样存在此风险，处理方式相同。
14.7.2 扩展中国剩余定理（EXCRT）¶
当 \(m_i\) 不一定两两互质时，不能直接使用 CRT。EXCRT 的方法是逐个合并方程。
合并两个方程：将 \(x \equiv a_1 \pmod{m_1}\) 和 \(x \equiv a_2 \pmod{m_2}\) 合并。
设 \(x = m_1 k + a_1\)，代入第二个方程：\(m_1 k + a_1 \equiv a_2 \pmod{m_2}\)，即 \(m_1 k \equiv a_2 - a_1 \pmod{m_2}\)。
用扩展欧几里得解此方程，得到新的 \(a\) 和 \(m = \text{lcm}(m_1, m_2)\)。
#include <iostream>
#include <vector>
using namespace std;
long long exgcd(long long a, long long b, long long& x, long long& y) {
if (b == 0) { x = 1; y = 0; return a; }
long long g = exgcd(b, a % b, y, x);
y -= a / b * x;
return g;
}
// 扩展 CRT: 合并 x ≡ a (mod m) 和当前解 x ≡ cur_a (mod cur_m)
// 返回: 是否有解; 合并后 cur_a, cur_m 更新
bool merge(long long& cur_a, long long& cur_m, long long a, long long m) {
long long x, y;
long long g = exgcd(cur_m, m, x, y);
// 方程: cur_m * x + m * y = g
// 需要 cur_m * t ≡ (a - cur_a) (mod m)
long long diff = a - cur_a;
if (diff % g != 0) return false; // 无解
// t = x * (diff / g) mod (m / g)
long long t = x * (diff / g) % (m / g);
if (t < 0) t += m / g;
// 更新: 新 a = cur_a + cur_m * t, 新 m = lcm(cur_m, m)
cur_a = cur_a + cur_m * t;
cur_m = cur_m / g * m; // lcm
cur_a = (cur_a % cur_m + cur_m) % cur_m;
return true;
}
// 扩展中国剩余定理
// 返回: 最小非负整数解，无解返回 -1
long long excrt(const vector<long long>& a, const vector<long long>& m) {
int k = a.size();
long long cur_a = a[0], cur_m = m[0];
for (int i = 1; i < k; i++) {
if (!merge(cur_a, cur_m, a[i], m[i])) return -1;
}
return cur_a;
}
int main() {
int k;
cin >> k;
vector<long long> a(k), m(k);
for (int i = 0; i < k; i++) cin >> a[i] >> m[i];
long long ans = excrt(a, m);
if (ans == -1) cout << "无解" << endl;
else cout << "x = " << ans << endl;
return 0;
}
| 方法 | 模数要求 | 时间复杂度 |
|---|---|---|
| CRT | 两两互质 | \(O(k \log M)\) |
| EXCRT | 无特殊要求 | \(O(k \log M)\) |
14.8 容斥原理¶
14.8.1 基本原理¶
容斥原理是计数中的经典方法。设 \(A_1, A_2, \ldots, A_n\) 是有限集合，则：
简记为：
graph TD
subgraph 容斥原理示意 - 三个集合
A["|A ∪ B ∪ C|"] --> B["= |A| + |B| + |C|"]
B --> C["- |A∩B| - |A∩C| - |B∩C|"]
C --> D["+ |A∩B∩C|"]
end14.8.2 经典应用：与 n 互质的数的个数¶
求 \([1, n]\) 中与 \(n\) 互质的数的个数，即欧拉函数 \(\varphi(n)\)。
设 \(n\) 的质因子为 \(p_1, p_2, \ldots, p_k\)，则利用容斥原理：
#include <iostream>
#include <vector>
using namespace std;
// 欧拉函数（利用质因数分解 + 容斥）
// 时间复杂度: O(sqrt(n))
long long eulerPhi(long long n) {
long long result = n;
for (long long i = 2; i * i <= n; i++) {
if (n % i == 0) {
// i 是 n 的质因子
while (n % i == 0) n /= i;
result = result / i * (i - 1); // result *= (1 - 1/p)
}
}
if (n > 1) result = result / n * (n - 1);
return result;
}
int main() {
long long n;
cin >> n;
cout << "phi(" << n << ") = " << eulerPhi(n) << endl;
return 0;
}
14.8.3 容斥原理的位运算实现¶
当需要枚举集合的所有子集进行容斥时，可以用二进制枚举。
#include <iostream>
#include <vector>
using namespace std;
// 容斥原理: 求 [1, n] 中至少能被给定质数之一整除的数的个数
// primes: 给定的质数列表
// n: 上界
// 时间复杂度: O(2^k * k)，k 为质数个数
long long inclusionExclusion(const vector<int>& primes, long long n) {
int k = primes.size();
long long ans = 0;
// 枚举所有非空子集（用二进制表示）
for (int mask = 1; mask < (1 << k); mask++) {
long long product = 1;
int bits = 0; // 子集大小
for (int i = 0; i < k; i++) {
if (mask & (1 << i)) {
bits++;
product *= primes[i];
}
}
// 奇数个集合取交 → 加；偶数个 → 减
if (bits & 1) {
ans += n / product;
} else {
ans -= n / product;
}
}
return ans;
}
int main() {
long long n;
int k;
cin >> n >> k;
vector<int> primes(k);
for (int i = 0; i < k; i++) cin >> primes[i];
cout << inclusionExclusion(primes, n) << endl;
return 0;
}
容斥原理的使用技巧
- 当集合数量 \(k\) 较小（通常 \(k \leq 20\)）时，可以直接二进制枚举 \(2^k\) 个子集
- 当 \(k\) 较大时，需要结合其他技巧（如按质因子分组）
- 容斥原理经常与 GCD、LCM、互质性等问题结合出题
- 计算欧拉函数时，直接用公式 \(n \prod (1 - 1/p_i)\) 比容斥更高效
14.9 博弈论¶
14.9.1 博弈论基础概念¶
组合博弈（Combinatorial Game）满足以下条件的两人博弈：
- 两人轮流操作
- 信息完全公开（完美信息）
- 无随机因素
- 有限步内结束，且没有平局
- 最后无法操作的一方输（Normal Play Convention）
P-position（必败态）：前一手的玩家（Previous player）获胜，即当前玩家必败。
N-position（必胜态）：后一手的玩家（Next player）获胜，即当前玩家必胜。
判定规则：
- 终态（无法操作）是 P-position
- 能到达至少一个 P-position 的状态是 N-position
- 所有后继状态都是 N-position 的状态是 P-position
14.9.2 Nim 游戏¶
规则：有 \(n\) 堆石子，每堆分别有 \(a_1, a_2, \ldots, a_n\) 个。两人轮流操作，每次从某一堆中取走至少一个石子。取走最后一个石子的人获胜。
Sprague-Grundy 定理（Nim 和）：
其中 \(\oplus\) 表示按位异或（XOR）。
#include <iostream>
#include <vector>
using namespace std;
// Nim 游戏判定
// 返回: true 表示先手必胜，false 表示先手必败
bool nimGame(int n, const vector<int>& a) {
int xorSum = 0;
for (int i = 0; i < n; i++) {
xorSum ^= a[i];
}
return xorSum != 0; // 异或和不为 0 则先手必胜
}
int main() {
int n;
cin >> n;
vector<int> a(n);
for (int i = 0; i < n; i++) cin >> a[i];
cout << (nimGame(n, a) ? "先手必胜" : "先手必败") << endl;
return 0;
}
14.9.3 SG 函数（Sprague-Grundy 函数）¶
对于更一般的博弈，每个状态 \(x\) 可以转移到若干后继状态。定义 SG 函数：
其中 mex（Minimum EXcluded value）是"最小不在集合中的非负整数"。
Sprague-Grundy 定理：多个独立子游戏的 SG 值等于各子游戏 SG 值的异或和。
#include <iostream>
#include <cstring>
#include <set>
#include <vector>
using namespace std;
const int MAXN = 10010;
int sg[MAXN]; // SG 函数值
bool vis[MAXN]; // 辅助数组，用于计算 mex
int moves[MAXN]; // 可行操作（如每次取 1,2,3 个石子）
int moveCnt; // 可行操作的数量
// 计算 SG 函数（记忆化搜索）
// x: 当前状态（剩余石子数）
int getSG(int x) {
if (sg[x] != -1) return sg[x]; // 已计算
set<int> nextSG; // 后继状态的 SG 值集合
for (int i = 0; i < moveCnt; i++) {
if (x >= moves[i]) {
nextSG.insert(getSG(x - moves[i]));
}
}
// 计算 mex
int mex = 0;
while (nextSG.count(mex)) mex++;
return sg[x] = mex;
}
// 多堆博弈，利用 SG 定理
bool sgGame(int n, const vector<int>& piles) {
int xorSum = 0;
for (int i = 0; i < n; i++) {
xorSum ^= getSG(piles[i]);
}
return xorSum != 0;
}
int main() {
memset(sg, -1, sizeof(sg));
sg[0] = 0; // 终态 SG 值为 0
// 示例: 每次可以取 1, 2, 或 3 个石子
moveCnt = 3;
moves[0] = 1; moves[1] = 2; moves[2] = 3;
int n;
cin >> n;
vector<int> piles(n);
for (int i = 0; i < n; i++) cin >> piles[i];
cout << (sgGame(n, piles) ? "先手必胜" : "先手必败") << endl;
return 0;
}
14.9.4 经典博弈模型¶
| 博弈类型 | 规则 | SG 值 | 判定 |
|---|---|---|---|
| 取石子（Nim） | 从任意一堆取任意个 | \(\text{SG}(x) = x\) | 异或和非 0 先手胜 |
| Bash Game | 从一堆取 \(1 \sim m\) 个 | \(\text{SG}(x) = x \bmod (m+1)\) | \(n \bmod (m+1) \neq 0\) 先手胜 |
| Wythoff Game | 两堆，取同数量或从一堆取 | SG 值与黄金分割有关 | 特征对 \((\lfloor k\phi \rfloor, \lfloor k\phi^2 \rfloor)\) |
| Fibonacci Game | 取 \(1 \sim 2k\) 个（\(k\) 为上次取的数） | 与 Fibonacci 数列相关 | \(n\) 非 Fibonacci 数时先手胜 |
博弈论的核心思路
- 找规律：先从小数据手算，尝试发现模式
- SG 函数：将复杂博弈分解为独立子游戏，计算各自的 SG 值
- 异或合并：多个独立子游戏的 SG 值取异或
- 特判终态：SG(终态) = 0
- 许多博弈题的难点不在于计算 SG，而在于识别博弈结构和找规律
14.10 算法对比与总结¶
graph TD
subgraph 数论算法选择指南
A{"需要解决什么问题?"}
A -->|"质数相关"| B["筛法 + 质因数分解"]
A -->|"GCD / 线性方程"| C["欧几里得 + 扩展欧几里得"]
A -->|"快速幂运算"| D["快速幂 / 矩阵快速幂"]
A -->|"模运算中的除法"| E["模逆元"]
A -->|"组合计数"| F["组合数预处理 / Lucas"]
A -->|"同余方程组"| G["CRT / EXCRT"]
A -->|"计数/互质"| H["容斥原理 / 欧拉函数"]
A -->|"博弈/决策"| I["SG 函数 / Nim"]
A -->|"巨大指数降幂"| J["欧拉定理 / 扩展欧拉定理"]
A -->|"⌊n/i⌋ 求和"| K["数论分块"]
A -->|"离散对数"| L["BSGS"]
end| 算法 | 适用场景 | 时间复杂度 | 空间复杂度 |
|---|---|---|---|
| 埃氏筛 / 欧拉筛 | 求范围内的质数 | \(O(n \log \log n)\) / \(O(n)\) | \(O(n)\) |
| 质因数分解 | 分解整数 | \(O(\sqrt{n})\) | \(O(\log n)\) |
| GCD / EXGCD | 最大公约数、线性方程 | \(O(\log n)\) | \(O(1)\) |
| 快速幂 | 幂运算 | \(O(\log b)\) | \(O(1)\) |
| 矩阵快速幂 | 线性递推加速 | \(O(d^3 \log n)\) | \(O(d^2)\) |
| 模逆元 | 模意义下除法 | \(O(\log p)\) / \(O(n)\) | \(O(1)\) / \(O(n)\) |
| 组合数 + Lucas | 组合计数 | \(O(1)\) 查询 / \(O(\log_p n)\) | \(O(n)\) / \(O(p)\) |
| CRT / EXCRT | 同余方程组 | \(O(k \log M)\) | \(O(1)\) |
| 容斥原理 | 集合计数 | \(O(2^k)\) | \(O(k)\) |
| SG 函数 | 博弈论 | \(O(n \cdot m)\) | \(O(n)\) |
| 欧拉定理 / 扩展欧拉定理 | 巨大指数降幂 | \(O(\sqrt{m} + \log b)\) | \(O(1)\) |
| 线性筛求 \(\varphi\) | 批量求欧拉函数 | \(O(n)\) | \(O(n)\) |
| 数论分块 | \(\sum \lfloor n/i \rfloor\) 类求和 | \(O(\sqrt{n})\) | \(O(1)\) |
| BSGS | 离散对数 \(a^x \equiv b \pmod{p}\) | \(O(\sqrt{p})\) | \(O(\sqrt{p})\) |
拓展方向
表中后四行的内容见本章 14.11 ~ 14.14 的进阶小节。此外，高斯消元可在 \(O(n^3)\) 内求解线性方程组与异或方程组（开关问题、线性基相关），属于本章的拓展内容，此处仅点到为止，遇到相关题目可查阅线性代数专题资料。
14.11 欧拉定理与扩展欧拉定理¶
14.11.1 从欧拉函数说起¶
欧拉函数 \(\varphi(n)\)：\([1, n]\) 中与 \(n\) 互质的数的个数（定义与容斥推导见 14.8.2）。由质因数分解 \(n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}\) 可得：
单点求值直接使用 14.8.2 中的 eulerPhi 函数（\(O(\sqrt{n})\) 试除实现）；需要批量求值时用 14.12 的线性筛。
14.11.2 欧拉定理¶
定理：若 \(\gcd(a, m) = 1\)，则
当 \(m\) 为质数 \(p\) 时 \(\varphi(p) = p - 1\)，欧拉定理退化为 14.5.2 的费马小定理，因此它是费马小定理在任意模数上的推广。
降幂推论：\(\gcd(a, m) = 1\) 时，指数可以直接对 \(\varphi(m)\) 取模：
14.11.3 扩展欧拉定理（任意底数的降幂）¶
当 \(\gcd(a, m) \neq 1\) 时上面的推论不再成立。例如 \(2^4 = 16 \equiv 0 \pmod{8}\)，但 \(2^{4 \bmod \varphi(8)} = 2^0 = 1\)。
扩展欧拉定理：对任意整数 \(a\)，只要 \(b \geq \log_2 m\)，就有
实现时通常使用更方便的判据「\(b \geq \varphi(m)\)」（该条件下结论同样成立），即分段处理：
典型套路——巨大指数：求 \(2^{10^{10^5}} \bmod p\) 这类题时，指数本身有十万位，任何整数类型都存不下。做法是把 \(b\) 当作十进制字符串逐位读入，一边读一边对 \(\varphi(m)\) 取模，并用一个布尔量记录「\(b\) 是否曾经 \(\geq \varphi(m)\)」，最后按上式决定是否补加 \(\varphi(m)\)，再用快速幂收尾。
例题：洛谷 P5091 【模板】扩展欧拉定理（https://www.luogu.com.cn/problem/P5091）
题意：给定 \(a, m, b\)，求 \(a^b \bmod m\)。其中 \(1 \leq a \leq 10^9\)，\(1 \leq m \leq 10^8\)，\(1 \leq b \leq 10^{2 \times 10^7}\)（\(b\) 以十进制字符串形式给出）。
#include <iostream>
#include <string>
using namespace std;
typedef long long ll;
// O(sqrt(n)) 求欧拉函数（同 14.8.2 的 eulerPhi）
ll eulerPhi(ll n) {
ll result = n;
for (ll i = 2; i * i <= n; i++) {
if (n % i == 0) {
while (n % i == 0) n /= i;
result = result / i * (i - 1);
}
}
if (n > 1) result = result / n * (n - 1);
return result;
}
// 快速幂
ll qpow(ll a, ll b, ll m) {
ll res = 1 % m;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
int main() {
ll a, m;
string b; // 指数 b 可能有几万位，以字符串读入
cin >> a >> m >> b;
ll phi = eulerPhi(m);
ll r = 0;
bool big = false; // 记录是否有 b >= phi(m)
for (char c : b) {
r = r * 10 + (c - '0');
if (r >= phi) { big = true; r %= phi; }
}
if (big) r += phi; // 扩展欧拉定理: b >= phi(m) 时指数用 b mod phi(m) + phi(m)
cout << qpow(a, r, m) << endl;
return 0;
}
验证记录（点击展开）
以上程序已编译运行验证：
- 输入
2 7 10（\(b \geq \varphi(7)=6\) 分支）：输出2，与 \(2^{10} = 1024 \equiv 2 \pmod 7\) 一致 - 输入
3 10 2（\(b < \varphi(10)=4\) 分支）：输出9 - 输入 \(a=2\)、\(m=10^9+7\)、\(b=10^{100}\)（101 位大指数）：输出
314344290，与 Pythonpow(2, 10**100, 10**9+7)对拍一致
扩展欧拉定理常见坑
- \(b < \varphi(m)\) 时绝不能加 \(\varphi(m)\)：此时直接计算 \(a^b\)，无条件补加会得到错误结果
- 「只取模不加 \(\varphi(m)\)」仅在 \(\gcd(a, m) = 1\) 时合法；写通用模板时统一按扩展欧拉定理的分段形式处理最稳妥
- 幂塔 \(a^{b^{c}} \bmod m\) 需要逐层递归降幂：第二层模数为 \(\varphi(m)\)、第三层为 \(\varphi(\varphi(m))\)……\(\varphi\) 迭代 \(O(\log m)\) 次后必然到 1，递归深度有限
- \(m = 1\) 时任何数模 1 都是 0，模板中
1 % m的写法可自然处理该边界
14.12 线性筛求欧拉函数¶
14.12.1 为什么能在线性筛中顺带求出 φ¶
单点 \(O(\sqrt{n})\) 求 \(\varphi\)，批量求 \([1, n]\) 全部取值就要 \(O(n \sqrt{n})\)，太慢。注意到 \(\varphi\) 是积性函数（\(\gcd(a,b)=1 \Rightarrow \varphi(ab) = \varphi(a)\varphi(b)\)），且 \(\varphi(p^k) = p^k - p^{k-1}\)。而 14.1.3 的欧拉筛恰好保证每个合数只被最小质因子筛一次，这正是积性函数递推需要的全部信息。
设 \(p = primes[j]\) 为当前用于筛的质数，\(x = i \cdot p\)，递推规则共三条：
| 情形 | 递推式 | 依据 |
|---|---|---|
| \(i\) 为质数 | \(\varphi(i) = i - 1\) | 质数与 \([1, i-1]\) 中所有数互质 |
| \(i \bmod p \neq 0\) | \(\varphi(ip) = \varphi(i)\,(p-1)\) | \(\gcd(i, p) = 1\)，用积性与 \(\varphi(p) = p-1\) |
| \(i \bmod p = 0\) | \(\varphi(ip) = \varphi(i) \cdot p\) | \(ip\) 与 \(i\) 质因子集合相同，公式 \(n\prod(1-\frac{1}{p_i})\) 中 \(\prod\) 部分不变，前面的 \(n\) 多乘一个 \(p\) |
14.12.2 模板代码¶
#include <iostream>
using namespace std;
const int MAXN = 10000010;
bool isComp[MAXN];
int primes[MAXN], pcnt;
int phi[MAXN]; // phi[i] = 欧拉函数值
// 线性筛同时求欧拉函数
// 时间复杂度: O(n)
void sievePhi(int n) {
phi[1] = 1;
for (int i = 2; i <= n; i++) {
if (!isComp[i]) {
primes[pcnt++] = i;
phi[i] = i - 1; // 质数 p: phi(p) = p - 1
}
for (int j = 0; j < pcnt && (long long)i * primes[j] <= n; j++) {
int x = i * primes[j];
isComp[x] = true;
if (i % primes[j] == 0) {
// primes[j] 是 i 的最小质因子: x 与 i 质因子集合相同
// 由公式 phi(n) = n * prod(1 - 1/p) 得 phi(x) = phi(i) * p
phi[x] = phi[i] * primes[j];
break;
} else {
// primes[j] 与 i 互质: 由积性 phi(x) = phi(i) * phi(p)
phi[x] = phi[i] * (primes[j] - 1);
}
}
}
}
int main() {
int n;
cin >> n;
sievePhi(n);
for (int i = 1; i <= n; i++) {
cout << "phi(" << i << ") = " << phi[i] << endl;
}
return 0;
}
14.12.3 正确性验证（前 20 项）¶
上面的程序已编译运行，输出的前 20 项如下，并与 \(O(\sqrt{n})\) 单点公式在 \([1, 1000]\) 范围内逐项对拍一致：
| \(n\) | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(\varphi(n)\) | 1 | 1 | 2 | 2 | 4 | 2 | 6 | 4 | 6 | 4 |
| \(n\) | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 |
|---|---|---|---|---|---|---|---|---|---|---|
| \(\varphi(n)\) | 10 | 4 | 12 | 6 | 8 | 8 | 16 | 6 | 18 | 8 |
14.12.4 例题：洛谷 P2158 [SDOI2008] 仪仗队¶
题目链接：https://www.luogu.com.cn/problem/P2158
题意：\(N \times N\) 的方阵（坐标 \((0,0)\) 至 \((N-1, N-1)\)）站满学生，观察者位于 \((0,0)\)。\((x, y)\) 处的学生可见当且仅当连线上没有其他人，即 \(\gcd(x, y) = 1\)。求可见学生总数。
分析：\((1,0)\)、\((0,1)\)、\((1,1)\) 三个方向各可见 1 人；其余可见点关于对角线对称，只需统计 \(x > y \geq 1\) 且 \(\gcd(x, y) = 1\) 的点数，即 \(\sum_{x=2}^{N-1} \varphi(x)\)。故答案为 \(3 + 2\sum_{x=2}^{N-1}\varphi(x)\)（\(N = 1\) 时答案为 0）。
#include <iostream>
using namespace std;
const int MAXN = 40010;
bool isComp[MAXN];
int primes[MAXN], pcnt;
int phi[MAXN];
// 线性筛同时求欧拉函数
void sievePhi(int n) {
phi[1] = 1;
for (int i = 2; i <= n; i++) {
if (!isComp[i]) {
primes[pcnt++] = i;
phi[i] = i - 1;
}
for (int j = 0; j < pcnt && (long long)i * primes[j] <= n; j++) {
int x = i * primes[j];
isComp[x] = true;
if (i % primes[j] == 0) {
phi[x] = phi[i] * primes[j];
break;
} else {
phi[x] = phi[i] * (primes[j] - 1);
}
}
}
}
int main() {
int n;
cin >> n;
if (n == 1) { // 只有自己一个点，看不到任何人
cout << 0 << endl;
return 0;
}
sievePhi(n - 1);
long long sum = 0;
for (int i = 2; i <= n - 1; i++) sum += phi[i];
// 3 = (0,1)、(1,0)、(1,1) 三个方向; 其余可见点按 phi 对称计数
cout << 3 + 2 * sum << endl;
return 0;
}
样例输入 4 输出 9，与题面一致（\(3 + 2(\varphi(2)+\varphi(3)) = 3 + 2 \times 3 = 9\)）；边界 1 输出 0，均已编译验证。
线性筛求 φ 的常见坑
phi[1] = 1必须手动初始化，否则涉及前缀和的题目从第一项就错i % primes[j] == 0分支中要先写phi再break，漏掉break既破坏线性复杂度又会错误覆盖i * primes[j]在int下可能溢出，比较时先转long long（模板已处理）- \(\sum_{i=1}^{n} \varphi(i) \approx \frac{3}{\pi^2} n^2\) 增长很快，前缀和务必用
long long存储
14.13 数论分块（整除分块）¶
14.13.1 问题与核心观察¶
问题：求 \(\sum_{i=1}^{n} \left\lfloor \dfrac{n}{i} \right\rfloor\)。朴素做法 \(O(n)\)，当 \(n \sim 10^{12}\) 时无法通过。
核心观察：\(\lfloor n/i \rfloor\) 至多只有 \(O(\sqrt{n})\) 个不同取值——\(i \leq \sqrt{n}\) 时 \(i\) 本身至多 \(\sqrt{n}\) 个；\(i > \sqrt{n}\) 时 \(\lfloor n/i \rfloor < \sqrt{n}\)，取值也至多 \(\sqrt{n}\) 种。且取值相同的 \(i\) 构成连续区间，可以整块跳过。
块的右端点：设当前块左端点为 \(l\)，记 \(v = \lfloor n/l \rfloor\)，则满足 \(\lfloor n/r \rfloor = v\) 的最大 \(r\) 为
以 \(n = 10\) 为例，分块过程如下：
| 块 \([l, r]\) | \(\lfloor 10/i \rfloor\) | 块内贡献 \((r-l+1) \cdot \lfloor 10/l \rfloor\) |
|---|---|---|
| \([1, 1]\) | 10 | 10 |
| \([2, 2]\) | 5 | 5 |
| \([3, 3]\) | 3 | 3 |
| \([4, 5]\) | 2 | 4 |
| \([6, 10]\) | 1 | 5 |
总和 \(10 + 5 + 3 + 4 + 5 = 27\)，只用 5 块就完成了 10 项求和。
14.13.2 模板代码¶
#include <iostream>
using namespace std;
typedef long long ll;
// 数论分块求 sum_{i=1}^{n} floor(n / i)
// 时间复杂度: O(sqrt(n))
ll divBlockSum(ll n) {
ll ans = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l); // [l, r] 内 floor(n/i) 的值都等于 n/l
ans += (r - l + 1) * (n / l);
}
return ans;
}
int main() {
ll n;
cin >> n;
cout << divBlockSum(n) << endl;
return 0;
}
14.13.3 例题：洛谷 P1403 [AHOI2005] 约数研究¶
题目链接：https://www.luogu.com.cn/problem/P1403
题意：设 \(f(i)\) 表示 \(i\) 的约数个数，求 \(\sum_{i=1}^{n} f(i)\)。
分析：交换统计顺序——不数「每个数有几个约数」，改数「每个数 \(d\) 作为约数出现几次」：\(d\) 是 \(d, 2d, 3d, \ldots\) 的约数，在 \([1, n]\) 内共出现 \(\lfloor n/d \rfloor\) 次。于是
直接套用上面的模板程序即可（原题 \(n \leq 10^6\)，\(O(n)\) 也能过，但数论分块可以应对 \(n \leq 10^{12}\) 的加强版数据）。
已编译验证：输入 3 输出 5（\(f(1)+f(2)+f(3) = 1+2+2\)，与题面样例一致）；输入 10 输出 27，与 14.13.1 的手算分块表吻合。
数论分块常见坑
- 右端点公式必须写成
r = n / (n / l)；循环条件l <= n保证了n / l非零，不会除零 - 贡献
(r - l + 1) * (n / l)两个因子都可能很大，务必用long long；带模意义下每项先取模再累加 - 和式带权时（如 \(\sum_i i \cdot \lfloor n/i \rfloor\)），块内权值和用等差数列求和 \(\frac{(l+r)(r-l+1)}{2}\) 计算
- 二维形式 \(\sum_i \lfloor n/i \rfloor \lfloor m/i \rfloor\) 的右端点取 \(\min\left(\lfloor n/\lfloor n/l \rfloor\rfloor,\ \lfloor m/\lfloor m/l \rfloor\rfloor\right)\)
14.14 BSGS（大步小步算法）¶
14.14.1 离散对数问题¶
问题：给定质数 \(p\) 与整数 \(a, b\)（\(\gcd(a, p) = 1\)），求最小非负整数 \(x\) 使得
这相当于模意义下的「取对数」，称为离散对数问题。由费马小定理，\(a^x \bmod p\) 以 \(p-1\) 为周期，因此若有解，必存在 \(x < p\)，只需在 \([0, p)\) 内搜索。
大步小步（Baby-Step Giant-Step）：设 \(t = \lceil \sqrt{p} \rceil\)，把 \(x\) 写成 \(x = i \cdot t - j\)（\(1 \leq i \leq t\)，\(0 \leq j < t\)），原方程变形为
- 小步：枚举 \(j \in [0, t)\)，把 \(b \cdot a^j \bmod p\) 存入哈希表（键为值、值为 \(j\)）
- 大步：枚举 \(i \in [1, t]\)，计算 \((a^t)^i \bmod p\) 并查表，命中即得 \(x = it - j\)
graph LR
subgraph BSGS 求解 a^x ≡ b mod p
A["t = ⌈√p⌉, 设 x = i·t − j"] --> B["小步: 存 b·a^j → j, 0 ≤ j < t"]
B --> C["大步: 枚举 (a^t)^i, 1 ≤ i ≤ t"]
C --> D{"命中哈希表?"}
D -->|是| E["x = i·t − j"]
D -->|否| F["无解"]
end14.14.2 模板代码¶
#include <iostream>
#include <cmath>
#include <unordered_map>
using namespace std;
typedef long long ll;
// 快速幂
ll qpow(ll a, ll b, ll m) {
ll res = 1 % m;
a %= m;
while (b > 0) {
if (b & 1) res = res * a % m;
a = a * a % m;
b >>= 1;
}
return res;
}
// BSGS: 求 a^x ≡ b (mod p) 的最小非负整数解
// 要求 p 为质数且 gcd(a, p) = 1, 无解返回 -1
// 时间复杂度: O(sqrt(p))
ll bsgs(ll a, ll b, ll p) {
a %= p; b %= p;
if (b == 1 % p) return 0; // a^0 = 1
ll t = (ll)ceil(sqrt((double)p));
// Baby step: 存 b * a^j -> j (0 <= j < t)
// 相同键后写入的 j 更大, 对应的解 x = i*t - j 更小
unordered_map<ll, ll> table;
ll cur = b;
for (ll j = 0; j < t; j++) {
table[cur] = j;
cur = cur * a % p;
}
// Giant step: 依次检查 a^(i*t) 是否命中表中的 b * a^j
ll at = qpow(a, t, p); // a^t
cur = at;
for (ll i = 1; i <= t; i++) {
auto it = table.find(cur);
if (it != table.end()) return i * t - it->second; // a^(i*t) = b * a^j
cur = cur * at % p;
}
return -1;
}
int main() {
ll p, b, n;
cin >> p >> b >> n; // 求最小的 l 使 b^l ≡ n (mod p)
ll ans = bsgs(b, n, p);
if (ans == -1) cout << "no solution" << endl;
else cout << ans << endl;
return 0;
}
14.14.3 例题：洛谷 P3846 [TJOI2007] 可爱的质数/【模板】BSGS¶
题目链接：https://www.luogu.com.cn/problem/P3846
题意：给定质数 \(p\) 与整数 \(b, n\)（\(2 \leq b, n < p < 2^{31}\)），求最小的非负整数 \(l\) 使 \(b^l \equiv n \pmod{p}\)；无解输出 no solution。
上面的模板程序就是本题的完整题解。
验证记录（点击展开）
以上程序已编译运行验证：
- 官方样例
5 2 3：输出3（\(2^3 = 8 \equiv 3 \pmod 5\)） 7 3 4：输出4（\(3^4 = 81 \equiv 4 \pmod 7\)）7 2 3：输出no solution（2 模 7 的幂只有 \(\{1, 2, 4\}\)）- 大数据
998244353 3 20190810：输出732811431，用 Pythonpow(3, 732811431, 998244353)反验证结果为20190810，正确
BSGS 常见坑
- \(b = 1\) 时答案为 \(x = 0\)，要在建表前特判（模板中
b == 1 % p的写法同时兼容 \(p = 1\)） - \(a \equiv 0 \pmod{p}\) 时不满足 \(\gcd(a, p) = 1\)：若 \(b \equiv 0\) 则 \(x = 1\)，否则无解，需单独特判
- 哈希表中相同的键要保留较大的 \(j\)（循环中直接赋值覆盖即可），这样求出的 \(x = it - j\) 才是最小解
unordered_map最坏可被构造数据卡到 \(O(n)\) 单次操作，比赛中可自定义哈希函数或改用「排序 + 二分」- 模数不是质数、或 \(\gcd(a, p) \neq 1\) 时需要扩展 BSGS（exBSGS），本节模板不适用
14.15 ACM Day9 数论 练习题¶
本节为 ACM Day9 数论专题的配套练习。题目覆盖 GCD 性质、组合数学、快速幂、质因数分解、扩展欧几里得、扩展 CRT 等核心知识点，难度由易到难分级排列。建议先独立思考 15-20 分钟再看解题提示。
使用指南
每道题的知识点标注对应本章各小节，可快速回溯相关理论。解题提示仅给出思路方向，具体实现需自行完成。完成所有练习后，尝试用本章 14.10 的算法选择指南总结每道题所属的算法类别。
题目总表¶
| 题号 | 平台 | 题目 | 难度 | 链接 | 完成 |
|---|---|---|---|---|---|
| 1920C | CF | Partitioning the Array | 1600 | https://codeforces.com/contest/1920/problem/C | - [ ] |
| 1957C | CF | How Does the Rook Move? | 1600 | https://codeforces.com/contest/1957/problem/C | - [ ] |
| 1985G | CF | D-Function | 1700 | https://codeforces.com/contest/1985/problem/G | - [ ] |
| 2039D | CF | Shohag Loves GCD | 1700 | https://codeforces.com/contest/2039/problem/D | - [ ] |
| 2050F | CF | Maximum modulo equality | 1700 | https://codeforces.com/contest/2050/problem/F | - [ ] |
| 2065G | CF | Skibidus and Capping | 1700 | https://codeforces.com/contest/2065/problem/G | - [ ] |
| 1995C | CF | Squaring | 1800 | https://codeforces.com/contest/1995/problem/C | - [ ] |
| abc240_d | AT | Strange Balls | 1600 | https://atcoder.jp/contests/abc240/tasks/abc240_d | - [ ] |
| abc303_d | AT | Shift vs. CapsLock | 1700 | https://atcoder.jp/contests/abc303/tasks/abc303_d | - [ ] |
| abc286_c | AT | Rotate and Palindrome | 1800 | https://atcoder.jp/contests/abc286/tasks/abc286_c | - [ ] |
| P1082 | 洛谷 | 【模板】同余方程 | 1700 | https://www.luogu.com.cn/problem/P1082 | - [ ] |
| P4777 | 洛谷 | 【模板】扩展中国剩余定理 | 1800 | https://www.luogu.com.cn/problem/P4777 | - [ ] |
知识点分布总览¶
| 知识点 | 对应章节 | 涉及题号 |
|---|---|---|
| GCD 性质 | 14.3 | 第 1、5、7 题 |
| 组合数学 | 14.6 | 第 2 题 |
| 快速幂 | 14.4 | 第 3 题 |
| 栈与模拟 | 基础数据结构 | 第 4 题 |
| 质因数分解 | 14.2 | 第 6 题 |
| DP 与分类讨论 | 动态规划 | 第 8 题 |
| 枚举与回文 | 字符串算法 | 第 9 题 |
| 容斥/博弈 | 14.8 / 14.9 | 第 10 题 |
| 扩展欧几里得 | 14.3 | 第 11 题 |
| 扩展 CRT | 14.7 | 第 12 题 |
一、基础难度（建议 30 分钟以内完成）¶
学习建议
基础题旨在巩固课堂所学的基本算法模板。做题时注意手写代码而非复制模板，重点理解每一步的推导逻辑。完成后尝试用不同的方法重新求解，加深对算法本质的理解。这些题目对应本章中 GCD（14.3）、组合数学（14.6）、快速幂（14.4）等基础节的内容。
| 序号 | 题号 | 题目名/简述 | 知识点 | 难度 | 解题提示 |
|---|---|---|---|---|---|
| 1 | CF 1920C | Partitioning the Array：给定数组，求满足特定条件的最大划分方案 | GCD 性质 | ⭐⭐ | 枚举划分块数 \(k\)（\(k\) 需满足 \(k \mid n\)），对每个 \(k\) 检查每块中对应位置差值的 GCD 是否大于 1。由于 \(n\) 的因子数量有限，总复杂度可控。关键观察：若所有块的对应位置差值的 \(\gcd > 1\)，则划分合法。 |
| 2 | CF 1957C | How Does the Rook Move?：棋盘上放置车，求方案数 | 组合数学 | ⭐⭐ | 每步操作相当于固定一行一列，剩下未被占的行列数递减 1 或 2（取决于车是否放在对角线上）。设剩余 \(m\) 行/列，递推 \(f(m) = f(m-1) + (m-1) \cdot f(m-2)\)。预处理阶乘和组合数辅助计算最终方案数。 |
| 3 | CF 1995C | Squaring：对数组元素反复开平方使其满足单调性 | 快速幂 | ⭐⭐ | 注意到 \(a^{2^k}\) 增长极快，当指数超过 30 左右时数值已远超题目范围。从前往后维护当前元素所需的最小平方次数，累加每步差值。若前缀全为 1 则后续无法增大，需特判无解。 |
| 4 | AtCoder abc240_d | Strange Balls：球入桶的消去模拟 | 栈与模拟 | ⭐⭐ | 使用栈维护连续相同球的段，每段记录(编号, 数量)。当某段数量恰好等于其编号时整段弹出。每次操作后输出栈中所有段的数量之和。注意边界：栈为空时直接压入新段。 |
二、中等难度（建议 45-60 分钟完成）¶
学习建议
中等题通常需要将多个知识点组合运用，或在基本算法的基础上做适当的变形。解题时先分析问题本质，识别可以套用的数学模型，再考虑边界情况和复杂度优化。这些题目综合了 14.2（质因数分解）、14.3（GCD 性质）、14.6（组合数学）以及基础 DP 的知识。
| 序号 | 题号 | 题目名/简述 | 知识点 | 难度 | 解题提示 |
|---|---|---|---|---|---|
| 5 | CF 1985G | D-Function：统计满足数字之和为 \(d\) 的倍数的区间 | GCD 与前缀和 | ⭐⭐⭐ | 数字之和函数 \(D(x)\) 满足 \(D(kx) = D(x)\)。利用前缀和差分思想，统计满足 \(D(l) \equiv D(r) \pmod{k}\) 的区间对数，再乘上 \(l\) 的可选方案数。对每个 \(k\) 独立计算后求和。注意模运算下的取值范围。 |
| 6 | CF 2039D | Shohag Loves GCD：从给定集合选元素构造数组，使 \(a_{\gcd(i,j)} \neq \gcd(a_i, a_j)\) 对所有 \(i < j\) 成立且总和最大 | 质因数分解 | ⭐⭐⭐ | 设 \(\Omega(i)\) 为 \(i\) 的质因子个数（含重数），可用线性筛（14.2.2 的最小质因子数组）预处理。将集合降序排序，令 \(a_i\) 取集合中第 \(\Omega(i)+1\) 大的元素：若 \(i \mid j\) 且 \(i \neq j\)，则 \(\Omega(i) < \Omega(j)\)，从而 \(a_i > a_j\)，约束自然满足。所需层数为 \(\lfloor \log_2 n \rfloor + 1\)，集合大小不足时输出 -1。 |
| 7 | CF 2050F | Maximum Modulo Equality：求区间内最大公共模数 | GCD 性质 | ⭐⭐⭐ | 关键观察：若区间内所有元素模 \(m\) 同余，则 \(m\) 必须整除任意相邻元素之差。因此答案等于区间差分数组的 GCD。预处理差分数组的 ST 表或 Sparse Table，\(O(1)\) 查询区间 GCD。特判区间长度为 1 时答案为 \(10^9\)。 |
| 8 | AtCoder abc303_d | Shift vs. CapsLock：大小写切换与移位的最优策略 | DP 与分类讨论 | ⭐⭐⭐ | 设 \(dp[i][0/1]\) 表示处理完第 \(i\) 个字符后 CapsLock 状态为关/开的最小代价。转移时需分类讨论：当前字符是否匹配当前大小写状态，以及是否按下 Shift 键。注意 Shift 和 CapsLock 代价的不对称性。 |
| 9 | AtCoder abc286_c | Rotate and Palindrome：旋转字符串后求最小回文修改代价 | 枚举与回文 | ⭐⭐⭐ | 枚举旋转偏移量 \(k = 0 \sim n-1\)，对每种旋转计算使字符串变为回文的代价。代价计算方法：比较对称位置 \((i, n-1-i)\)，统计不匹配数乘以替换代价 \(A\) 或交换代价 \(B\)，取较小值。总时间 \(O(n^2)\)。 |
三、进阶难度（建议 60-90 分钟完成）¶
学习建议
进阶题往往涉及算法的深层性质或多个算法的综合运用。如果卡住，建议回顾本章 14.3（扩展欧几里得）和 14.7（扩展 CRT）的推导过程，理解合并方程的核心思想。这类题在正式比赛中是区分选手实力的关键，务必反复练习直至熟练。
| 序号 | 题号 | 题目名/简述 | 知识点 | 难度 | 解题提示 |
|---|---|---|---|---|---|
| 10 | CF 2065G | Skibidus and Capping：与整数约束相关的容斥/博弈问题 | 容斥/博弈 | ⭐⭐⭐⭐ | 仔细分析题目中的约束条件，将其转化为集合计数或博弈判定问题。用容斥原理排除不合法的情况，或通过 SG 函数分析必胜必败态。参考本章 14.8（容斥原理）和 14.9（博弈论）的内容，注意边界条件和大数取模。 |
| 11 | 洛谷 P1082 | 【模板】同余方程：求 \(ax \equiv b \pmod{m}\) 的最小正整数解 | 扩展欧几里得 | ⭐⭐⭐⭐ | 直接套用扩展欧几里得模板求解。步骤：(1) 用 exgcd 求 \(\gcd(a, m)\)；(2) 若 \(\gcd(a, m) \nmid b\) 则方程无解；(3) 否则特解 \(x_0 = x \cdot (b / g)\)，通解为 \(x_0 + k \cdot (m/g)\)；(4) 对 \(m/g\) 取模调整到 \([1, m/g]\) 范围内得到最小正整数解。 |
| 12 | 洛谷 P4777 | 【模板】扩展中国剩余定理：模数不互质的同余方程组 | 扩展 CRT | ⭐⭐⭐⭐ | 逐个合并方程：将 \(x \equiv a_1 \pmod{m_1}\) 与 \(x \equiv a_2 \pmod{m_2}\) 合并。具体做法：设 \(x = m_1 t + a_1\)，代入第二个方程得 \(m_1 t \equiv a_2 - a_1 \pmod{m_2}\)，用 exgcd 解出 \(t\)，新模数为 \(\text{lcm}(m_1, m_2)\)。注意溢出问题，乘法运算建议使用 __int128 或快速乘避免爆 long long（见 14.7.1 的 warning）。 |
四、练习策略¶
高效刷题建议
- 分层推进：先完成基础题建立信心，再攻克中等和进阶题。每层完成后进行总结回顾
- 限时训练：基础题每道 20 分钟，中等题每道 40 分钟，进阶题每道 60 分钟。超时后先看提示再思考，不要死磕
- 总结归纳：每做完一道题，回顾使用了哪些知识点，在笔记本上记录关键技巧和易错点
- 代码模板化：将 exgcd、快速幂、组合数预处理等常用代码整理为模板，但务必理解每一行的含义而非死记硬背
- 复盘错题：做错或超时的题目标记为"待复习"，一周后重新做一遍检验掌握程度
推荐练习顺序¶
循序渐进路线
第一轮（基础巩固）：第 1 题 → 第 3 题 → 第 4 题 → 第 2 题 这四道题分别练习 GCD 性质、快速幂、栈模拟和组合数学，覆盖本章最基础的四个知识点。建议在 2 小时内全部完成。
第二轮（综合运用）：第 7 题 → 第 5 题 → 第 8 题 → 第 9 题 → 第 6 题 中等难度题目需要将基础算法与问题建模能力结合。重点体会 GCD 性质的灵活运用（第 5、7 题）和 DP 状态设计的思路（第 8 题）。
第三轮（深度突破）：第 11 题 → 第 12 题 → 第 10 题 进阶题对应本章的核心难点。先练习扩展欧几里得的模板题（P1082），再挑战扩展 CRT（P4777），最后综合运用容斥与博弈解决第 10 题。
错误排查清单¶
| 症状 | 可能原因 | 对应解法 |
|---|---|---|
| 扩展欧几里得答案为负 | 未做取模调整 | 解出后对 \(m/\gcd\) 取模并补正 |
| 快速幂结果为 0 | 取模后 base 变为 0 | 检查 base 是否被模数整除 |
| 扩展 CRT 合并后溢出 | 乘法爆 long long |
使用 __int128 或 (a*b)%m 快速乘 |
| 组合数结果为 0 | 预处理范围不足 | 确保阶乘数组大小覆盖最大查询值 |
| GCD 递归栈溢出 | 参数过大或为 0 | 使用 __gcd 或改为迭代写法 |
常见易错点
- 溢出问题：GCD、快速幂、扩展 CRT 中的乘法运算容易溢出
long long范围，必要时使用__int128或快速乘 - 取模细节：负数取模需调整为正数，公式中先除后乘避免精度问题；注意模数是否为质数决定了能否使用费马小定理
- 边界特判：\(\gcd(a, 0) = a\)、\(0^0 = 1\)（竞赛中通常如此约定）、模数为 1 等特殊情况不要遗漏
- 时间复杂度：注意算法在最坏情况下的复杂度是否满足题目限制，必要时做常数优化或更换算法思路