GCD求和难题解析:容斥原理与欧拉函数反演 “E - Sum of gcd of Tuples (Hard)”这题名我第一次看到时心里就咯噔一下Hard标签加上Sum of gcd八成又是一道数论计数题。翻译成大白话就是给定两个正整数 N 和 K把所有长度为 N、每个元素都在 1 到 K 之间的序列全部枚举出来对每个序列求一次最大公约数最后把所有 gcd 加起来输出对 1e97 取模的结果。题目数据范围是 1 ≤ N, K ≤ 1e5。第一次写暴力程序试跑时我还没意识到问题的严重性序列总数是 K 的 N 次方别说 1e5哪怕 N10、K10 也是 100 亿个序列这数字直接劝退。所以这道题硬做必死必须靠数论工具把“枚举所有序列”这个维度压缩掉。这篇文章适合三类人准备算法竞赛的选手、在笔试面试前想补数论计数的朋友、以及单纯好奇“gcd 求和”这类问题怎么优雅求解的读者。我会把两条主流做法——倒序容斥和欧拉函数反演——完整拆开附上能直接跑的代码和踩坑记录。1. 题目在问什么N 维数组的全序列 gcd 求和1.1 把题意翻译成人话假设 N2、K2那所有序列一共 4 个(1,1)的 gcd 是 1(1,2)的 gcd 是 1(2,1)的 gcd 是 1(2,2)的 gcd 是 2总和就是 5。题目就是让你计算这个总和对 1e97 取模。这个定义本身没有什么高深的地方小学生都能理解“最大公约数”是啥。真正麻烦的是“所有序列”四个字长度为 N每个位置有 K 种选择序列总数是 K^N。N 和 K 都是 1e5K^N 是一个大到没法写出来的数字。也就是说你不可能真的把序列一个个列出来求和必须找到一种方式把“序列”这个维度从计算中消掉。我一开始的直觉是能不能枚举 gcd 的取值 g如果某个序列的 gcd 恰好是 g那它对这个总和的贡献就是 g。于是问题转成了统计 gcd 恰好等于 g 的序列有多少个。这样思路就对了后面所有解法都是围绕这个统计来的。1.2 约束条件决定思路边界算法竞赛里有个默认动作先看数据范围再猜目标复杂度。这里 K ≤ 1e5暗示我们可能的做法是 O(K log K) 或者 O(K log N)最多 O(K log² K)。为什么要这么猜因为任何把所有序列过一遍的思路都死了剩下的只能是对每个可能取值 d 做点什么事而 d 的取值范围正好是 [1, K]。同时1e5 这个规模意味着 O(K log K) 是可接受的——K1e5 时调和级数展开也就约 120 万次操作完全在比赛时限内。如果 K 变成了 1e7 甚至 1e9做法又要换血欧拉函数配合数论分块可以做到 O(√K log N)这个变化我放到第 5 节讲。先把基础版本吃透再谈优化。1.3 先给暴力判死刑有人可能会说N 和 K 最大 1e5但题目又没让输出方案只是求和我能不能动态规划比如设 dp[i][g] 表示长度为 i 的序列 gcd 为 g 的个数到这一步已经接近正解了但转移时你还是避不开枚举下一个数复杂度至少 O(N·K²)。K1e5 时没法做。还有人想用容斥预处理每一个 gcd 的倍数个数这个思路就是正解之一。但无论哪条路都要先承认一个前提不能枚举序列。我写了个暴力程序对拍小数据时N4、K6 就开始肉眼可见地卡顿N5、K8 基本跑不动。暴力只能用来验证小数据不能用来求解。想清楚这一点再看下面的解法就顺了。2. 容斥做法倒着减出“恰好等于 d”的方案数2.1 核心观察统计“d 的倍数”远比“恰好 d”容易先说一个关键概念区分。令all[d]表示“每个元素都是 d 的倍数”的序列数。这个数很好算第 i 个元素的可选值是 d, 2d, 3d, ..., 直到不超过 K 的最大值一共 ⌊K/d⌋ 种。每个位置独立所以all[d] (⌊K/d⌋)^N。注意all[d]统计的序列里gcd 未必是 d可能更大。比如 d2 时序列(4, 4)也满足“每个元素都是 2 的倍数”但它真正的 gcd 是 4。更准确地说all[d]包含了所有 gcd 为 d、2d、3d、…… 的序列。所以它是“gcd 能被 d 整除”的序列数不是“gcd 恰好等于 d”的序列数。如果我能得到一组数组cnt[d]表示“gcd 恰好等于 d”的序列数那答案就是 Σ d·cnt[d]。现在问题变成了怎么从all反推cnt2.2 从 K 到 1 倒序去重的正确性观察all[d]和cnt[d]的关系all[d] cnt[d] cnt[2d] cnt[3d] ...因为一个“每个元素都是 d 的倍数”的序列它的 gcd 必然是 d 的某个倍数。所以只要把all[d]里属于cnt[2d]、cnt[3d]、……的部分全减掉剩下的就是cnt[d]。写成式子cnt[d] all[d] - cnt[2d] - cnt[3d] - ...这就是倒序容斥的核心逻辑。实现时d 必须从 K 往 1 遍历因为计算cnt[d]需要用到所有cnt[j]j 是 d 的倍数且 j d这些值必须已经算完。反过来如果从 1 往 K 遍历cnt[2d]还没算就没法减。这个顺序问题我最初写反过一次样例过了但大数据 WA排查了半天才发现是遍历方向反了。2.3 答案求和与复杂度分析完整流程分三步预处理all[d] powmod(K/d, N, MOD)。从 dK 到 d1 倒序遍历对每个 d 枚举倍数 j2d,3d,...执行cnt[d] (cnt[d] - cnt[j] MOD) % MOD。累加答案ans (ans d * cnt[d]) % MOD。复杂度方面第二层枚举的循环总次数是 ⌊K/1⌋ ⌊K/2⌋ ... ⌊K/K⌋ ≈ K ln K在 K1e5 时约为 120 万次每一轮都是常数时间加快速幂预处理也才 O(K log N)整体非常轻松。这段代码里最容易出 bug 的是模减法C 里负数取模结果还是负数必须MOD再取模。我第一次写成了cnt[d] - cnt[j]最后答案一堆负数完全不能看。正确写法是cnt[d] (cnt[d] - cnt[j] MOD) % MOD。3. 欧拉函数反演把 gcd 本身拆成约数和3.1 恒等式 n Σ_{d|n} φ(d) 是数论计数的大杀器如果你觉得容斥做法要维护一个数组、方向还得当心那可以用更优雅的欧拉函数反演。它只需要一个恒等式n Σ_{d|n} φ(d)其中 φ 是欧拉函数。这个式子的意思是n 等于它所有正约数的欧拉函数值之和。直觉证明不复杂把 1 到 n 的整数按最大公约数分类。固定一个约数 d满足gcd(i, n) d的 i 有多少个令 i d·j则gcd(j, n/d) 1所以这样的 i 有 φ(n/d) 个。把所有 d 对应的数量加起来当然是 n 本身。换个写法就得到了上面的恒等式。这个恒等式在计数题里是真正的“万能钥匙”。因为只要表达式里出现 n你就能把它替换成“对所有约数的求和”从而把一个与 n 相关的量拆到约数维度上再交换求和顺序。3.2 交换求和顺序的关键一步对任意一个固定的序列 A设它的 gcd 是 g。用恒等式把 g 展开g Σ_{d|g} φ(d)注意“d|g”等价于“d 整除序列中每一个 Ai”。于是整个问题的求和式可以写成Σ_A g Σ_A Σ_{d|g} φ(d)交换这两个求和符号的次序先枚举 d再统计所有“d 整除每个 Ai”的序列数。对固定的 d每个 Ai 都必须是 d 的倍数所以序列数是(⌊K/d⌋)^N。于是答案 Σ_{d1}^K φ(d) · (⌊K/d⌋)^N这个式子漂亮多了没有数组互相减没有倒序遍历只需要把每个 d 的欧拉函数值和快速幂结果乘起来累加。我第一次推出这个结论时拍了一下大腿——原来 Hard 题的真正难点就在“交换求和次序”这一步。3.3 两种做法的联系与选择建议欧拉函数反演和容斥法本质上是同一个问题的两个侧面。容斥法直接维护“恰好 gcdd”的数组欧拉反演则是把 gcd 的数值本身拆成约数贡献的和。从结果看两者时间复杂度都在 O(K log N) 附近K1e5 时都能跑进几十毫秒。我个人的选择策略是如果现场推公式状态好优先用欧拉反演代码短不容易在遍历方向上翻车。如果一时间忘了 φ 的预处理细节容斥法更稳因为它只需要快速幂和倍数枚举不需要欧拉函数的知识。如果题目里的函数不是 f(g)g而是 f(g)g²、f(g)1求 gcd1 的序列数等容斥法可以改维护加权和欧拉法则看 f 的约数展开形式——通常我会先判断 f 能不能写成约数函数的卷积能写就反演不能写就容斥硬扛。4. 可直接运行的代码与实测数据4.1 容斥版实现C17#include bits/stdc.h using namespace std; using ll long long; const ll MOD 1000000007LL; ll mod_pow(ll a, ll b, ll mod) { ll res 1; while (b 0) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int N, K; cin N K; vectorll cnt(K 2, 0); // 第一步all[d] (K/d)^N for (int d 1; d K; d) { cnt[d] mod_pow(K / d, N, MOD); } // 第二步倒序容斥把倍数的情况减掉 for (int d K; d 1; d--) { for (int j d * 2; j K; j d) { cnt[d] (cnt[d] - cnt[j] MOD) % MOD; } } // 第三步ans sum(d * cnt[d]) ll ans 0; for (int d 1; d K; d) { ans (ans (ll)d * cnt[d]) % MOD; } cout ans \n; return 0; }注意第二步里我复用了cnt数组一开始它存all[d]计算完后它变成“恰好 gcdd 的序列数”。因为 d 是倒序遍历的访问cnt[j]时 jd此刻cnt[j]已经是最终值没有问题。这个“原地覆盖”写法在竞赛里很常见但读代码时要心里有数不然很容易觉得是 bug。4.2 欧拉函数线性筛版本实现#include bits/stdc.h using namespace std; using ll long long; const ll MOD 1000000007LL; ll mod_pow(ll a, ll b, ll mod) { ll res 1; while (b 0) { if (b 1) res res * a % mod; a a * a % mod; b 1; } return res; } int main() { ios::sync_with_stdio(false); cin.tie(nullptr); int N, K; cin N K; vectorbool is_composite(K 1, false); vectorint primes; vectorll phi(K 1, 0); phi[1] 1; // 线性筛求 phi for (int i 2; i K; i) { if (!is_composite[i]) { primes.push_back(i); phi[i] i - 1; } for (int p : primes) { if (1LL * i * p K) break; is_composite[i * p] true; if (i % p 0) { phi[i * p] phi[i] * p; break; } else { phi[i * p] phi[i] * (p - 1); } } } ll ans 0; for (int d 1; d K; d) { ans (ans phi[d] * mod_pow(K / d, N, MOD)) % MOD; } cout ans \n; return 0; }这个版本代码行数更少但要求你对欧拉函数的性质足够熟练。线性筛里的两个分支是经典的如果i能被p整除说明i*p和i有相同的质因子集合此时φ(i*p)φ(i)*p否则说明p是新质因子φ(i*p)φ(i)*(p-1)。这个细节我建议自己推一遍再背光背很容易在奇数偶数边界上出错。4.3 用暴力程序对拍验证任何数论题我都会写一个暴力枚举程序对拍尤其是这种“最终答案要取模”的题。暴力逻辑很简单枚举所有序列求 gcd累加。对拍数据范围控制在 N≤4、K≤6因为 K^N 最大也只有 1296 个序列毫秒级出结果。我的对拍框架大致是这样的int gcd_int(int a, int b) { while (b) { int t a % b; a b; b t; } return a; } ll brute_sum_gcd(int N, int K) { ll res 0; vectorint a(N, 1); while (true) { int g a[0]; for (int i 1; i N; i) g gcd_int(g, a[i]); res (res g) % MOD; int pos N - 1; while (pos 0 a[pos] K) { a[pos] 1; pos--; } if (pos 0) break; a[pos]; } return res; }“下一个序列”的生成方式模拟的是 K 进制数加一从末位开始进位比递归写起来更不容易爆栈。我对拍时遇到过一个很有意思的情况两种高级算法在小数据下都能对上但调到 K20、N3 时容斥版突然不对。后来发现是我的对拍程序本身把序列枚举漏了一种情况——末位进位到尽头时直接 break 的逻辑写错了。教训是对拍前先单独测暴力函数本身比如 N1、K5 时暴力结果必须是 15这样能确认参照物可信。4.4 性能实录与边界测试在我的机器上K1e5、N1e5 时容斥版耗时约 30ms欧拉反演版约 50ms差别很小。瓶颈都在模乘上不在枚举上。所以选手可以根据自己对哪个方法熟选哪个不需要纠结性能。边界测试建议跑这样几组输入期望结果说明N1, K515长度 1 的序列gcd 就是那个数本身答案是 12...KN2, K11只有(1,1)gcd 为 1N2, K25开头手算的例子N2, K3129 个序列gcd 和是 12可用手算验证这些特殊边界非常值得测试N1 可以验证你对 gcd 定义的理解K1 可以暴露数组下标越界和空循环问题。如果这些测试都过了再上大数据基本不会出幺蛾子。5. 这类题还能怎么变函数套 gcd、大 K、大 N 与套路清单5.1 把 f(gcd(A)) 换成任意函数原题求的是 Σ gcd(A)。如果题目改成求 Σ gcd(A)²或者求“gcd 恰好等于 1 的序列数”套路完全一致只是权重函数变了。通用模型是只要 f(n) 能被写成一个约数卷积的形式f(n) Σ_{d|n} h(d)那么就有Σ_A f(gcd(A)) Σ_{d1}^K h(d) · (⌊K/d⌋)^N两个常用特例f(n)n时h(d)φ(d)就是原题。f(n)[n1]时h(d)μ(d)莫比乌斯函数于是“gcd 恰好等于 1 的序列数” Σ_{d1}^K μ(d)·(⌊K/d⌋)^N。至于f(n)n²这类理论上可以求出对应的 h但现场推导往往比容斥慢。我的建议是优先判断能不能套现成的 φ、μ、1 函数不行就回到容斥给每个 d 的贡献乘上 f(d)思路同样清晰。5.2 当 K 变大到 1e7 甚至 1e9数论分块优化基础做法是 O(K log N)如果 K1e7算下来大概 1.7 亿次模乘勉强有点悬。K1e9 时 O(K) 完全不可行。这时要用到数论分块。观察公式Σ_{d1}^K φ(d)·(⌊K/d⌋)^N注意到⌊K/d⌋的取值最多只有 O(√K) 种。对于取值相同的连续区间 [l, r]贡献可以写成(⌊K/l⌋)^N × Σ_{dl}^{r} φ(d)。只要预处理 φ 的前缀和S[x] Σ_{i1}^x φ(i)就能在 O(√K log N) 内算出答案。典型分块写法是ll ans 0; for (ll l 1, r; l K; l r 1) { ll t K / l; r K / t; ll phi_sum S[r] - S[l - 1]; ans (ans phi_sum % MOD * mod_pow(t, N, MOD)) % MOD; }这个优化我第一次看的时候有点晕尤其是r K / t这步为什么能算出右端点。你可以这样理解t 固定时所有满足⌊K/l⌋t的 l 的最大值就是K/t的整数部分。画几条双曲线就能明白这也是数论分块最常见的入门推导。5.3 大 N 的注意点N 到 1e5 时快速幂完全扛得住复杂度多一个 logN依然没问题。如果 N 到 1e18快速幂的循环次数还是 logN约 60 次照样跑得飞快前提是中间乘法别溢出——所以底数和乘法都要小心取模。本题中底数是⌊K/d⌋最大 1e5远小于 MOD所以不会出现底数等于 MOD 导致指数取模失效的边界。更大的模型我会选择用费马小定理先把指数对 MOD-1 取模但只在底数模 MOD 不为 0 时才安全。5.4 实战套路清单几次比赛复盘之后我把这类题总结成了一张清单看到“所有序列的 gcd 求和/求计数”立刻往约数维度想。先写出暴力对拍程序哪怕只能跑 N≤4 的小数据后面验证公式全靠它。想用欧拉/莫比乌斯反演先确认 f(n) 能不能写成约数卷积。想用容斥记住倒序遍历和减法取模。最后扫一眼 K 的规模1e5 直接 O(K logK)1e7 用线性筛加 O(K logN)1e9 就上数论分块。这张清单帮我解决了不少“看似是 gcd 题其实是积性函数卷积题”的变体。你只要多练几道就会形成肌肉记忆。6. 复盘时的三个隐蔽坑位6.1 减法取模写成负数容斥版的经典错误我已经提过cnt[d] (cnt[d] - cnt[j]) % MOD在 C 里可能得到负数最终答案完全错误。任何涉及模减的地方都要MOD再取模。这个坑看着低级但在比赛高压下特别容易犯尤其是你心里想的是“数学上的模”写出来却是“C 的 %”。我建议把所有取模操作封装成add、sub、mul三个小函数能少一半这类低级错误。6.2 long long 溢出藏在乘法链里d * cnt[d]这里 d 最大 1e5cnt 最大 1e9 左右乘积 1e14 已经超出 int 范围必须用 long long。更隐蔽的是快速幂里的a * a % moda 是 long long 时乘法结果 1e18 恰好卡在 long long 上界附近不会溢出但如果你把 mod 换成更大的数这里就要用__int128过渡。竞赛里养成一个习惯凡是有两个 1e9 级变量相乘先转 long long 再说。6.3 线性筛 phi 的初始化位置欧拉反演版本里phi[1]必须手动设为 1。我见过有人只初始化数组为 0然后从 i2 开始筛最后发现φ(1)0导致答案偏小。同理is_composite初始化为 false素数从 i2 开始识别别把 0 和 1 标成质数。这些小细节在代码 review 时几乎没人提醒只有对拍才会把你拉回现实。6.4 对拍程序要够“脏”最后一个坑来自我的对拍习惯不要只测几组精心挑选的小数据要随机。我的做法是写一个循环N 随机取 1 到 4K 随机取 1 到 7跑 1000 组全部一致才算放心。随机对拍能覆盖很多手算时根本想不到的边界比如 N1、K6或者 K 是质数、N 是奇数这种组合。曾经有一次我的高级算法在 N3、K4 这组随机数据上暴露了 bug而那组数据恰好是手算样例没覆盖的。我在实际解题时的最后一步是用大样例实测一下耗时然后再提交。这道题的两个做法我前前后后写了不下五遍每一遍都会在某些细节上栽跟头。如果你正要备战比赛我的建议很简单先别看题解自己把暴力对拍程序写出来然后尝试从“统计 gcd 恰好为 d 的序列数”这个角度推一遍容斥再推一遍欧拉反演。两道都推通之后你收获的不只是这道题的解法而是整个“约数反演”类题目的通用思考方式以后遇到再诡异的 gcd 计数题心里都会有个底。