第 83 章 整除分块与数论进阶¶
配套例题：BISHI41 【模板】整除分块、BISHI42 余数求和、 BISHI39 【模板】Pollard-Rho 算法、BISHI40 数组取精 来源：S3 day10《数论与组合数学》（周尚彦）第 17–20、28–33 页； day9 resource《数论》（高铭鸿）第 4–6 页、《数论》（邹雨恒）第 21–26、80–106、133–138 页、 《数论与信息学竞赛》（胡渊鸣）第 19–21、38 页、《数论》（李有为）第 16–25 页 前置：80-数论基础、82-欧拉函数与欧拉降幂
这一章把三件「进阶但常考」的事讲清楚：
- 整除分块——把 \(\sum_{i=1}^n f(\lfloor n/i\rfloor)\) 从 \(O(n)\) 降到 \(O(\sqrt n)\)；
- Miller-Rabin / Pollard-Rho——大数的素性判定与分解；
- 莫比乌斯反演——把 \([\gcd(a,b)=1]\) 这种条件拆成可交换求和的形式。
Python 在这一章有一处巨大优势：Miller-Rabin 和 Pollard-Rho 在 C++ 里的一大半代码
都是在跟 long long 溢出搏斗（龟速乘、__int128、Montgomery 模乘），
而 Python 的 int 任意精度，这些全都不需要写。
83.1 整除分块（数论分块）¶
核心事实¶
高铭鸿的课件把它列成「有关整除的一些小性质」的第一条：
\(\lfloor n/i \rfloor\) 的不同结果数是 \(O(\sqrt n)\) 的。
证明：分两段看。 - \(i \le \sqrt n\)：这样的 \(i\) 至多 \(\sqrt n\) 个，所以贡献至多 \(\sqrt n\) 个不同值； - \(i > \sqrt n\)：此时 \(\lfloor n/i \rfloor < \sqrt n\)，取值只能落在 \([0, \sqrt n)\)， 同样至多 \(\sqrt n\) 个。
合计至多 \(2\sqrt n\) 个不同值。\(\square\)
周尚彦的课件给的提示是同一件事的另一个说法：
\(k \times \lfloor n/k \rfloor \le n\)，其中必有 \(k \le \sqrt n\) 或 \(\lfloor n/k \rfloor \le \sqrt n\)。
块的右端点¶
引理：设 \(v = \lfloor n/l \rfloor\)（\(v \ge 1\)），则使 \(\lfloor n/r \rfloor = v\) 的最大 \(r\) 是 $\(r = \left\lfloor \frac{n}{v} \right\rfloor\)$
证明：\(\lfloor n/r \rfloor = v \iff v \le n/r < v+1 \iff \frac{n}{v+1} < r \le \frac{n}{v}\)。 因为 \(r\) 是整数，最大的 \(r\) 就是 \(\lfloor n/v \rfloor\)。\(\square\)
于是可以「整块整块地跳」：l = 1 开始，每次算 v = n // l、r = n // v，
处理 \([l, r]\) 这一整段，然后 l = r + 1。
模板一：整除分块¶
# [片段] 模板：整除分块，枚举所有 (l, r, v) 三元组
def blocks(n, limit=None):
"""生成 (l, r, v)：区间 [l, r] 内 n//i 恒等于 v。共 O(sqrt n) 段。"""
# 上界必须收到 min(n, limit)：一旦 l > n 就有 v = n//l == 0，下一行会除零
m = n if limit is None else min(n, limit)
l = 1
while l <= m:
v = n // l # 本段的公共商，由左端点决定
# floor(n/r) == v 等价于 n/(v+1) < r <= n/v，所以最大的 r 就是 floor(n/v)。
# 这一行是整除分块的全部：把「逐个 i」变成「整段 [l, r] 一起处理」
r = min(m, n // v) # ★ 右端点要对 limit 取 min
yield l, r, v
l = r + 1 # 跳到下一段的左端点，商必然变小
def sum_floor(n):
"""Σ_{i=1..n} floor(n/i)，O(sqrt n)。"""
total = 0
l = 1
while l <= n:
v = n // l # l <= n 保证 v >= 1，n // v 不会除零
r = n // v # 同值区间 [l, r] 的右端点
total += v * (r - l + 1) # 这一段有 r-l+1 项，每项都等于 v
l = r + 1
return total
两个常见变体：
| 求和式 | 每块的贡献 |
|---|---|
| \(\sum \lfloor n/i \rfloor\) | \(v \times (r-l+1)\) |
| \(\sum i \cdot \lfloor n/i \rfloor\) | \(v \times \frac{(l+r)(r-l+1)}{2}\) |
| \(\sum f(i)\lfloor n/i\rfloor\) | \(v \times (F(r)-F(l-1))\)，\(F\) 是 \(f\) 的前缀和 |
| 二维 \(\sum \lfloor n/i\rfloor\lfloor m/i\rfloor\) | 右端点取 \(\min(n//(n//l),\ m//(m//l))\) |
v = n // l时v一定 \(\ge 1\)（因为 \(l \le n\)），所以n // v不会除零。 但若循环上界是 \(m > n\)（比如 BISHI42 里 \(n > k\) 的情形）， 就必须先把上界收到 \(\min(n, m)\)，否则n // l == 0会ZeroDivisionError。
一个更快的等价写法：数格点¶
对 \(\sum_{i=1}^{n}\lfloor n/i\rfloor\) 这个最经典的式子，还有一个常数更小的写法。
注意到 \(\lfloor n/i \rfloor\) 就是「\(i\) 的倍数中不超过 \(n\) 的个数」，所以
即双曲线 \(ij=n\) 下方的格点数。以 \(s = \lfloor\sqrt n\rfloor\) 为界做容斥 （这就是数论里的「双曲线法 / Dirichlet 双曲线技巧」）：
证明：把格点分成「\(i \le s\)」和「\(j \le s\)」两块。 每块各有 \(\sum_{i=1}^{s}\lfloor n/i\rfloor\) 个格点（第二块按 \(j\) 数）。 两块的交集是 \(\{(i,j): i\le s, j\le s\}\) 里满足 \(ij\le n\) 的格点—— 而 \(i,j\le s\) 时 \(ij \le s^2 \le n\) 自动成立，所以交集恰是 \(s\times s\) 的完整正方形。 并集的元素恰好覆盖全部满足 \(ij\le n\) 的格点（若 \(i>s\) 且 \(j>s\) 则 \(ij>n\)）。 由容斥即得。\(\square\)
为什么在 Python 里更划算：循环体只有一次整除加一次加法，
可以直接写成 sum(n // i for i in range(1, s + 1)) 让求和落到 C 层；
而整除分块的循环体有 4–5 个操作且无法向量化。
83.2 约数个数与约数和¶
设 \(n = \prod p_i^{a_i}\)：
（推导见 80 章。陈许旻的课件把这两条列在「因子个数」一节。）
前缀和：用整除分块¶
证明：交换求和顺序。 \(\sum_{i\le n} d(i) = \sum_{i\le n}\sum_{d\mid i} 1 = \sum_{d\le n}\#\{i\le n : d\mid i\} = \sum_{d\le n}\lfloor n/d\rfloor\)。 约数和同理，只是每个约数 \(d\) 贡献 \(d\) 而不是 1。\(\square\)
这解释了 BISHI41 为什么值得单独出成模板题：\(\sum\lfloor n/i\rfloor\) 既是格点计数， 也是约数个数的前缀和，还是整除分块的最小实例。
求 \(1..n\) 每个数的约数个数 / 约数和¶
不用分解，直接倍数枚举（调和级数 \(O(n\log n)\)）：
# [片段] 模板：O(n log n) 求 1..n 的约数个数表与约数和表
def divisor_tables(n):
d = [0] * (n + 1) # 约数个数
s = [0] * (n + 1) # 约数和
# 不做分解，改成「枚举约数 i，去更新它的所有倍数 j」：
# 内层次数是 n/i，全部加起来是调和级数 n·ln n
for i in range(1, n + 1):
for j in range(i, n + 1, i): # i 是 j 的约数
d[j] += 1 # j 多了一个约数 i
s[j] += i # 约数和累加这个约数本身
return d, s
Python 现实性：内层总迭代次数是 \(n\ln n\)，\(n=10^6\) 时约 \(1.4\times10^7\) 次， 约 4 秒——危险。\(n \le 2\times10^5\) 才稳。 若只需要约数个数，可以改用线性筛（\(d\) 是积性函数），但常数同样不小。 能用整除分块只求前缀和的，就不要把整张表建出来。
枚举单个数的全部约数¶
# [片段] 模板：O(sqrt n) 枚举 n 的全部约数
import math
def divisors(n):
"""返回 n 的全部约数（无序）。O(sqrt n)。"""
res = []
# 约数成对出现：i 与 n//i 中必有一个不超过 sqrt(n)，所以只枚举小的那半
for i in range(1, math.isqrt(n) + 1):
if n % i == 0:
res.append(i) # 小的那个
if i != n // i: # ★ 完全平方数时别把 sqrt 算两次
res.append(n // i) # 配对的大的那个
return res
83.3 Miller-Rabin 素性测试¶
试除法是 \(O(\sqrt n)\)，\(n = 10^{18}\) 时不可行；询问次数多时连 \(n=10^{12}\) 都撑不住。
从费马测试到 Miller-Rabin¶
费马测试：若 \(n\) 是质数且 \(\gcd(a,n)=1\)，则 \(a^{n-1}\equiv1\pmod n\)。 反过来用：随机取 \(a\)，若 \(a^{n-1}\not\equiv1\)，则 \(n\) 一定是合数。
胡渊鸣的课件指出了它的问题：
反例：\(2^{340}\equiv 1 \pmod{341}\)。反例确实不多, 但是不能忽略。
（\(341 = 11\times31\)。更糟的是 Carmichael 数，如 \(561\)，对所有与之互质的底数都通过费马测试。）
加强的依据是二次探测定理：
二次探测定理：\(p\) 是奇质数，若 \(x^2 \equiv 1 \pmod p\)，则 \(x \equiv \pm1 \pmod p\)。
证明：\(p \mid (x-1)(x+1)\)，\(p\) 是质数所以 \(p\mid(x-1)\) 或 \(p\mid(x+1)\)。\(\square\)
于是把 \(n-1\) 写成 \(d\cdot2^s\)（\(d\) 为奇数），从 \(a^d\) 开始不断平方， 最终要到达 \(a^{n-1}\equiv1\)。如果 \(n\) 是质数，这条平方链第一次变成 1 之前的那一项必须是 \(-1\)。
Miller-Rabin 判定：对底数 \(a\)，\(n\) 通过测试当且仅当 $\(a^{d}\equiv1 \pmod n \quad\text{或}\quad \exists\, 0\le r<s,\ a^{d\cdot2^r}\equiv-1\pmod n\)$ 否则 \(n\) 必为合数。
单个底数的误判概率 \(\le 1/4\)，\(k\) 个随机底数误判概率 \(\le 4^{-k}\)。
确定性底数：竞赛里不要用随机¶
已知的确定性结论（可以直接背）：
| \(n\) 的上界 | 底数集合 |
|---|---|
| \(3.2\times10^9\)（\(2^{32}\) 级） | \(\{2, 3, 5, 7\}\) |
| \(3.47\times10^{12}\) | \(\{2, 3, 5, 7, 11, 13\}\) |
| \(3.4\times10^{14}\) | \(\{2,3,5,7,11,13,17\}\) |
| \(3.3\times10^{24}\)（含全部 \(2^{64}\)） | 前 13 个质数 \(\{2,3,\ldots,41\}\) |
| \(< 2^{64}\)（更快） | \(\{2, 325, 9375, 28178, 450775, 9780504, 1795265022\}\)（7 个） |
（胡渊鸣课件提到的 \(\{2,3,7,61,24251\}\) 在 \(10^{16}\) 内只有唯一强伪素数 \(46856248255981\)， 所以那组底数不能直接用，要额外特判——用上表更省心。）
为什么竞赛里坚持确定性底数？ 随机底数意味着同一份代码同一组数据可能这次 AC 下次 WA， 调试时无法复现。确定性底数让程序完全可预测。
模板二：Miller-Rabin¶
# [片段] 模板：Miller-Rabin，确定性底数，适用于 n < 3.3e24
_SMALL = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41)
def is_prime(n):
"""确定性 Miller-Rabin。Python 的 pow(a,d,n) 是 C 层模幂，不需要龟速乘。"""
if n < 2:
return False
# 先用小质数试除：既筛掉绝大多数合数，也保证后面 a % n != 0
for p in _SMALL:
if n % p == 0:
return n == p # 小质数本身算质数，其余倍数直接判负
# 把 n-1 的因子 2 全部提出来，得到 n - 1 = d * 2^s（d 为奇数）。
# 二次探测就沿着 a^d, a^(2d), a^(4d), ... 这条平方链检查
d = n - 1
s = 0
while d % 2 == 0: # n - 1 = d * 2^s，d 为奇数
d //= 2
s += 1
for a in _SMALL: # n < 3.3e24 时这 13 个底数是确定性的
x = pow(a, d, n) # ★ 内置快速幂，C 实现
if x == 1 or x == n - 1:
continue # 链首就是 ±1，这个底数不能证伪，换下一个
# 还剩 s-1 次平方机会：最后一次平方得到的是 a^(n-1)，不必再看
for _ in range(s - 1):
x = x * x % n # ★ Python 不会溢出，C++ 这里要 __int128
if x == n - 1:
break # 出现 -1，符合质数应有的形态
else:
# for 的 else 只在没 break 时执行：整条链都没出现 -1，
# 却最终会平方到 1，违反二次探测定理
return False # 平方链里没出现 -1 -> 合数
return True
for ... else的用法：else分支只在循环没有break时执行。 这里正好表达「整条平方链都没出现 \(-1\)」。 见 10-条件与循环。
Python 相对 C++ 的三处便宜：
| C++ 要做的事 | Python |
|---|---|
\(n>10^{18}\) 时 a*b%n 溢出 → 龟速乘 / __int128 / Montgomery |
a * b % n 直接写 |
| 手写快速幂 | pow(a, d, n) |
| 随机数生成器、种子管理 | 确定性底数，一个随机数都不用 |
83.4 Pollard-Rho 分解¶
试除法分解是 \(O(\sqrt n)\)；Pollard-Rho 期望 \(O(n^{1/4})\)。 李有为的课件一句带过：「用 Pollard rho 可以做到 \(O(n^{1/4})\) 分解」。
原理（生日悖论）：取伪随机函数 \(f(x) = (x^2+c)\bmod n\)，从 \(x_0\) 迭代得到序列。 设 \(p\) 是 \(n\) 的一个质因子，那么这个序列模 \(p\) 的周期期望是 \(O(\sqrt p) = O(n^{1/4})\)。 一旦出现 \(x_i \equiv x_j \pmod p\) 但 \(x_i \ne x_j\)， \(\gcd(|x_i-x_j|, n)\) 就是 \(n\) 的一个非平凡因子。
用 Floyd 判圈（龟兔）或 Brent 的批量取 gcd 来找这样的 \((i,j)\)。
模板三：Pollard-Rho（Brent 变体）¶
# [片段] 模板：Pollard-Rho 完整分解，配合上面的 is_prime 使用
import math
import random
def _pollard(n):
"""返回 n 的一个非平凡因子。要求 n 是合数且 n 不是 2 的倍数。"""
if n % 2 == 0:
return 2
# 外层 while：一组 (c, 起点) 可能只找出平凡因子，失败就换一组重来
while True:
c = random.randrange(1, n)
f = lambda x: (x * x + c) % n # 伪随机迭代函数，c 是这一轮的随机偏移
x = y = random.randrange(0, n)
d = 1
# Brent：把多步的差乘在一起，批量取一次 gcd，省掉大量 gcd 调用
while d == 1:
prod = 1
for _ in range(128): # 128 步为一批，批越大 gcd 调用越少
x = f(x) # 龟：走一步
y = f(f(y)) # 兔：走两步
if x == y:
break # 龟兔相遇，本轮迭代已进入环，无解
prod = prod * abs(x - y) % n
if prod == 0: # 乘出 0，退回逐步取 gcd
break
# 一次 gcd 覆盖整批：只要批内任一个 |x-y| 与 n 有公因子，就会被抓到
d = math.gcd(prod, n)
if x == y or prod == 0:
break
if 1 < d < n: # 非平凡因子才算成功，否则换 c 重来
return d
def factorize_big(n):
"""大整数完整分解，返回 {p: e}。期望 O(n^(1/4))。"""
res = {}
stack = [n] # 用显式栈代替递归，不受递归深度限制
while stack:
v = stack.pop()
if v == 1:
continue # 分解到 1 就没有质因子了
if is_prime(v): # 递归边界靠 Miller-Rabin 判定
res[v] = res.get(v, 0) + 1 # 同一质因子可能被拆出多次，累加重数
continue
d = _pollard(v) # 拆出一个非平凡因子，两半继续分解
stack.append(d)
stack.append(v // d)
return res
什么时候真的需要 Pollard-Rho？ \(n \le 10^{12}\) 时试除只要 \(10^6\) 次循环，Python 下 0.3 秒， 单组数据用试除就够。只有「\(T\) 组 × \(n \ge 10^{14}\) 要分解」才必须上 Pollard-Rho。 牛客这道 BISHI39 虽然挂着 Pollard-Rho 的名字，实际只要 Miller-Rabin（见例题）。
83.5 莫比乌斯函数与反演¶
定义¶
邹雨恒的课件一句话点破本质：
本质上是一种容斥原理所产生的系数。
\(\mu\) 是积性函数，\(\mu(p) = -1\)，\(\mu(p^a) = 0\ (a\ge2)\)。
最关键的一条式子¶
证明：\(n=1\) 时显然。\(n>1\) 时设 \(n\) 的不同质因子有 \(k\) 个（\(k\ge1\)）。 含平方因子的 \(d\) 贡献 0，所以只需考虑「不同质因子的子集」形成的 \(d\)， 选出 \(j\) 个质因子的 \(d\) 有 \(\binom kj\) 个，各贡献 \((-1)^j\)。于是 $\(\sum_{d\mid n}\mu(d) = \sum_{j=0}^{k}\binom kj(-1)^j = (1-1)^k = 0 \qquad \square\)$
这条式子的用法（邹雨恒课件的「小总结」）：
注意 \(\sum_{d\mid n}\mu(d)=[n=1]\)，就能把 \([\gcd(a,b)=1]\) 转化为 \([d\mid\gcd(a,b)]\)， 整除某两个数最大公约数就好做多了。
具体地：
最后一步再套整除分块就是 \(O(\sqrt N)\)。 「莫比乌斯 + 整除分块」是这类题的标准组合拳。
莫比乌斯反演¶
定理：若 \(F(n) = \sum_{d\mid n} f(d)\)，则 \(f(n) = \sum_{d\mid n}\mu\!\left(\frac nd\right)F(d)\)。
证明： $\(\sum_{d\mid n}\mu\left(\frac nd\right)F(d) = \sum_{d\mid n}\mu\left(\frac nd\right)\sum_{e\mid d} f(e) = \sum_{e\mid n} f(e)\sum_{d:\ e\mid d\mid n}\mu\left(\frac nd\right)\)$ 内层令 \(d = e t\)，则 \(t \mid \frac ne\)，内层和为 \(\sum_{t\mid n/e}\mu\left(\frac{n/e}{t}\right) = [n/e=1]\)。 于是只剩 \(e=n\) 一项，即 \(f(n)\)。\(\square\)
与欧拉函数的对偶¶
同一套手法也可以用 \(\varphi\)（见 82 章 的 \(\sum_{d\mid n}\varphi(d)=n\)）。 邹雨恒课件的经典推导：
记忆法：\([\gcd=1]\) 用 \(\mu\)，\(\gcd\) 本身用 \(\varphi\)。
模板四：线性筛莫比乌斯函数¶
# [片段] 模板：线性筛同时求 μ 与素数表，O(n)
def sieve_mu(n):
"""返回 (primes, mu)，mu[i] = μ(i)，i = 0..n。"""
mu = [0] * (n + 1)
mu[1] = 1 # μ(1) = 1，空的质因子集合
primes = []
is_comp = bytearray(n + 1) # 1 字节一格，比 list 省 8 倍内存
for i in range(2, n + 1):
if not is_comp[i]: # 没被任何更小的质数筛掉 -> i 是质数
primes.append(i)
mu[i] = -1 # μ(p) = -1
# 用 i 乘上「不超过 i 的最小质因子」的那些质数去筛，保证每个合数
# 只被它的最小质因子筛掉一次，总次数是 O(n)
for p in primes:
v = i * p
if v > n:
break # 越界，后面的 p 只会更大
is_comp[v] = 1 # v 的最小质因子就是 p
if i % p == 0:
# p 已经是 i 的质因子，v = i*p 含 p²，μ(v) = 0；
# 同时这个 break 保证再往后的质数不会重复筛到同一个合数
mu[v] = 0 # p² | v -> μ = 0
break
mu[v] = -mu[i] # 多一个不同质因子，符号取反
return primes, mu
83.6 狄利克雷卷积入门¶
定义：\((f * g)(n) = \sum_{d\mid n} f(d)\, g\!\left(\frac nd\right)\)。
几个基本函数：
| 记号 | 定义 |
|---|---|
| \(\varepsilon(n)\) | \([n=1]\)（卷积单位元） |
| \(\mathbf1(n)\) | 恒为 1 |
| \(\operatorname{id}(n)\) | \(n\) |
| \(\mu, \varphi, d, \sigma\) | 见前文 |
必须记住的几条卷积恒等式：
| 恒等式 | 等价的经典式子 |
|---|---|
| \(\mu * \mathbf1 = \varepsilon\) | \(\sum_{d\mid n}\mu(d)=[n=1]\) |
| \(\varphi * \mathbf1 = \operatorname{id}\) | \(\sum_{d\mid n}\varphi(d)=n\) |
| \(\mathbf1 * \mathbf1 = d\) | 约数个数 |
| \(\operatorname{id} * \mathbf1 = \sigma\) | 约数和 |
| \(\varphi = \mu * \operatorname{id}\) | 由前两条 + 结合律推出 |
性质（邹雨恒课件）：
- 狄利克雷卷积满足交换律、结合律，对加法有分配律；
- 两个积性函数的狄利克雷卷积仍是积性函数；
- 莫比乌斯反演就是「两边卷 \(\mu\)」：\(F = f * \mathbf1 \iff f = F * \mu\)。
为什么这套语言有用：把「\(\sum_{d\mid n}\) 型恒等式」变成代数运算之后， 推导可以像解方程一样做。比如要证 \(\varphi = \mu * \operatorname{id}\)： \(\varphi * \mathbf1 = \operatorname{id}\)，两边卷 \(\mu\)， 得 \(\varphi * (\mathbf1 * \mu) = \operatorname{id} * \mu\)，即 \(\varphi * \varepsilon = \mu * \operatorname{id}\)。\(\square\)
进阶：杜教筛能在 \(O(n^{2/3})\) 内求 \(\mu\)、\(\varphi\) 的前缀和， 核心是邹雨恒课件最后给的式子 \(G(n) = \sum_{i=1}^{n} F(\lfloor n/i\rfloor)\)—— 它本身又是整除分块。本教程不展开。
83.7 例题¶
BISHI41 【模板】整除分块（中等）¶
给定 \(1 \le n \le 10^{12}\)，求 \(\sum_{i=1}^{n}\lfloor n/i\rfloor\)。 题面见 BISHI41 原题（牛客）。
两条路都是 \(O(\sqrt n)\)：整除分块，或 83.1 末尾的数格点容斥。
本题解取后者，因为循环体更短、可以整个交给 C 层的 sum。
import sys
from math import isqrt
n = int(sys.stdin.buffer.read().split()[0])
s = isqrt(n) # 精确整数开方；n 到 1e12 时 int(n ** 0.5) 可能差 1
# 数格点 (i, j) 且 i*j <= n：两块各算一次，重叠的 s*s 减掉。
# 生成器交给内置 sum，1e6 次迭代全在 C 层，比整除分块的循环体更省
print(2 * sum(n // i for i in range(1, s + 1)) - s * s)
验算 \(n=10\)：\(s=3\)，\(\sum_{i=1}^{3}\lfloor10/i\rfloor = 10+5+3=18\)，\(2\times18-9=27\)； 手算 \(10+5+3+2+2+1+1+1+1+1=27\) ✓
三个坑：
math.isqrt而非int(n ** 0.5)。\(n\) 到 \(10^{12}\) 时浮点开方可能差 1， 容斥的正方形边长就错了，答案会偏差 \(O(\sqrt n)\)；- 容斥减掉的是 \(s^2\)（重叠的正方形），不是 \(s\)；
- \(n=1\) 时 \(s=1\)，\(2\times1-1=1\) ✓，边界自然成立。
规模：\(\sqrt{10^{12}} = 10^6\) 次 C 层迭代，约 0.1–0.2 秒。
答案量级 \(n\ln n \approx 2.8\times10^{13}\)，C++ 要 long long，Python 无忧。
题解见 solutions/BISHI41.py。
BISHI42 余数求和（中等）¶
给定 \(1 \le n, k \le 10^9\)，求 \(\sum_{i=1}^{n}(k \bmod i)\)。 题面见 BISHI42 原题（牛客）。
第一步永远是把取模拆开（周尚彦课件的原话：「按照取模定义式把原式改写」）：
当 \(i > k\) 时 \(\lfloor k/i\rfloor = 0\)，所以右边那个和只需算到 \(m=\min(n,k)\)。 再用整除分块，每块贡献 \(v\cdot\frac{(l+r)(r-l+1)}{2}\)。
import sys
def main():
n, k = map(int, sys.stdin.buffer.read().split()[:2])
m = min(n, k) # i > k 的项 floor(k/i) = 0，不必进循环
total = 0 # 累加 Σ i * floor(k/i)
l = 1
while l <= m:
v = k // l # l <= m <= k 保证 v >= 1，下一行不会除零
# 由 floor(k/r) = v 的最大 r 是 floor(k/v)；再对 m 取 min 防止越出求和上界
r = min(m, k // v) # [l, r] 内 floor(k/i) 恒等于 v
# 这一段的贡献是 v 乘上 l..r 的等差和；`// 2` 必须是整除，
# 用 `/ 2` 会转成浮点，5e17 量级直接丢精度
total += v * (l + r) * (r - l + 1) // 2
l = r + 1 # 跳到下一段
# k mod i = k - i*floor(k/i)，对 i = 1..n 求和即 n*k - total
print(n * k - total)
main()
验算样例 \(n=10, k=5\)：\(nk=50\)； \(\sum i\lfloor5/i\rfloor = 1\cdot5+2\cdot2+3\cdot1+4\cdot1+5\cdot1 = 5+4+3+4+5=21\)； \(50-21=29\) ✓（手算 \(0+1+2+1+0+5+5+5+5+5=29\)）。
四个坑：
- \(n\) 可能大于 \(k\)。此时 \(i\in(k,n]\) 的每项都等于 \(k\)， 这部分被 \(nk\) 和「\(\sum\) 只算到 \(\min(n,k)\)」自动吸收，不要重复加；
- 右端点必须对 \(m\) 取 min，否则 \(i\) 会越界到 \(n\) 以外（且 \(l>k\) 时
k//l == 0会除零）； - 等差数列求和用整数
// 2，不要用/ 2（Python 3 的/是浮点， \(5\times10^{17}\) 量级会丢精度，见 23-浮点与科学计数法）； - 答案量级最大约 \(nk/2 \approx 5\times10^{17}\)，C++ 必须
long long。
规模：\(O(\sqrt k)\approx6.3\times10^4\) 次循环，毫秒级。
题解见 solutions/BISHI42.py。
BISHI39 【模板】Pollard-Rho算法（中等）¶
\(T \le 10^5\) 组，每组给 \(1 \le x \le 10^{12}\)，判断是否为质数，输出
Yes/No。 题面见 BISHI39 原题（牛客）。
题目名字是个陷阱：它挂着 Pollard-Rho 的名字，但问的只是素性判定， 根本不需要把 \(x\) 分解出来。核心是 Miller-Rabin。 （Pollard-Rho 是「找非平凡因子」的算法，它内部还要靠 Miller-Rabin 判断递归边界。）
\(x \le 10^{12} < 3.47\times10^{12}\)，所以底数集合 \(\{2,3,5,7,11,13\}\) 是确定性的。
import sys
SMALL = (2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
53, 59, 61, 67, 71, 73, 79, 83, 89, 97)
# 对 n < 3.47e12 而言，这组底数的 Miller-Rabin 是确定性的
BASES = (2, 3, 5, 7, 11, 13)
def is_prime(n):
if n < 2:
return False # 1 和 0 都不是质数，且会让下面的逻辑失效
# 100 以内的小素数试除：绝大多数合数在这里就出局，主循环只留给大素数与大半素数。
# 顺带保证了 n > 97，底数 a 不会等于 n 本身（否则 a mod n = 0，判定失真）
for p in SMALL:
if n % p == 0:
return n == p # 小素数本身算质数，其余倍数直接判负
d = n - 1 # 把 n-1 写成 d * 2^s，d 为奇数
s = 0
while d % 2 == 0:
d //= 2
s += 1
for a in BASES:
x = pow(a, d, n) # 内置快速幂，C 实现
if x == 1 or x == n - 1:
continue # 链首已是 ±1，这个底数证伪不了
for _ in range(s - 1): # 还剩 s-1 次平方；最后一次得到的是 a^(n-1)
x = x * x % n # Python 大整数不溢出，C++ 这里要 __int128
if x == n - 1:
break # 出现 -1，符合质数应有的形态
else:
return False # 没有出现 -1，判定为合数
return True
def main():
data = sys.stdin.buffer.read().split()
t = int(data[0])
cache = {} # T 到 1e5 时重复询问很常见，记住算过的答案
out = []
for i in range(1, t + 1):
x = int(data[i])
r = cache.get(x)
if r is None: # 用 None 而不是假值判断，答案本身是字符串
r = "Yes" if is_prime(x) else "No"
cache[x] = r
out.append(r)
sys.stdout.write("\n".join(out) + "\n") # 一次性输出，10 万行逐行 print 会被 IO 拖死
main()
两处压常数的手段（\(T=10^5\) 时都必要）：
- 先用 100 以内的小素数试除。绝大多数合数在这里就被筛掉， 只有「大素数 / 大半素数」才真正走进 Miller-Rabin 的主循环；
- 用字典缓存已算过的 \(x\)。\(T\) 到 \(10^5\) 时重复询问很常见。
四个坑：
- \(x=1\) 不是质数，必须特判（\(n-1=0\) 会让 Miller-Rabin 的逻辑失效）；
- 底数 \(a\) 可能等于 \(n\) 本身（\(n=2,3,5,7,11,13\)），此时 \(a \bmod n = 0\)， 所以要先把小素数直接判 Yes 再进主流程；
- 输出是
Yes/No（首字母大写），不是YES/NO； - 不需要随机数——确定性底数让同一输入多次运行结果必然一致。
Python 版的代码比 C++ 短一半：C++ 在 \(x\sim10^{12}\) 时
x * x % n就已经越过long long（\(10^{24}\)），必须写__int128或龟速乘； Python 直接写x * x % n。这是本章开头说的那处「巨大优势」。
题解见 solutions/BISHI39.py。
BISHI40 数组取精（中等）¶
两个长度 \(n \le 10^5\) 的正整数序列 \(A, B\)。求一个下标子集 \(P\)， \(|P| \le \lfloor n/2\rfloor+1\)，且 \(2\sum_{i\in P}a_i > \sum a_i\)、\(2\sum_{i\in P}b_i > \sum b_i\) （两边都严格过半）。输出任意一组解。 题面见 BISHI40 原题（牛客）。
这题其实不是数论题，大纲把它排在本章，但真正的考点是配对贪心构造。 放在这里当作「数论章里的一道换脑子的题」也无妨——它的证明手法（构造单射） 在数论推导里同样常见。
构造：
- 把下标按 \(a\) 从大到小排序，得 \(o_1, o_2, \ldots, o_n\)；
- 必选 \(o_1\)（\(a\) 最大的那个）；
- 把 \(o_2..o_n\) 两两分组 \((o_2,o_3), (o_4,o_5), \ldots\)，每组选 \(b\) 更大的那个；
- 若 \(n\) 为偶数，\(o_n\) 落单，也一并选上。
选出的个数 \(=1+\lfloor(n-1)/2\rfloor\ (+1\text{ 若 }n\text{ 偶}) = \lfloor n/2\rfloor+1\)，正好卡满上限。
为什么 \(a\) 那边一定过半？ 把每个「未选中」的元素映射到「前一组里选中的元素」 （第 1 组的未选者映射到 \(o_1\)）。由于按 \(a\) 降序排，前一组的任意元素的 \(a\) 都不小于后一组的任意元素，所以这是一个单射且每对都满足「选中 \(\ge\) 未选中」。 选中集合比未选中集合多至少一个元素（最后一组的选中者没有被映射到）， 而 \(a_i \ge 1 > 0\)，所以 \(\sum_{\text{选中}} > \sum_{\text{未选中}}\)，即 $2\sum_{\text{选中}} > $ 总和。
\(b\) 那边更直接：每组里选的就是 \(b\) 大的那个，再加上白送的 \(o_1\)（\(b\ge1>0\)）， 同样严格过半。\(\square\)
import sys
def main():
data = sys.stdin.buffer.read().split()
n = int(data[0])
a = [int(x) for x in data[1:n + 1]] # 前 n 个是 A，紧接着 n 个是 B
b = [int(x) for x in data[n + 1:2 * n + 1]]
# 排序键取负号实现降序，排的是下标而不是值——最后要输出的是下标
order = sorted(range(n), key=lambda i: -a[i]) # 按 a 降序
pick = [order[0]] # a 最大的必选
i = 1
# 从第二名起两两一组：组内 a 的差距已被降序压住，于是可以放心按 b 取舍
while i + 1 < n: # 成对处理，组内取 b 大者
x, y = order[i], order[i + 1]
pick.append(x if b[x] >= b[y] else y)
i += 2
if i < n: # n 为偶数时最后一个落单
pick.append(order[i]) # 落单的也收下，仍不超过上限
# 输出下标是 1-based，所以统一 +1
sys.stdout.write("%d\n%s\n" % (len(pick),
" ".join(str(p + 1) for p in pick)))
main()
五个坑：
- 是严格过半（\(2\Sigma >\) 总和），「恰好一半」不算；
- \(n=1\) 时答案就是 \(\{1\}\)（\(\lfloor1/2\rfloor+1=1\)），\(2a_1>a_1\) 成立；
- \(n=2\) 时两个都要选，此时 \(2\Sigma = 2\times\)总和 \(>\) 总和；
- 排序键必须是 \(a\)，组内取舍必须看 \(b\)，反过来就证不出来；
- 输出下标是 1-based，且答案不唯一，本地必须 special judge
（本仓库的
solutions/_spj/BISHI40.py就是干这个的）。
题解见 solutions/BISHI40.py。
83.8 本章速查¶
公式¶
| 名称 | 内容 |
|---|---|
| 整除分块取值数 | \(\lfloor n/i\rfloor\) 只有 \(O(\sqrt n)\) 种取值 |
| 块右端点 | \(v=\lfloor n/l\rfloor\) 时最大的 \(r = \lfloor n/v\rfloor\) |
| 格点容斥 | \(\sum_{i\le n}\lfloor n/i\rfloor = 2\sum_{i\le s}\lfloor n/i\rfloor - s^2\)，\(s=\lfloor\sqrt n\rfloor\) |
| 取模拆分 | \(k\bmod i = k - i\lfloor k/i\rfloor\) |
| 约数个数前缀和 | \(\sum_{i\le n} d(i) = \sum_{i\le n}\lfloor n/i\rfloor\) |
| 约数和前缀和 | \(\sum_{i\le n}\sigma(i) = \sum_{i\le n} i\lfloor n/i\rfloor\) |
| 二次探测 | \(p\) 质数、\(x^2\equiv1 \Rightarrow x\equiv\pm1 \pmod p\) |
| Miller-Rabin | \(n-1=d2^s\)，检查 \(a^d\equiv1\) 或 \(\exists r,\ a^{d2^r}\equiv-1\) |
| 莫比乌斯 | \(\sum_{d\mid n}\mu(d)=[n=1]\) |
| 莫比乌斯反演 | \(F=f*\mathbf1 \iff f=F*\mu\) |
| \(\gcd=1\) 计数 | \(\sum_{d}\mu(d)\lfloor N/d\rfloor\lfloor M/d\rfloor\) |
| \(\gcd\) 求和 | \(\sum_{d}\varphi(d)\lfloor N/d\rfloor\lfloor M/d\rfloor\) |
| 卷积恒等式 | \(\mu*\mathbf1=\varepsilon\)，\(\varphi*\mathbf1=\operatorname{id}\)，\(\mathbf1*\mathbf1=d\)，\(\varphi=\mu*\operatorname{id}\) |
确定性 Miller-Rabin 底数¶
| \(n <\) | 底数 |
|---|---|
| \(3.2\times10^9\) | 2, 3, 5, 7 |
| \(3.47\times10^{12}\) | 2, 3, 5, 7, 11, 13 |
| \(3.3\times10^{24}\) | 前 13 个质数 |
Python 取舍¶
| 场景 | 做法 |
|---|---|
| \(\sum\lfloor n/i\rfloor\) | 格点容斥 + sum(生成器)，循环下沉到 C |
| 一般整除分块 | while l <= m: v = n//l; r = min(m, n//v); ... |
| 上界 \(m>n\) 时 | 先 m = min(n, m)，否则 n//l == 0 除零 |
| 平方根 | math.isqrt，永远 |
| 等差求和 | (l+r)*(r-l+1)//2，整数除法 |
| 大数模乘 | 直接 a*b%n，不需要龟速乘 / __int128 |
| 素性判定 | pow(a,d,n) + 确定性底数，不用随机 |
| \(n\le10^{12}\) 单次分解 | 试除法就够，别急着上 Pollard-Rho |
| 线性筛 \(\mu\)/\(\varphi\) | \(n \le 2\times10^6\)（纯 Python 双重循环） |
看到什么 → 想到什么¶
| 题面特征 | 第一反应 |
|---|---|
| 式子里有 \(\lfloor n/i\rfloor\) 或 \(n\bmod i\) | 整除分块（先把 mod 拆成 \(k-i\lfloor k/i\rfloor\)） |
| \(\sum\sum[\gcd(a,b)=1]\) | 莫比乌斯：\(\sum\mu(d)\lfloor N/d\rfloor\lfloor M/d\rfloor\) |
| \(\sum\sum\gcd(a,b)\) | 欧拉：\(\sum\varphi(d)\lfloor N/d\rfloor\lfloor M/d\rfloor\) |
| \(T\) 大 + \(n\) 大的素性判定 | Miller-Rabin，确定性底数 |
| \(n\ge10^{14}\) 要分解 | Pollard-Rho |
| 「\(\sum_{d\mid n}\) 型恒等式」推导 | 狄利克雷卷积语言 |
| 「选一半多一点使两个和都过半」 | 排序 + 配对贪心（BISHI40） |