数论与数学基础
第三十五章：数论与数学基础
你写了一个加密算法，需要计算 2^1000000007 mod 998244353。如果用循环乘一百万次，程序跑到宇宙热寂也跑不完。但用快速幂，50 次乘法就够了。你需要判断一个数是否为素数——对于 10^18 量级的数，试除法需要 10^9 次除法，而 Miller-Rabin 算法只需要几十次模幂运算。
数论不是数学家的玩具。它是密码学的基石（RSA 依赖大素数分解的困难性）、哈希函数的理论保证（取模运算的均匀分布性）、竞赛算法的核心工具（组合计数离不开模逆元）。本章从最基础的 GCD 出发，逐步构建一个完整的数论工具箱。
Level 1 · 你需要知道的
1.1 最大公约数（GCD）与欧几里得算法
问题定义：给定两个正整数 a 和 b，求它们的最大公约数 gcd(a, b)。
核心定理：gcd(a, b) = gcd(b, a mod b)，当 b = 0 时，gcd(a, 0) = a。
这个定理为什么成立？设 d = gcd(a, b)，则 d|a 且 d|b。因为 a mod b = a - ⌊a/b⌋·b，而 d|a 且 d|b，所以 d|(a mod b)。反方向同理，d|b 且 d|(a mod b)，则 d|(a mod b + ⌊a/b⌋·b) = d|a。因此 gcd(a, b) 的公约数集合与 gcd(b, a mod b) 的公约数集合完全相同，最大公约数自然相等。
def gcd(a: int, b: int) -> int:
"""欧几里得算法求最大公约数
时间复杂度: O(log(min(a, b)))
空间复杂度: O(1) 迭代版本
"""
while b:
a, b = b, a % b
return a
# 递归版本（更直观但有栈溢出风险）
def gcd_recursive(a: int, b: int) -> int:
if b == 0:
return a
return gcd_recursive(b, a % b)
复杂度分析：欧几里得算法的迭代次数不超过 O(log(min(a, b)))。这是因为每两步迭代后，较大的数至少减半。具体地，如果 a > b，则 a mod b < a/2（证明：若 b ≤ a/2，则 a mod b < b ≤ a/2；若 b > a/2，则 a mod b = a - b < a/2）。因此最多经过 2·log₂(min(a, b)) 次迭代即可终止。
最坏情况：当 a 和 b 是相邻的斐波那契数时，欧几里得算法达到最坏情况。例如 gcd(F_{n+1}, F_n) 需要 n-1 次迭代，因为 F_{n+1} mod F_n = F_{n-1}。这也从侧面证明了复杂度是 O(log n)，因为 F_n ≈ φ^n / √5，所以 n ≈ log_φ(F_n)。
1.2 最小公倍数（LCM）
定义：lcm(a, b) 是能同时被 a 和 b 整除的最小正整数。
与 GCD 的关系：lcm(a, b) = a · b / gcd(a, b)
为什么这个公式成立？将 a 和 b 做质因数分解：a = p₁^α₁ · p₂^α₂ · ... ，b = p₁^β₁ · p₂^β₂ · ...。那么 gcd(a, b) = ∏ pᵢ^min(αᵢ, βᵢ)，lcm(a, b) = ∏ pᵢ^max(αᵢ, βᵢ)。由于 min(x, y) + max(x, y) = x + y，所以 gcd(a, b) · lcm(a, b) = a · b。
def lcm(a: int, b: int) -> int:
"""最小公倍数
注意：先除后乘，防止溢出（在有限精度语言中）
"""
return a // gcd(a, b) * b
# 多个数的 LCM
def lcm_multiple(numbers: list) -> int:
from functools import reduce
return reduce(lcm, numbers)
常见错误：直接计算 a * b / gcd(a, b) 可能导致整数溢出（在 C/C++/Java 中）。正确做法是先除后乘：a / gcd(a, b) * b。Python 的大整数不存在溢出问题，但在竞赛中使用 C++ 时必须注意。
1.3 素数判断
试除法：一个大于 1 的自然数 n，如果不能被 2 到 √n 之间的任何整数整除，则 n 是素数。
为什么只需要检查到 √n？如果 n = a · b 且 a ≤ b，则 a ≤ √n。所以如果 n 有一个大于 1 小于 n 的因子，那么必有一个因子不超过 √n。
def is_prime(n: int) -> bool:
"""试除法判断素数
时间复杂度: O(√n)
"""
if n < 2:
return False
if n < 4:
return True
if n % 2 == 0 or n % 3 == 0:
return False
# 所有大于 3 的素数都可以写成 6k±1 的形式
i = 5
while i * i <= n:
if n % i == 0 or n % (i + 2) == 0:
return False
i += 6
return True
6k±1 优化的原理：所有整数可以写成 6k, 6k+1, 6k+2, 6k+3, 6k+4, 6k+5 之一。其中 6k, 6k+2, 6k+4 能被 2 整除，6k+3 能被 3 整除。所以大于 3 的素数只可能是 6k+1 或 6k+5（即 6k-1）的形式。这让我们跳过了 2/3 的候选因子。
1.4 埃拉托斯特尼筛法（Sieve of Eratosthenes）
问题：找出 1 到 n 之间的所有素数。
思想：从 2 开始，把每个素数的所有倍数标记为合数。未被标记的就是素数。
def sieve_of_eratosthenes(n: int) -> list:
"""埃氏筛法
时间复杂度: O(n log log n)
空间复杂度: O(n)
"""
if n < 2:
return []
is_prime = [True] * (n + 1)
is_prime[0] = is_prime[1] = False
# 只需筛到 √n
i = 2
while i * i <= n:
if is_prime[i]:
# 从 i² 开始标记（i*2, i*3, ..., i*(i-1) 已经被更小的素数标记过了）
for j in range(i * i, n + 1, i):
is_prime[j] = False
i += 1
return [i for i in range(2, n + 1) if is_prime[i]]
为什么从 i² 开始标记？当我们处理素数 p 时，p·2, p·3, ..., p·(p-1) 这些数已经在处理 2, 3, ..., p-1 时被标记过了。例如，当 p=5 时，5·2=10（已被2标记），5·3=15（已被3标记），5·4=20（已被2标记）。所以第一个需要新标记的是 5·5=25。
复杂度分析：总操作次数为 ∑{p≤n, p是素数} n/p = n · ∑{p≤n} 1/p。由 Mertens 定理（Mertens, 1874），素数的倒数之和 ∑_{p≤n} 1/p = ln(ln n) + M（M 是 Meissel-Mertens 常数 ≈ 0.2615）。因此总复杂度为 O(n log log n)。
1.5 线性筛（欧拉筛）
埃氏筛的问题是一个合数可能被多次标记。例如 30 会被 2, 3, 5 各标记一次。线性筛保证每个合数只被其最小质因子标记一次，从而实现严格 O(n) 的时间复杂度。
def linear_sieve(n: int) -> list:
"""线性筛（欧拉筛）
时间复杂度: O(n)
空间复杂度: O(n)
核心思想: 每个合数 m 只被它的最小质因子 lpf(m) 标记一次
"""
is_prime = [True] * (n + 1)
primes = []
for i in range(2, n + 1):
if is_prime[i]:
primes.append(i)
for p in primes:
if i * p > n:
break
is_prime[i * p] = False
# 关键: 当 p | i 时停止
# 因为后续 i * p' (p' > p) 的最小质因子是 p（因为 p|i），
# 应该留给 (i*p'/p) * p 来标记
if i % p == 0:
break
return primes
为什么 if i % p == 0: break 保证线性？假设 i = p · k，那么对于 primes 列表中 p 之后的任何素数 q（q > p），i · q = p · k · q。这个数的最小质因子是 p（因为 p|i，所以 p|i·q），而不是 q。如果我们继续标记 i · q，就会导致这个数被标记两次——一次是现在（错误地认为 q 是最小质因子），一次是未来遍历到 k · q 时（正确地由 p 标记）。所以我们必须在 p|i 时停止。
线性筛 vs 埃氏筛的实际选择：
| 特性 | 埃氏筛 | 线性筛 |
|---|---|---|
| 时间复杂度 | O(n log log n) | O(n) |
| 常数因子 | 小（缓存友好） | 大（分支预测差） |
| 实际速度（n=10⁸） | 约相同 | 约相同 |
| 附加信息 | 无 | 可得最小质因子 |
| 代码复杂度 | 简单 | 稍复杂 |
实践中对于 n ≤ 10⁸，两者速度接近。线性筛的优势在于它同时给出了每个数的最小质因子，这对于快速分解质因数非常有用。
1.6 常见错误与陷阱
错误 1：GCD 的参数顺序
# 错误：认为第一个参数必须大于第二个
gcd(3, 12) # 正确！第一步变成 gcd(12 % 3... 不对)
# 实际上 gcd(3, 12) → gcd(12, 3) 是错的
# 正确过程: gcd(3, 12) → a=3, b=12 → a,b = 12, 3%12=3... 不对
# 让我重新算: gcd(3, 12): a=3, b=12 → a,b = b, a%b = 12, 3
# 等等，不对。a=3, b=12, 新的 a=b=12, 新的 b=a%b=3%12=3
# gcd(12, 3) → a=12, b=3 → a,b = 3, 0 → return 3
# 其实无论顺序如何，第一步就会自动交换。gcd(3,12) 第一步得到 gcd(12, 3)
欧几里得算法不需要保证 a > b。如果 a < b，第一步 a % b = a（因为 a < b），所以变成 gcd(b, a)，自动交换了顺序。
错误 2：LCM 溢出
# C++ 中的危险写法（Python 无此问题）
# long long lcm = a * b / gcd(a, b); // a*b 可能溢出！
# 正确写法
# long long lcm = a / gcd(a, b) * b; // 先除后乘
错误 3：筛法中的数组越界
# 错误：range 应该到 n+1
is_prime = [True] * n # 应该是 n+1
# 访问 is_prime[n] 会越界
错误 4：把 1 当作素数
1 既不是素数也不是合数。这是约定俗成的定义，原因是唯一分解定理要求分解的唯一性——如果 1 是素数，那么 6 = 2 × 3 = 1 × 2 × 3 = 1 × 1 × 2 × 3...就不唯一了。
Level 2 · 它是怎么运行的
2.1 快速幂取模（Modular Exponentiation）
问题：计算 a^n mod m，其中 n 可能非常大（如 10^18）。
思想：利用二进制展开，将 n 次乘法减少到 O(log n) 次。
a^13 = a^(1101₂) = a^8 · a^4 · a^1
每次平方使指数翻倍，同时检查 n 的当前最低位是否为 1。
def power_mod(base: int, exp: int, mod: int) -> int:
"""快速幂取模
计算 base^exp % mod
时间复杂度: O(log exp)
空间复杂度: O(1)
"""
result = 1
base %= mod
while exp > 0:
# 如果当前位为 1，乘上当前的 base
if exp & 1:
result = result * base % mod
# base 自乘，对应指数翻倍
base = base * base % mod
exp >>= 1
return result
执行过程示例：计算 2^13 mod 1000
| 步骤 | exp (二进制) | base | result | 操作 |
|---|---|---|---|---|
| 初始 | 1101 | 2 | 1 | - |
| 1 | 110 | 4 | 2 | exp&1=1, result*=base; base*=base |
| 2 | 11 | 16 | 2 | exp&1=0; base*=base |
| 3 | 1 | 256 | 32 | exp&1=1, result*=base; base*=base |
| 4 | 0 | 65536 | 8192 | exp&1=1, result*=base; 结束 |
验证：2^13 = 8192，8192 mod 1000 = 192。等等，让我重新算：result 在第3步应该是 216=32，第4步 32256=8192。8192 mod 1000 = 192。正确。
为什么每步都取模？如果不取模，中间结果会指数增长——2^(10^18) 的位数比宇宙中的原子还多。每步取模保证中间结果始终小于 mod²（两个小于 mod 的数相乘），这样才能在有限精度内完成计算。
2.2 扩展欧几里得算法（Extended Euclidean Algorithm）
问题：给定 a, b，找整数 x, y 使得 ax + by = gcd(a, b)。
贝祖定理（Bézout's Identity）：对任意整数 a, b，存在整数 x, y 使得 ax + by = gcd(a, b)。这是扩展欧几里得算法的理论基础。
def extended_gcd(a: int, b: int) -> tuple:
"""扩展欧几里得算法
返回 (g, x, y) 使得 ax + by = g = gcd(a, b)
时间复杂度: O(log(min(a, b)))
"""
if b == 0:
return a, 1, 0 # a*1 + 0*0 = a
g, x1, y1 = extended_gcd(b, a % b)
# 由 b*x1 + (a%b)*y1 = g
# 即 b*x1 + (a - ⌊a/b⌋*b)*y1 = g
# 即 a*y1 + b*(x1 - ⌊a/b⌋*y1) = g
x = y1
y = x1 - (a // b) * y1
return g, x, y
推导过程详解：
设我们已经知道 gcd(b, a%b) = g，且 b·x₁ + (a%b)·y₁ = g。
由于 a%b = a - ⌊a/b⌋·b，代入得：
- b·x₁ + (a - ⌊a/b⌋·b)·y₁ = g
- b·x₁ + a·y₁ - ⌊a/b⌋·b·y₁ = g
- a·y₁ + b·(x₁ - ⌊a/b⌋·y₁) = g
对比 a·x + b·y = g，得 x = y₁，y = x₁ - ⌊a/b⌋·y₁。
迭代版本（避免递归栈溢出）：
def extended_gcd_iterative(a: int, b: int) -> tuple:
"""扩展欧几里得算法 - 迭代版本"""
old_r, r = a, b
old_s, s = 1, 0
old_t, t = 0, 1
while r != 0:
quotient = old_r // r
old_r, r = r, old_r - quotient * r
old_s, s = s, old_s - quotient * s
old_t, t = t, old_t - quotient * t
# old_r = gcd, old_s = x, old_t = y
return old_r, old_s, old_t
2.3 模逆元（Modular Inverse）
定义：a 的模 m 逆元是满足 a·x ≡ 1 (mod m) 的 x，记作 a⁻¹ mod m。
存在条件：a 的模 m 逆元存在当且仅当 gcd(a, m) = 1（即 a 与 m 互质）。
方法一：扩展欧几里得
由 ax + my = gcd(a, m) = 1，得 ax ≡ 1 (mod m)，所以 x 就是 a 的模 m 逆元。
def mod_inverse(a: int, m: int) -> int:
"""用扩展欧几里得求模逆元
前提: gcd(a, m) = 1
返回 x 使得 a*x ≡ 1 (mod m)
"""
g, x, _ = extended_gcd(a, m)
if g != 1:
raise ValueError(f"{a} 和 {m} 不互质，逆元不存在")
return x % m # 确保返回正数
方法二：费马小定理（当 m 是素数时）
费马小定理：若 p 是素数且 gcd(a, p) = 1，则 a^(p-1) ≡ 1 (mod p)。
因此 a · a^(p-2) ≡ 1 (mod p)，即 a⁻¹ ≡ a^(p-2) (mod p)。
def mod_inverse_fermat(a: int, p: int) -> int:
"""用费马小定理求模逆元（仅当 p 是素数时有效）
a^(-1) ≡ a^(p-2) (mod p)
"""
return power_mod(a, p - 2, p)
两种方法的比较：
| 方法 | 适用条件 | 时间复杂度 | 优势 |
|---|---|---|---|
| 扩展欧几里得 | gcd(a, m) = 1 | O(log m) | 通用，m 不必是素数 |
| 费马小定理 | m 是素数 | O(log m) | 代码简洁，常数小 |
2.4 组合数计算与 Lucas 定理
组合数 C(n, k) 的计算：
方法一（小规模）：使用杨辉三角 / Pascal's Triangle
def comb_pascal(n: int, k: int) -> int:
"""杨辉三角法计算组合数
时间复杂度: O(nk)
空间复杂度: O(k)
"""
if k > n - k:
k = n - k # C(n,k) = C(n,n-k), 取较小的 k
dp = [0] * (k + 1)
dp[0] = 1
for i in range(1, n + 1):
for j in range(min(i, k), 0, -1):
dp[j] = dp[j] + dp[j - 1]
return dp[k]
方法二（模素数）：预处理阶乘及其逆元
class CombMod:
"""支持 O(1) 查询组合数 mod p 的预处理类
预处理时间: O(n)
查询时间: O(1)
"""
def __init__(self, n: int, mod: int):
self.mod = mod
self.fact = [1] * (n + 1)
self.inv_fact = [1] * (n + 1)
# 预处理阶乘
for i in range(1, n + 1):
self.fact[i] = self.fact[i - 1] * i % mod
# 预处理阶乘逆元
# 先求 n! 的逆元，然后倒推
self.inv_fact[n] = power_mod(self.fact[n], mod - 2, mod)
for i in range(n - 1, -1, -1):
self.inv_fact[i] = self.inv_fact[i + 1] * (i + 1) % mod
def comb(self, n: int, k: int) -> int:
"""计算 C(n, k) mod p"""
if k < 0 or k > n:
return 0
return self.fact[n] * self.inv_fact[k] % self.mod * self.inv_fact[n - k] % self.mod
Lucas 定理：当 p 是素数且 n, k 可能很大（但 p 较小）时：
C(n, k) mod p = C(n mod p, k mod p) · C(n // p, k // p) mod p
def lucas(n: int, k: int, p: int) -> int:
"""Lucas 定理计算 C(n, k) mod p
要求: p 是素数
适用: n, k 很大但 p 较小（如 p ≤ 10⁵）
时间复杂度: O(p + log_p(n) · log p)
"""
# 预处理 0..p-1 的阶乘和逆阶乘
fact = [1] * p
for i in range(1, p):
fact[i] = fact[i - 1] * i % p
inv_fact = [1] * p
inv_fact[p - 1] = power_mod(fact[p - 1], p - 2, p)
for i in range(p - 2, -1, -1):
inv_fact[i] = inv_fact[i + 1] * (i + 1) % p
def small_comb(n, k):
if k < 0 or k > n:
return 0
return fact[n] * inv_fact[k] % p * inv_fact[n - k] % p
result = 1
while n > 0 or k > 0:
ni, ki = n % p, k % p
result = result * small_comb(ni, ki) % p
n //= p
k //= p
return result
Lucas 定理的直觉：把 n 和 k 看作 p 进制数，C(n, k) mod p 等于每一位上的组合数之积。这类似于十进制下"逐位"处理的思想，但适用于组合数在模素数下的计算。
2.5 快速幂的变体与应用
矩阵快速幂：用于加速线性递推关系。
def matrix_mult(A: list, B: list, mod: int) -> list:
"""矩阵乘法 mod"""
n = len(A)
m = len(B[0])
k = len(B)
C = [[0] * m for _ in range(n)]
for i in range(n):
for j in range(m):
for l in range(k):
C[i][j] = (C[i][j] + A[i][l] * B[l][j]) % mod
return C
def matrix_power(M: list, n: int, mod: int) -> list:
"""矩阵快速幂
计算 M^n mod mod
用于加速 Fibonacci 等线性递推
"""
size = len(M)
# 单位矩阵
result = [[1 if i == j else 0 for j in range(size)] for i in range(size)]
while n > 0:
if n & 1:
result = matrix_mult(result, M, mod)
M = matrix_mult(M, M, mod)
n >>= 1
return result
def fibonacci_fast(n: int, mod: int = 10**9 + 7) -> int:
"""O(log n) 计算第 n 个 Fibonacci 数
[F(n+1), F(n)] = [[1,1],[1,0]]^n · [1, 0]
"""
if n <= 1:
return n
M = [[1, 1], [1, 0]]
result = matrix_power(M, n, mod)
return result[0][1]
应用场景：
- Fibonacci 数列第 10^18 项（mod p）
- 线性递推 f(n) = a₁f(n-1) + a₂f(n-2) + ... + aₖf(n-k) 的快速求解
- 图上恰好走 k 步的路径数（邻接矩阵的 k 次幂）
2.6 素数分布与素数定理
素数定理（Hadamard 和 de la Vallée-Poussin，1896年独立证明）：
π(n) ~ n / ln(n)
其中 π(n) 是不超过 n 的素数个数。更精确地：
lim_{n→∞} π(n) / (n / ln n) = 1
实际应用：这告诉我们在 [1, n] 中大约有 n/ln(n) 个素数。例如：
- n = 10⁶: 约 72,382 个素数（精确 78,498）
- n = 10⁹: 约 48,254,942 个素数（精确 50,847,534）
素数定理的实际意义：在 n 附近，平均每 ln(n) 个连续整数中有一个素数。这意味着：
- 如果你需要找一个 1024 位的素数（RSA），平均只需测试约 1024·ln(2) ≈ 710 个随机奇数
Level 3 · 规范怎么定义的
3.1 欧几里得算法的历史（公元前 300 年）
欧几里得算法记载于欧几里得的《几何原本》（Elements）第七卷（命题一和命题二），约公元前 300 年。这使它成为人类历史上最古老的非平凡算法之一——比中国的《九章算术》（约公元 200 年）中的更相减损术还早约 500 年。
原始表述：欧几里得使用的不是取模运算，而是反复相减。给定两条线段 AB 和 CD（AB > CD），反复从 AB 中减去 CD 的长度，直到剩余部分小于 CD，然后交换角色继续。这在数学上等价于取模运算（a mod b = a - ⌊a/b⌋·b），但取模版本收敛更快。
更相减损术：《九章算术》中记载的方法与欧几里得的减法版本本质相同："可半者半之，不可半者副置分母子之数以少减多更相减损求其等也"。它先尝试同时除以 2，然后用较小数反复减较大数。
Lamé 定理（1844）：Gabriel Lamé 证明了欧几里得算法的步数不超过 5 倍的输入位数。具体地，如果 a > b ≥ 1，则欧几里得算法的除法次数不超过 5·log₁₀(b) + 1。这是算法分析历史上最早的复杂度证明之一。
Lamé 的证明核心：他证明了使欧几里得算法需要 n 步的最小输入对恰好是相邻 Fibonacci 数 (F_{n+2}, F_{n+1})。因为 F_n ≈ φ^n/√5（φ = (1+√5)/2 ≈ 1.618），所以步数 n ≈ log_φ(b) = log₁₀(b) / log₁₀(φ) ≈ 2.078·log₁₀(b)，给出上界 5·log₁₀(b)。
3.2 费马小定理
定理（Pierre de Fermat，1640年书信中提出，Euler 1736年首次证明）：
若 p 是素数，a 是不被 p 整除的整数，则 a^(p-1) ≡ 1 (mod p)。
等价形式：对任意整数 a，a^p ≡ a (mod p)。
证明（用群论语言）：考虑模 p 的非零剩余类 {1, 2, ..., p-1}，它在乘法下构成一个 p-1 阶群（记为 (Z/pZ)*）。由拉格朗日定理，群中每个元素的阶整除群的阶，即 ord(a) | (p-1)。因此 a^(p-1) = (a^ord(a))^((p-1)/ord(a)) = 1^((p-1)/ord(a)) = 1。
组合证明（不用群论）：考虑 a, 2a, 3a, ..., (p-1)a 这 p-1 个数模 p 的值。它们两两不同余（否则 ia ≡ ja (mod p) 意味着 (i-j)a ≡ 0 (mod p)，但 p∤a 且 0 < |i-j| < p，矛盾）。所以它们恰好是 1, 2, ..., p-1 的某个排列。因此：
a · 2a · 3a · ... · (p-1)a ≡ 1 · 2 · 3 · ... · (p-1) (mod p)
a^(p-1) · (p-1)! ≡ (p-1)! (mod p)
由于 (p-1)! 与 p 互质（(p-1)! 的所有因子都小于 p），两边约去得 a^(p-1) ≡ 1 (mod p)。
应用：
- 求模逆元：a⁻¹ ≡ a^(p-2) (mod p)
- 素性测试（Fermat 测试）：若 a^(n-1) ≢ 1 (mod n)，则 n 一定不是素数
- 注意：存在 Carmichael 数（如 561 = 3·11·17），它们是合数但对所有 gcd(a, n)=1 的 a 都满足 a^(n-1) ≡ 1 (mod n)
3.3 中国剩余定理（CRT）
定理（《孙子算经》，约公元 3-5 世纪）：
设 m₁, m₂, ..., mₖ 两两互质，则同余方程组
x ≡ a₁ (mod m₁) x ≡ a₂ (mod m₂) ... x ≡ aₖ (mod mₖ)
在模 M = m₁·m₂·...·mₖ 下有唯一解。
历史背景：《孙子算经》中的原始问题是："今有物不知其数，三三数之余二，五五数之余三，七七数之余二，问物几何？"（答：23）。这是 x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7) 的联立。
构造性证明：
设 M = m₁·m₂·...·mₖ，Mᵢ = M/mᵢ。
因为 gcd(Mᵢ, mᵢ) = 1（Mᵢ 由其他所有 mⱼ 的乘积组成，与 mᵢ 互质），所以存在 Mᵢ 的模 mᵢ 逆元 tᵢ，使得 Mᵢ · tᵢ ≡ 1 (mod mᵢ)。
解为：x = Σᵢ aᵢ · Mᵢ · tᵢ (mod M)
验证：对 x 模 mⱼ，除了 j = i 的项（贡献 aⱼ · 1 = aⱼ），其他所有项的 Mᵢ 都含因子 mⱼ，模 mⱼ 为 0。
def chinese_remainder_theorem(remainders: list, moduli: list) -> int:
"""中国剩余定理
输入: remainders = [a1, a2, ..., ak]
moduli = [m1, m2, ..., mk]（两两互质）
输出: x 使得 x ≡ ai (mod mi) 对所有 i 成立
时间复杂度: O(k · log(max(mi)))
"""
M = 1
for m in moduli:
M *= m
x = 0
for ai, mi in zip(remainders, moduli):
Mi = M // mi
# 求 Mi 的模 mi 逆元
_, ti, _ = extended_gcd(Mi, mi)
x += ai * Mi * ti
return x % M
实际例子：
# 孙子算经问题
# x ≡ 2 (mod 3), x ≡ 3 (mod 5), x ≡ 2 (mod 7)
result = chinese_remainder_theorem([2, 3, 2], [3, 5, 7])
print(result) # 23
# 验证: 23 % 3 = 2 ✓, 23 % 5 = 3 ✓, 23 % 7 = 2 ✓
CRT 在实际系统中的应用：
-
RSA 加速（CRT-RSA）：RSA 解密时，不直接计算 c^d mod n（n 是 2048 位），而是分别计算 c^d mod p 和 c^d mod q（各 1024 位），再用 CRT 合并。这快约 4 倍（两个 1024 位的幂运算比一个 2048 位的快很多）。
-
大整数乘法：选择多个小素数 p₁, p₂, ...，分别在模 pᵢ 下做乘法，最后用 CRT 恢复结果。
-
竞赛中的多模数方法：当答案可能很大时，分别对 3 个大素数取模，再用 CRT 恢复精确值。
3.4 欧拉定理与欧拉函数
欧拉函数 φ(n)：小于 n 且与 n 互质的正整数个数。
def euler_phi(n: int) -> int:
"""欧拉函数
φ(n) = n · ∏(1 - 1/p)，其中 p 是 n 的所有不同质因子
时间复杂度: O(√n)
"""
result = n
p = 2
while p * p <= n:
if n % p == 0:
# 去除所有因子 p
while n % p == 0:
n //= p
result -= result // p # result *= (1 - 1/p)
p += 1
if n > 1: # n 本身是质因子
result -= result // n
return result
欧拉定理（Euler，1763）：若 gcd(a, n) = 1，则 a^φ(n) ≡ 1 (mod n)。
这是费马小定理的推广：当 n = p 是素数时，φ(p) = p-1，退化为费马小定理。
欧拉函数的性质：
- φ(1) = 1
- φ(p) = p - 1（p 是素数）
- φ(p^k) = p^k - p^(k-1) = p^(k-1)·(p-1)
- φ(mn) = φ(m)·φ(n) 当 gcd(m,n) = 1（积性函数）
- Σ_{d|n} φ(d) = n（n 的所有因子的欧拉函数之和等于 n）
3.5 Miller-Rabin 素性测试
对于大数（如 10^18），试除法太慢。Miller-Rabin 测试是一种概率素性测试，基于以下观察：
前提：若 p 是奇素数，写 p-1 = 2^s · d（d 是奇数）。对任意 1 < a < p：
- 要么 a^d ≡ 1 (mod p)
- 要么存在 0 ≤ r < s 使得 a^(2^r · d) ≡ -1 (mod p)
（这来源于费马小定理 a^(p-1) ≡ 1 和"平方根只有 ±1"的性质）
如果一个数 n 对某个 a 不满足上述条件，则 n 一定是合数。
def miller_rabin(n: int, witnesses: list = None) -> bool:
"""Miller-Rabin 素性测试
对于 n < 3,317,044,064,679,887,385,961,981，
使用 witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37] 即可确定性判断。
对于 n < 3.3×10^24，可以确定性判断（无概率误差）。
"""
if n < 2:
return False
if n == 2 or n == 3:
return True
if n % 2 == 0:
return False
# 写 n-1 = 2^s * d
s, d = 0, n - 1
while d % 2 == 0:
s += 1
d //= 2
if witnesses is None:
# 对 n < 3.3×10^24 确定性的 witness 集合
witnesses = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37]
for a in witnesses:
if a >= n:
continue
x = power_mod(a, d, n)
if x == 1 or x == n - 1:
continue
for _ in range(s - 1):
x = x * x % n
if x == n - 1:
break
else:
return False # n 是合数
return True # n （极大概率）是素数
确定性结果：Gary Miller（1976）证明了在广义黎曼假设（GRH）下，只需测试 a = 2, 3, ..., 2(ln n)² 即可确定性判断。实际中，对于 64 位整数（n < 2^64），固定的 12 个 witness 就够了。
Level 4 · 边界与陷阱
4.1 面试题：计数质数 (LeetCode #204)
题目：给定整数 n，返回小于 n 的所有素数的数量。
分析：这道题直接用埃氏筛。注意是"小于 n"而不是"不超过 n"。
def count_primes(n: int) -> int:
"""LeetCode 204: 计数质数
时间复杂度: O(n log log n)
空间复杂度: O(n)
"""
if n <= 2:
return 0
is_prime = [True] * n
is_prime[0] = is_prime[1] = False
i = 2
while i * i < n:
if is_prime[i]:
# 标记 i 的所有倍数
for j in range(i * i, n, i):
is_prime[j] = False
i += 1
return sum(is_prime)
面试跟进问题：
- "能否只用 O(√n) 空间？" —— 分段筛（Segmented Sieve）
- "如果 n = 10^10 呢？" —— Meissel-Lehmer 算法，O(n^(2/3)) 时间
- "如何判断单个大数是否为素数？" —— Miller-Rabin
优化技巧：可以用 bitset 代替 bool 数组节省 8 倍内存。Python 中可以使用 bytearray：
def count_primes_optimized(n: int) -> int:
"""优化版本：只处理奇数"""
if n <= 2:
return 0
if n <= 3:
return 1
# 只存奇数: sieve[i] 代表 2i+3 是否是素数
size = (n - 3) // 2 + 1
sieve = bytearray(b'\x01') * size # True
i = 0
while (2 * i + 3) ** 2 < n:
if sieve[i]:
p = 2 * i + 3
# p² 对应的索引
start = (p * p - 3) // 2
sieve[start::p] = bytearray(len(sieve[start::p]))
i += 1
return sum(sieve) + 1 # +1 for the prime 2
4.2 面试题：快乐数 (LeetCode #202)
题目：写一个算法来判断一个数 n 是不是快乐数。"快乐数"定义：对于一个正整数，每一次将该数替换为它每个位上数字的平方和，然后重复这个过程直到这个数变为 1，也可能是无限循环但始终变不到 1。如果可以变为 1，那么这个数就是快乐数。
数论分析：关键洞察是过程一定会进入循环（因为对于 d 位数，下一步的值不超过 81d，所以值的范围是有限的）。而且不快乐的数最终一定会进入包含 4 的循环：4 → 16 → 37 → 58 → 89 → 145 → 42 → 20 → 4。
def is_happy(n: int) -> bool:
"""LeetCode 202: 快乐数
方法一: Floyd 环检测
时间复杂度: O(log n)（实际上是 O(243 · 3) = O(1)）
空间复杂度: O(1)
"""
def digit_square_sum(num):
total = 0
while num:
num, digit = divmod(num, 10)
total += digit * digit
return total
slow = n
fast = digit_square_sum(n)
while fast != 1 and slow != fast:
slow = digit_square_sum(slow)
fast = digit_square_sum(digit_square_sum(fast))
return fast == 1
def is_happy_v2(n: int) -> bool:
"""方法二: 直接检查是否遇到 4"""
def digit_square_sum(num):
total = 0
while num:
num, digit = divmod(num, 10)
total += digit * digit
return total
while n != 1 and n != 4:
n = digit_square_sum(n)
return n == 1
为什么收敛很快？对于 n ≤ 999（3位数），下一步最大是 9²+9²+9² = 243。对于 n ≤ 9999999（7位数），下一步最大是 7×81 = 567。所以无论初始值多大，几步之内就会降到 1000 以内，然后在有限集合中循环。
4.3 面试题：Pow(x, n) (LeetCode #50)
题目：实现 pow(x, n)，即计算 x 的 n 次幂函数。
关键陷阱：
- n 可以是负数 → x^(-n) = 1 / x^n
- n 可以是 -2^31（取绝对值会溢出 int32）
- x = 0 且 n < 0 → 无定义（但题目保证不会这样）
def my_pow(x: float, n: int) -> float:
"""LeetCode 50: Pow(x, n)
快速幂，处理负指数和边界情况
时间复杂度: O(log |n|)
空间复杂度: O(1)
"""
if n == 0:
return 1.0
# 处理负指数
if n < 0:
x = 1.0 / x
n = -n
result = 1.0
current = x
while n > 0:
if n & 1:
result *= current
current *= current
n >>= 1
return result
面试中的讨论点：
- 为什么不能用
n = abs(n)？因为 abs(-2^31) 在 32 位整数中溢出（int32 范围是 [-2^31, 2^31 - 1]） - Python 没有这个问题（大整数），但在 C++/Java 中需要先转 long
- 递归版本更易理解但有栈开销
4.4 面试题：阶乘后的零 (LeetCode #172)
题目：给定一个整数 n，返回 n! 结果尾数中零的数量。
数论分析：尾零来自因子 10 = 2 × 5。在 n! 中，2 的因子个数远多于 5 的因子个数（每隔一个偶数贡献至少一个 2，但每隔 5 个数才贡献一个 5）。所以尾零数量 = n! 中 5 的因子总数。
Legendre 公式：n! 中素数 p 的幂次为 Σ_{i=1}^{∞} ⌊n/p^i⌋
def trailing_zeroes(n: int) -> int:
"""LeetCode 172: 阶乘后的零
计算 n! 末尾零的个数 = n! 中因子 5 的个数
时间复杂度: O(log₅ n)
空间复杂度: O(1)
"""
count = 0
power_of_5 = 5
while power_of_5 <= n:
count += n // power_of_5
power_of_5 *= 5
return count
# 更简洁的写法
def trailing_zeroes_v2(n: int) -> int:
count = 0
while n >= 5:
n //= 5
count += n
return count
直觉解释：
- ⌊n/5⌋ 计算有多少个 5 的倍数（每个贡献一个 5）
- ⌊n/25⌋ 计算有多少个 25 的倍数（每个额外贡献一个 5）
- ⌊n/125⌋ 计算有多少个 125 的倍数（再额外贡献一个 5）
- ...
例如 n = 100：⌊100/5⌋ + ⌊100/25⌋ + ⌊100/125⌋ = 20 + 4 + 0 = 24。所以 100! 末尾有 24 个零。
进阶问题："n! 在 k 进制下末尾零的个数" —— 对 k 做质因数分解，对每个质因子用 Legendre 公式计算，取最小值。
4.5 竞赛中的常见数论陷阱
陷阱 1：取模运算中的负数
# Python 中 -7 % 3 = 2（总是非负），这是正确的数学定义
# 但 C/C++ 中 -7 % 3 = -1（保留符号），需要手动处理
# C++ 中安全的取模: ((a % m) + m) % m
陷阱 2：模运算不满足除法分配律
# 错误: (a / b) % m ≠ (a % m) / (b % m)
# 正确: (a / b) % m = (a % m) * (b 的模逆元) % m
# 例: (10 / 2) % 7 = 5
# 错误做法: (10 % 7) / (2 % 7) = 3 / 2 = 1.5 ??
# 正确做法: (10 % 7) * mod_inverse(2, 7) % 7 = 3 * 4 % 7 = 12 % 7 = 5
# (因为 2 * 4 = 8 ≡ 1 (mod 7)，所以 2⁻¹ ≡ 4 (mod 7))
陷阱 3：快速幂中 base 为 0 或 mod 为 1
def power_mod_safe(base, exp, mod):
if mod == 1:
return 0 # 任何数 mod 1 = 0
if base == 0:
return 0 if exp > 0 else 1 # 0^0 通常定义为 1
# ... 正常快速幂
陷阱 4：组合数中 k > n
# C(n, k) 当 k > n 时应返回 0，不要让阶乘算出负数
def safe_comb(n, k, mod):
if k < 0 or k > n:
return 0
# ... 正常计算
陷阱 5：线性筛中的数组大小
在竞赛中，n = 10^7 时素数大约有 620,000 个。primes 数组开 n/ln(n) ≈ n/16 的大小就够了，但为安全起见通常开 n/10 或直接用动态数组。
4.6 综合应用：模运算下的完整工具箱
下面是竞赛/面试中数论问题的完整模板：
class NumberTheoryToolbox:
"""数论工具箱 - 竞赛/面试万用模板"""
def __init__(self, mod: int = 10**9 + 7):
self.mod = mod
def power(self, base: int, exp: int) -> int:
"""快速幂"""
result = 1
base %= self.mod
while exp > 0:
if exp & 1:
result = result * base % self.mod
base = base * base % self.mod
exp >>= 1
return result
def inverse(self, a: int) -> int:
"""模逆元（mod 必须是素数）"""
return self.power(a, self.mod - 2)
def factorial_and_inverse(self, n: int) -> tuple:
"""预处理阶乘和阶乘逆元"""
fact = [1] * (n + 1)
for i in range(1, n + 1):
fact[i] = fact[i - 1] * i % self.mod
inv_fact = [1] * (n + 1)
inv_fact[n] = self.power(fact[n], self.mod - 2)
for i in range(n - 1, -1, -1):
inv_fact[i] = inv_fact[i + 1] * (i + 1) % self.mod
return fact, inv_fact
def comb(self, n: int, k: int, fact: list, inv_fact: list) -> int:
"""O(1) 组合数查询"""
if k < 0 or k > n:
return 0
return fact[n] * inv_fact[k] % self.mod * inv_fact[n - k] % self.mod
@staticmethod
def gcd(a: int, b: int) -> int:
while b:
a, b = b, a % b
return a
@staticmethod
def lcm(a: int, b: int) -> int:
return a // NumberTheoryToolbox.gcd(a, b) * b
@staticmethod
def extended_gcd(a: int, b: int) -> tuple:
if b == 0:
return a, 1, 0
g, x, y = NumberTheoryToolbox.extended_gcd(b, a % b)
return g, y, x - (a // b) * y
@staticmethod
def crt(remainders: list, moduli: list) -> int:
"""中国剩余定理"""
M = 1
for m in moduli:
M *= m
result = 0
for a, m in zip(remainders, moduli):
Mi = M // m
_, t, _ = NumberTheoryToolbox.extended_gcd(Mi, m)
result += a * Mi * t
return result % M
4.7 面试应答策略
当面试中遇到数论问题时：
-
先确认边界：n 的范围决定算法选择
- n ≤ 10⁶: 埃氏筛 / 暴力
- n ≤ 10⁷: 线性筛
- n ≤ 10^18: 快速幂 / Miller-Rabin / Pollard Rho
-
模运算三原则：
- 加减乘可以随时取模
- 除法要转成乘逆元
- 每步取模防溢出
-
常见套路识别：
- "结果对 10^9+7 取模" → 需要模逆元、快速幂
- "求组合数" → 预处理阶乘 + Lucas
- "最大公约数" → 欧几里得
- "同余方程" → 扩展欧几里得 / CRT
-
时间复杂度速记：
- gcd: O(log n)
- 快速幂: O(log n)
- 埃氏筛 [1,n]: O(n log log n)
- 线性筛 [1,n]: O(n)
- 试除判素: O(√n)
- Miller-Rabin: O(k · log²n)，k 是 witness 数量