莫比乌斯反演与数论分块：从互质数对到杜教筛实战

第一次真正被莫比乌斯函数卡住是在一个求互质数对个数的题上给定 n 和 m求 1≤i≤n、1≤j≤m 中 gcd(i,j)1 的数对数量。我老老实实写了个双重循环暴力小数据过得挺开心一交上去 T 得连输出的机会都没有。后来有人丢过来一行式子说预处理一下莫比乌斯函数加个整除分块就完了。那一刻我才意识到数论题里有一道看不见的分水岭分水岭之前比的是代码能力分水岭之后比的是在纸上把求和式推干净的能力。莫比乌斯函数、莫比乌斯反演就是这条分水岭上最关键的两块石头。这篇东西我打算按自己这几年刷题的顺序来写先把莫比乌斯函数本身的定义和性质捋顺再把反演的两种形式讲明白然后从最简单的 gcd 计数一路推到YY的GCD那种需要换元加预处理的模型最后补上多组询问、数据范围爆炸时用到的杜教筛以及我自己踩过的一堆坑。全文的代码以 C 为主思路是通用的。如果你刚学完数论分块、欧拉筛和狄利克雷卷积这篇基本能把你从看得懂题解带到自己能推出来如果你已经会写套路题中间关于换元和预处理那几段也许能帮你省下几次重写的时间。1. 莫比乌斯函数到底是个什么东西1.1 从容斥讲起它的定义其实很好记莫比乌斯函数通常记作 μ(n)定义只有三条但每一条都有明确的用途。第一μ(1)1这是整个体系的锚点所有反演公式的最后一步都要靠它收口。第二如果 n 含有平方因子也就是存在某个素数 p 使得 p² 整除 n那么 μ(n)0。第三如果 n 是 k 个互不相同的素数的乘积也就是 np₁p₂…p_k那么 μ(n)(-1)^k。换句话说μ 的取值只有 -1、0、1 三种值域极小这一点非常重要因为它意味着前缀和、卷积这类操作的开销都可以压得很低。我第一次看到这个定义时的第一反应是这也太随意了吧。但你把几条放在一起看就会发现它其实是在描述n 的素因子集合这个对象平方因子相当于集合里出现了重复元素直接判零剩下干净的情况就看集合大小是奇是偶。用组合数学的话说μ(n) 就是在给 n 的因子集合做一次带符号的计数。为什么需要带符号因为后面所有反演的推导本质都是容斥原理的代数化表达而容斥的符号正是由多加一次还是少减一次决定的。这里有个容易被忽略的细节μ 是积性函数。也就是说当 gcd(a,b)1 时μ(ab)μ(a)μ(b)。这条性质是它能被线性筛的核心原因。积性函数有一整套通用筛法只要你能在i 与最小质因子互质和i 能被最小质因子整除两种情况下分别用已求出的值推出新值就能 O(n) 地把整张表求出来。μ 的两种情况特别好推互质时直接乘 -1被整除时因为出现了重复素因子而直接归零。1.2 那条最关键的求和性质真正让 μ 从一个奇怪的函数变成解决问题的工具的是下面这条恒等式∑_{d | n} μ(d) [n 1]右边的方括号表示当 n 等于 1 时取 1否则取 0。也就是说n 的所有因子的 μ 值加起来除了 n1 之外全部抵消干净。我建议你一定要亲手验证几个数n1 时只有 d1和为 1n2 时 d 取 1、2μ(1)μ(2)1-10n6 时 d 取 1、2、3、6和为 1-1-110n4 时 d 取 1、2、4和为 1-100。验完这四个你基本就能体会到它的作用了。这条性质的用途在于把等于 1变成一个可以展开的求和。在数论题里我们经常要处理 gcd(i,j)1 这种条件直接统计很难但如果用这个恒等式把它换成对 gcd 的所有因子求和那求和号就能交换位置内层变成简单的整除计数复杂度立刻降下来。这就是莫比乌斯反演最核心的一次思维转换后面的所有题目都是在这个骨架上加装饰。顺便说一句这条性质用狄利克雷卷积的写法会非常简洁μ * 1 ε其中 1 是恒为 1 的常函数ε 是只在 1 处取 1 的单位函数。如果你已经熟悉卷积那一套记号用这个写法推导会快很多如果不熟用上边的展开式慢慢算也完全没问题结果是一样的。1.3 三种常见写法与线性筛求法μ 的求法我见过三种写法各有用武之地。第一种是暴力质因数分解对每个数单独分解复杂度 O(n√n)只适合 n 很小或者只要求单个值的时候用。第二种是类似埃氏筛的做法枚举每个素数把它的倍数对应位置的值更新一遍具体做法是先把所有 μ 初始化为 1对每个素数 p遍历其倍数把 μ 乘上 -1同时对 p² 的倍数直接置零。这个做法复杂度 O(n log log n)写起来短比赛时如果时间紧可以直接用。第三种就是线性筛复杂度严格 O(n)而且顺带能筛出一堆别的积性函数是我最常用的写法。线性筛的代码结构几乎是模板化的我贴一份后面所有题都基于这个const int N 10000005; int primes[N / 10], cnt; int mu[N]; bool vis[N]; void sieve(int n) { mu[1] 1; for (int i 2; i n; i) { if (!vis[i]) { primes[cnt] i; mu[i] -1; // 素数的 μ 值为 -1 } for (int j 1; j cnt 1LL * i * primes[j] n; j) { vis[i * primes[j]] true; if (i % primes[j] 0) { mu[i * primes[j]] 0; // 出现平方因子 break; // 保证每个合数只被最小质因子筛掉 } mu[i * primes[j]] -mu[i]; // 新增一个不同的素因子符号翻转 } } }这段代码里有几个点值得说清楚。mu[i * primes[j]] 0必须写在break之前因为如果先跳出循环再赋值就漏掉了。break的条件i % primes[j] 0是线性筛的精髓它保证每个合数只被它的最小质因子筛到一次从而把复杂度压到 O(n)。另外 1LL 那个乘法是防止 i 和 primes[j] 都很大时 int 溢出虽然通常 N 在 1e7 量级不会出问题但养成习惯没坏处。提示如果你的数据范围是 1e7 甚至更大vis用 bool 数组会更省内存但如果要反复清空用 vector 动态开可能更灵活。我一直用的是全局静态数组因为 1e7 的 bool 只有 10MB 左右完全可以接受。2. 莫比乌斯反演的两种面孔2.1 约数形式F(n) ∑_{d|n} f(d)莫比乌斯反演的标准叙述有两种第一种叫约数形式如果两个函数满足F(n) ∑_{d | n} f(d)那么反过来就有f(n) ∑_{d | n} μ(d) · F(n / d)用卷积写就是 F f * 1所以 f F * μ。这个方向的推导非常直接把 F(n/d) 按定义展开交换求和顺序里面正好冒出 ∑ μ 的那条恒等式等于把不需要的项全部消掉。我建议你在纸上把这一步走一遍走完之后你就理解了为什么反演这个名字其实是逆运算的意思——1 这个函数在卷积意义下是可逆的而它的逆恰好就是 μ。这条形式在实际题目里出现的频率很高尤其是在已知每个数的倍数的某个统计量反推每个数本身的量这类问题中。比如某些计数问题会让你先求所有倍数的和再反推出恰好等于某个数的方案数那就是约数形式的直接应用。2.2 倍数形式F(n) ∑_{n|d} f(d)第二种叫倍数形式也是我在刷题时用得最多的一种。它的表述是如果F(n) ∑_{n | d} f(d)那么f(n) ∑_{n | d} μ(d / n) · F(d)注意这里所有求和都是对 n 的倍数进行的。这个形式和前面的约数形式看起来只是把整除方向翻转了一下但实际使用场景差别很大。约数形式适合从因子往上传的问题倍数形式适合从倍数往下推的问题。为什么我说倍数形式更常用因为绝大多数 gcd 计数题的推导路径是这样的先把条件 gcd(i,j)k 转成 gcd(i/k, j/k)1再用 μ 的求和性质把等于 1展开成对因子的求和。展开之后求和变量 d 是 gcd 的因子而 d 的倍数恰好对应着能被 d 整除的那些 i、j。你会发现这整个过程天然地和倍数关系绑定在一起所以倍数形式用起来更顺手。这也是为什么很多题解直接写[gcd(i,j)1] ∑_{d|gcd(i,j)} μ(d)而不去套用公式化的反演叙述——它是同一条性质的另一种表达。2.3 为什么我更推荐用卷积和指标替换来理解学到这里我想说一个个人观点反演那两条公式其实不用死记因为记错了也很容易自查出来。我更推荐你把注意力放在两件事上。第一件是狄利克雷卷积的基本对应关系μ * 1 εφ * 1 id也就是 ∑_{d|n} φ(d) nid * μ φ。这三条撑起了数论函数的大半壁江山记住它们比记住两个反演公式有用得多。第二件是指标替换这个操作看到 [条件] 这种方括号就想办法把它换成求和换成求和之后就能交换顺序看到求和就想办法找整除结构找到整除结构就能分块加速。这两件事一旦形成条件反射你在纸上推式子会快很多。我自己的经验是遇到新题先问三个问题这个式子里的恰好等于多少能不能换成倍数求和换完之后 d 的范围和 n、m 是什么关系交换求和顺序之后内层能不能写成某个数被多少个 i 整除这种简单形式这三个问题答完大多数题目的方向就定了。剩下的是技术细节比如预处理什么、分块怎么写。需要提醒的是反演公式本身对 f 和 F 的定义域有隐含要求一般默认定义在正整数上且求和是有限的。在竞赛场景里这通常不是问题因为 n、m 都是有界的。但如果你在论文或工程里用最好确认一下收敛性条件别把有限和的直觉直接搬到无限和上。3. 从零推一道经典题互质数对计数3.1 问题的引入和朴素做法先把题目定清楚给定 T 组询问每组给 n 和 m求∑_{i1}^{n} ∑_{j1}^{m} [gcd(i, j) 1]如果 n、m 只有几百双重循环加一个 gcd 就完事了复杂度 O(nm log)。但题目里 n、m 通常给到 1e5 甚至 1e7询问组数上千那暴力必然超时。这时候需要的就是把统计问题变成数论函数求和问题。3.2 关键转换[gcd 1] ∑ μ(d)整个推导的第一步是把条件改写[gcd(i, j) 1] ∑_{d | gcd(i, j)} μ(d)这一步的依据就是 1.2 节那条恒等式把 n 换成 gcd(i,j) 即可。改写之后原式变成∑_{i1}^{n} ∑_{j1}^{m} ∑_{d | gcd(i, j)} μ(d)第二步是交换求和顺序。注意 d 整除 gcd(i,j) 等价于 d 同时整除 i 和 j所以d | gcd(i,j)这个条件可以直接拆成d | i 且 d | j。这样一来求和可以重排成先枚举 d∑_{d1}^{min(n,m)} μ(d) · (∑_{i1}^{n} [d | i]) · (∑_{j1}^{m} [d | j])内层的两个求和是纯计数1 到 n 里 d 的倍数有 ⌊n/d 个于是最终形式清爽得让人舒服∑_{d1}^{min(n,m)} μ(d) · ⌊n / d⌋ · ⌊m / d⌋这就是那道题的核心式子。你会发现整个过程没有用到任何高深的技巧只有两步把条件换成求和把求和换成计数。我第一次推出来的时候有一种原来就这么回事的感觉但真到自己上手推卡壳是常事因为第二步的交换顺序需要你对整除关系和求和范围都足够熟悉。3.3 整除分块加速式子推出来了但直接算还是 O(min(n,m)) 一次询问组数一多照样超时。这时候要用整除分块也叫数论分块。它的依据是⌊n/d⌋ 这个函数随着 d 增大是分段常数的只有 O(√n) 个不同的取值。所以我们可以把取值相同的 d 归成一段一段一段地算。具体做法是当前左端点 l找到最大的 r 使得 ⌊n/l⌋ ⌊n/r⌋ 且 ⌊m/l⌋ ⌊m/r⌋这个 r 就是r min(n / (n / l), m / (m / l))一段之内的 ⌊n/d⌋ 和 ⌊m/d 都不变只有 μ 的前缀和在变所以只要维护 μ 的前缀和一段的贡献就是(sum[r] - sum[l-1]) * (n/l) * (m/l)然后跳到 r1。整体复杂度从 O(n) 降到 O(√n) 每组询问。配合 O(n) 的线性筛预处理总复杂度约 O(n T√n)1e5 的数据范围下非常宽裕。这里有个必须注意的点写r min(n/(n/l), m/(m/l))之前一定要保证 n ≤ m否则如果 n m当 d 超过 m 时 ⌊m/d⌋ 会变成 0虽然数学上没问题但代码里的除法步骤要小心而且循环上界应该取 min(n,m)不然会多算。我一般的写法是进函数先if (n m) swap(n, m);然后循环for (int l 1, r; l n; l r 1)。3.4 代码实现与复杂度把上面几步拼起来完整实现大概是这样int mu[N], sum[N]; // sum 为 μ 的前缀和 // 先调用 sieve(N) 求 mu再求前缀和 long long solve(int n, int m) { if (n m) swap(n, m); long long ans 0; for (int l 1, r; l n; l r 1) { r min(n / (n / l), m / (m / l)); ans 1LL * (sum[r] - sum[l - 1]) * (n / l) * (m / l); } return ans; }前缀和的预处理就是for (int i 1; i N; i) sum[i] sum[i - 1] mu[i];注意sum要用 int 就够因为 μ 只取 -1、0、1范围在 [−N, N] 内N 到 1e7 也不会溢出 int。但ans必须用 long long因为(n/l) * (m/l) * 区间长度这个量级在 n、m 取 1e7 时能到 1e21早就超过 int 了。这个坑我踩过不止一次程序本地测试小数据全对一交大数据就 WA查了半天才想起来是溢出。注意整除分块写法里n / (n / l)这个除法如果忘记先把 n、m 排序或者循环上界写成了 m都会出现除零或者越界。这是我在帮别人 debug 时见得最多的一类错误务必留意。4. 进阶模型换元之后的预处理技巧4.1 从前置问题到gcd 是质数把基础题做熟之后很快就会碰到它的变体最经典的一个是求∑_{i1}^{n} ∑_{j1}^{m} [gcd(i, j) 是质数]这道题比上一道难在是质数这个条件不是一个简单的整除判定而是对 gcd 的取值本身提要求。直接套前面的做法会发现 μ 那一步没问题但多出来的质数约束没办法直接吸收进求和式里于是需要做一次换元。第一步的换元思路是枚举质数 p把所有 gcd 等于 p 的贡献单独拎出来。因为 gcd(i,j)p 等价于 gcd(i/p, j/p)1 且 p 同时整除 i、j所以可以写成∑_{p 为质数} ∑_{i1}^{n/p⌋} ∑_{j1}^{⌊m/p} [gcd(i, j) 1]内层就是上一节已经解决的基础问题代入它的结论∑_{p 为质数} ∑_{d1}^{⌊min(n,m)/p⌋} μ(d) · ⌊n/(pd)⌋ · m/(pd)4.2 换元 T p·d到了这一步两层求和里出现了 p 和 d 的乘积形式不统一没法直接用整除分块。常规做法是令 T p·d把求和的主变量换成 T∑_{T1}^{min(n,m)} ⌊n/T⌋ · ⌊m/T⌋ · (∑_{p | T, p 为质数} μ(T/p))后半部分那个括号里的东西你可以把它看成是一个新的数论函数 f(T) ∑_{p|T} μ(T/p)其中求和只对 T 的质因子进行。一旦把它单独定义出来整个式子就回到了和基础题完全一样的形式只是 μ(d) 被换成了 f(T)。这意味着整除分块那一套可以原封不动地搬过来只要提前把 f 的前缀和算好。这就是这道题的全部难点换元 识别出可预处理的部分。4.3 f 数组的预处理技巧f 的计算有一个很漂亮的技巧用枚举质数的方式做// primes 是筛出来的质数表cnt 是数量 for (int j 1; j cnt; j) { int p primes[j]; for (int i 1; 1LL * i * p N; i) { f[i * p] mu[i]; } }这段代码的意思很直白对每个质数 p遍历它的所有倍数 ip把 μ(i) 累加到 f[ip] 上。这样枚举下来f[T] 恰好收集到了所有满足 T p·i 的项也就是所有质因子对应的 μ(T/p)。复杂度是多少每个质数枚举 N/p 次总数是 N·∑(1/p)约等于 N ln ln N在 1e7 的量级下也就是几千万次加法完全可以接受。接下来把 f 求前缀和然后在整除分块里替换掉原来 μ 前缀和的位置代码几乎不用改long long solve(int n, int m) { if (n m) swap(n, m); long long ans 0; for (int l 1, r; l n; l r 1) { r min(n / (n / l), m / (m / l)); ans 1LL * (sf[r] - sf[l - 1]) * (n / l) * (m / l); } return ans; }前后对比一下你会发现所谓的进阶题很多时候只是把分母上的一个简单函数换成了一个需要额外预处理的函数。真正的工作量在推式子和设计预处理方案上分块那一层反而是复用的。提示f 的值可能为负数因为它来自 μ 的累加所以sf前缀和数组要用 long long 或者确认不会溢出。另外 f[T] 的具体数值并不重要重要的是它的前缀和能分段累加所以不需要做任何归一化处理。5. 数据范围爆炸时杜教筛求 μ 前缀和5.1 从线性筛的局限说起线性筛的时间复杂度是 O(N)空间也是 O(N)当 N 在 1e7 以内时毫无压力。但当题目把 n 给到 1e9、1e10 甚至更大时我们既没法筛出这么长的表也没法存下来。这种情况下整除分块里的那一段区间可能长达数十亿而每段需要的是 μ 的区间和也就是前缀和。如果能快速算出任意一个 S(x) ∑_{i1}^{x} μ(i)问题就解决了。杜教筛就是干这个的。5.2 杜教筛的推导与记忆化实现推导的出发点还是那条恒等式 ∑_{d|n} μ(d) [n1]。对 1 到 n 求和∑_{k1}^{n} ∑_{d|k} μ(d) 1左边交换求和顺序枚举 d 和它的倍数 k d·m∑_{d1}^{n} μ(d) · ⌊n/d⌋ 1再把 ⌊n/d⌋ 展开成计数⌊n/d⌋ ∑_{m1}^{⌊n/d⌋} 1代入后∑_{m1}^{n} ∑_{d1}^{n/m⌋} μ(d) ∑_{m1}^{n} S(⌊n/m⌋) 1把 m1 那一项单独拎出来它就是 S(n)于是S(n) 1 - ∑_{m2}^{n} S(⌊n/m⌋)这个递推式就是杜教筛的全部内容。它的妙处在于右边求和里的 ⌊n/m 只有 O(√n) 种不同取值所以可以把它们分块每块递归地算一次 S 的值。为了不重复计算用一个哈希表或者数组把已经算出的 S 值存起来。实际实现里通常有一个阈值比如 n ≤ 1e6 直接查线性筛出来的前缀和超过阈值才走递归加记忆化。unordered_mapint, long long memo; long long S(int n) { if (n LIMIT) return pre[n]; // pre 是线性筛的前缀和 auto it memo.find(n); if (it ! memo.end()) return it-second; long long ans 1; for (int l 2, r; l n; l r 1) { r n / (n / l); ans - 1LL * (r - l 1) * S(n / l); } memo[n] ans; return ans; }实测下来这样一个记忆化加递归的版本在 n 到 1e10 时算一次的总调用次数在 1e5 量级非常快。需要注意的是memo用哈希表会有常数开销如果追求更快可以用两个数组分别处理小值和大值的映射或者直接开一个大小 1e5 左右的桶加手写哈希。如果不是被卡常数卡到极限unordered_map 完全够用。注意杜教筛的递推式里ans初始值是 1减号对应的是把 m1 那一项移过去。我第一次写的时候把符号弄反了本地点几个小数据还看不出问题因为小值走的是pre分支只有大数据才暴露。写完之后一定要用几个中等大小的值对比线性筛结果别只测小数据。6. 常见问题与排查技巧实录6.1 高频错误速查表我把自己和周围人踩过的坑整理成了一张表基本上每次有人问为什么我答案不对原因都在这几条里症状可能原因排查方法小数据对、大数据错答案变量或前缀和用了 int 溢出检查 ans 是否 long long估算最大量级整除分块死循环或除零没先保证 n ≤ m或 r 计算用了 min 之外的上界进函数先 swap循环上界取 min部分询问答案偏小预处理上界没取到最大 n数组开小了先读入所有询问取最大值再筛到该值倍数形式的 f 数组整体偏大枚举质数时循环条件写成了 i ≤ N改为 1LL * i * p ≤ N线性筛结果错乱μ 值写在 break 之后或素数标记漏了逐个数打印 μ与手算对照多组数据超时每组询问都重新筛一遍预处理放循环外只做一次关于预处理上界这一条我再展开说两句。很多题是 T 组询问每组给不同的 n、m如果你的筛法只筛到第一组的 n那后面几组就会越界访问。正确做法是先全部读入取所有 n 的最大值作为筛的上界。这个操作只需要多存一个询问数组代价很小但能避免一类非常隐蔽的错误。6.2 我的实操心得第一推式子的时候一定要写下来不要在脑子里算。我在纸上推的时候会习惯性标注每一步的求和范围尤其是交换顺序之后范围的边界最容易出错。比如原来 d 的上界是 min(n,m)换元之后 T 的上界还是 min(n,m)但内层对 i 的限制变成了 i ≤ n/T这些都要标清楚不然分块的时候很容易把区间取错。第二代码和推式子要分开写。我一般先把式子完整推完确认无误再动手写代码。中间如果发现推不下去就退回去检查是不是条件转换那一步错了。反过来如果代码写出来答案不对先别怀疑代码回头看看式子有没有问题——这个顺序能省很多时间因为式子错了的话代码怎么调都是错的。第三学会用暴力对拍。我的习惯是写两版程序一版是 O(nm) 的暴力一版是优化后的解然后用随机数据对拍几百组。这一步在莫比乌斯反演里特别重要因为推式子涉及大量符号和边界肉眼很难看出错。对拍能抓出 90% 以上的实现错误剩下的 10% 通常是大数据范围和下溢问题。第四关于整除分块的实现我的建议是把它封装成一个函数模板。因为你会发现从最简单的互质数对到复杂的质数 gcd分块那一层的代码几乎一字不改改的只是前缀和数组的名字。把它抽出来之后写新题的时候只要专注于推式子和预处理效率能提高一大截。我自己维护了一个小模板库里面就有block_sum(l, r, func)这类封装用起来很顺手。第五遇到恰好等于这类条件时第一时间想能不能用 μ 展开。这个反射建立起来之后很多题的入手点就自动出现了。反过来说如果题目里的条件是小于等于或者某个范围那要优先考虑转换或者容斥而不是硬套 μ因为 μ 处理的是精确的整除结构。最后分享一个小习惯每次做完一道莫比乌斯反演的题我会把核心式子和换元思路记到笔记里不写代码只写推导的关键几行。积累到十几道之后你会发现题目之间的差别其实很小无非是多乘一个函数、换成质数限制、数据范围加大要杜教筛这几种变体。把这些变体都归类清楚再遇到新题时你的第一反应就不是这题没做过而是它属于哪一类。这个转变我觉得比多刷几十道题都值。