GCD求和进阶：从倒序容斥到欧拉函数O(K)解法

说起E - Sum of gcd of Tuples (Hard)这道题我第一次看到它时其实是在某个面向竞赛选手的在线题库里当时它被放在一场比赛的 E 号位还带着一个 Hard 后缀明摆着是来拉开差距的。题目本身很短给定正整数 N 和 K求所有长度为 N、每个元素都在 1 到 K 之间的整数序列(a1, a2, ..., aN)它们的最大公约数之和是多少。N 和 K 都可能到 1e5也就是说序列总数是 K 的 N 次方这是一个连边都摸不着的大数暴力枚举完全没有可能。这篇文章我要从最直观的思考路径开始先手算一个小样例建立感觉然后给出两种解法第一种是 O(K log K) 的倒序容斥第二种是用 φ 函数在 O(K) 内直接算出答案。两种都会给完整代码、复杂度推导和易错点。如果你正在刷数论计数题或者想彻底搞懂“gcd 求和”这一类套路这篇应该能帮你省下不少时间。1. 先看懂题目从枚举到“能被 d 整除”的转化1.1 手算小样例把题意吃透拿 N2、K3 来做例子。所有长度为 2、元素在 1 到 3 之间的序列一共有 3^29 个我们完全可以手写出来(1,1) 的 gcd 是 1(1,2) 的 gcd 是 1(1,3) 的 gcd 是 1(2,1) 的 gcd 是 1(2,2) 的 gcd 是 2(2,3) 的 gcd 是 1(3,1) 的 gcd 是 1(3,2) 的 gcd 是 1(3,3) 的 gcd 是 3全部加起来是 12。这个例子虽然小但已经暴露出了解这类问题的关键gcd 的取值只能是 1 到 K 之间的整数所以“求和”本质上可以拆成“先统计每个 d 作为最大公约数出现了多少次再乘 d 累加”。比如这里 gcd1 出现了 7 次gcd2 出现 1 次gcd3 出现 1 次答案就是 1×7 2×1 3×1 12。把目标从“直接求和”改成“统计每个 d 的出现次数”事情就顺畅多了。真正要解决的问题只剩下一个如何不枚举序列快速得到“最大公约数恰好等于 d”的序列数量。只要把这个统计出来后面的累加只是顺手的事。1.2 关键一步定义“至少是 d”和“恰好是 d”在 1 到 K 之间d 的倍数一共有floor(K/d)个。如果要求序列里每个数都能被 d 整除那么每一位都有floor(K/d)种选择由乘法原理可知这样的序列数量是F[d] floor(K/d) ^ N我把 F[d] 称为“至少是 d 的序列数”因为“每个数都能被 d 整除”只说明 d 一定能整除最终 gcd并不代表最终 gcd 恰好等于 d。想一想一个序列 (2,4,6)它的每个数都能被 1 和 2 整除但 gcd 是 2不是 4 也不是 6 的什么别的倍数。所以在 F[d] 里会混入那些实际 gcd 是 2d、3d、4d 的序列。这时候就有一个经典思路再设 G[d] 表示“最大公约数恰好等于 d”的序列数。既然一个“每个数都能被 d 整除”的序列其真实 gcd 只可能是 d 的某个倍数那么用分类的思想很容易写出G[d] F[d] - G[2d] - G[3d] - G[4d] - ...直到超过 K右边的每一项都代表“实际上 gcd 是 2d、3d、4d”的序列它们原本都被包含在 F[d] 里但它们的 gcd 不是 d所以要一个一个减掉。最终答案就是所有d × G[d]的和。这个模型就是整个问题的地基后面两种解法都是从这个式子出发的。2. 解法一O(K log K) 的倒序容斥2.1 容斥关系是怎么推出来的从式子G[d] F[d] - G[2d] - G[3d] - ...可以立刻得到一个很重要的实现结论计算 G[d] 需要用到所有比自己大的倍数位置的 G 值。所以枚举 d 的顺序必须从 K 往 1 走而不能从 1 往 K 走。很多初学者第一次写这道题时习惯从小到大枚举结果在算 G[1] 的时候发现 G[2]、G[3] 还是未知数整个程序只能给出错误答案。举一个具体的例子。假设 K6计算顺序是这样的从 d6 开始6 的倍数超过 K 的不存在所以 G[6] F[6] (6/6)^N 1^N 1。d5同理 G[5] F[5]。d42×48 大于 K没有倍数要减G[4] F[4]。d3需要减掉 G[6]所以 G[3] F[3] - G[6]。d2需要减 G[4]、G[6]G[2] F[2] - G[4] - G[6]。d1需要减 G[2]、G[3]、G[4]、G[5]、G[6]。从大到小算每一步需要的右侧值都已经提前算好了非常直观。这里也不需要引入所谓的“容斥系数”对每个 d 直接把所有倍数的“恰好数量”减一遍就行因为一个序列的真实 gcd 是唯一的它只会落在某一个 G[md] 中绝不会被重复减去。2.2 这个“减倍数”的过程为什么不是指数级可能会有读者担心内层循环要从 2d 一直枚举到 K每个 d 都要做一遍看起来不是 O(K²) 吗其实不会。对于 d1要枚举 2,3,4,...,K 大约 K 个倍数对于 d2要枚举 4,6,8,... 大约 K/2 个对于 d3大约 K/3 个。把所有 d 的枚举量加起来是K/2 K/3 K/4 ... ≈ K × (1/2 1/3 ... 1/K) ≈ K log K这是调和级数求和所以总复杂度是 O(K log K)而不是 O(K²)。再加上每组F[d]需要一次快速幂复杂度里再带一个 O(log N)完全能跑 1e5 的数据。2.3 完整C代码与复杂度#include bits/stdc.h using namespace std; using ll long long; const ll MOD 1000000007; ll qpow(ll a, ll b) { 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 F(K 1), G(K 1); for (int d 1; d K; d) { F[d] qpow(K / d, N); } for (int d K; d 1; --d) { ll g F[d]; // 枚举 d 的倍数跳过 d 本身 for (int x d d; x K; x d) { g (g - G[x] MOD) % MOD; } G[d] g; } ll ans 0; for (int d 1; d K; d) { ans (ans d * G[d]) % MOD; } cout ans \n; return 0; }代码量不多但核心就一句话把“至少能被 d 整除”的序列数减去“gcd 恰好是 2d、3d、4d”的序列数剩下的就是“gcd 恰好为 d”的序列数。时间复杂度 O(K log K K log N)空间复杂度 O(K)。我实际测试下来K1e5 的情况下整个过程很轻松瓶颈反而不在复杂度而在取模和边界处理上。提示内层循环我写的是for (int x d d; x K; x d)而不是常见的for (int j 2; j * d K; j)。后一种写法在 K 较大的时候j * d可能会整数溢出换成x d最稳妥。3. 解法二用 φ 函数把答案压缩成一行公式3.1 把容斥展开后发生了什么倒序容斥能跑但还不够优雅。把式子G[d] F[d] - G[2d] - G[3d] - ...继续往下展开你会发现每一个 G[d] 最终都可以表示成若干个 F 的线性组合。真正把这些组合系数总结出来的工具是莫比乌斯函数 μμ(1) 1如果 n 含有平方因子μ(n) 0如果 n 是 k 个互不相同质数的乘积μ(n) (-1)^k有了 μG[d] 的展开式可以写成G[d] Σ_{m1}^{K/d} μ(m) * F[m*d]把答案Ans Σ_{d1}^K d * G[d]代进去然后交换 d 和 t 的求和顺序令 t m*d就能整理出Ans Σ_{t1}^K floor(K/t)^N * ( Σ_{d|t} d * μ(t/d) )右边那个内层求和恰好是一个著名的数论恒等式Σ_{d|n} d * μ(n/d) φ(n)所以答案一下子变成了一个非常漂亮的形式Ans Σ_{t1}^K φ(t) * floor(K/t)^N也就是说这道 Hard 题最后的答案本质上只是对 1 到 K 的 φ 函数做一次加权求和。看到这个结论再回头看整个题目你会有一种“原来如此”的感觉gcd 计数问题最终几乎都落在欧拉函数和莫比乌斯函数这对搭档上。3.2 证明 sum_{d|n} d*μ(n/d)φ(n)这个恒等式看起来很神奇但其实证明不难。设 f(n) Σ_{d|n} d * μ(n/d)。首先 f 是积性函数只要 gcd(a,b)1就有 f(ab)f(a)f(b)因为 μ 是积性函数这个和式也在互质拆分下保持乘积。接下来只需要看素数幂的情形。设 np^e那么 d 只可能是 p^0, p^1, ..., p^e。逐项代入当 dp^e 时n/d1μ(1)1贡献 p^e当 dp^(e-1) 时n/dpμ(p)-1贡献 -p^(e-1)当 d 是更小的 p 的幂时n/d 里含有 p² 因子μ(n/d)0贡献 0所以 f(p^e) p^e - p^(e-1) p^(e-1)(p-1)这正好是 φ(p^e)。由积性扩展到任意正整数 n恒等式就成立了。这个证明没有高深的东西却把“为什么答案里会出现 φ”解释得很清楚。3.3 线性筛 φ 的实现要求 1 到 K 的 φ 函数最无脑的方法是埃氏筛vectorint phi(K 1); for (int i 1; i K; i) phi[i] i; for (int i 2; i K; i) { if (phi[i] i) { // i 是质数 for (int j i; j K; j i) { phi[j] - phi[j] / i; } } }原理是φ(n) n * Π(1 - 1/p)对每个质因子 p 做一次phi[j] - phi[j] / p。这个写法简单直观复杂度 O(K log log K)在本题范围内完全够用。如果想追求极致可以用欧拉筛在线性时间内算出所有 φvectorint primes; vectorint phi(K 1); vectorbool notPrime(K 1, false); phi[1] 1; for (int i 2; i K; i) { if (!notPrime[i]) { primes.push_back(i); phi[i] i - 1; } for (int p : primes) { long long v 1LL * i * p; if (v K) break; notPrime[v] true; if (i % p 0) { // 此时 p 是 i 的最小质因子 phi[v] phi[i] * p; break; } else { phi[v] phi[i] * (p - 1); } } }欧拉筛的核心是“每个合数只被它的最小质因子筛掉一次”当i % p 0时立刻 break保证复杂度严格为 O(K)。如果这里不 break同一个合数会被多个质因子重复筛虽然结果可能还是对的但常数和复杂度都会变差。3.4 O(K) 求答案有了 φ 数组最终代码短到令人舒适ll ans 0; for (int t 1; t K; t) { ans (ans phi[t] * qpow(K / t, N)) % MOD; } cout ans \n;复杂度是 O(K K log N)空间 O(K)。相比倒序容斥这个版本不仅少了一层调和级数枚举代码上也简洁很多。从竞赛实战的角度看如果你熟悉两个恒等式和线性筛这几乎是“默写题”。注意这里的qpow(K / t, N)里K / t是整除结果就是 1 到 K 中 t 的倍数的个数和前面 F[t] 的定义完全一致。别把它写成浮点数除法否则精度会出问题。4. 两种方案怎么选复杂度、场景与题目变式4.1 两种解法的对比表对比维度倒序容斥φ 函数公式核心思路用“至少可除”推出“恰好等于”用数论恒等式直接压缩求和时间复杂度O(K log K K log N)O(K K log N)空间复杂度O(K)O(K)代码量中等更短能否输出每个 gcd 的个数分布能G 数组本身就是分布不能只给出总和对数学知识的要求只需要容斥思想需要 φ 函数与恒等式对 K 很大的敏感性1e6 轻松1e7 开始变慢1e6 轻松筛法同样受限于 K如果只是要求总和φ 公式无疑是首选。但如果你在做一个交互式查询系统或者题目要求你输出 gcd 恰好等于每个 d 的序列个数那就必须保留容斥求 G 数组的写法因为 φ 公式把 G 的信息已经压缩进加权和里了没法复原出每个 d 的独立数量。4.2 什么时候必须回到容斥我见过不少人学会了 φ 公式之后觉得容斥写法没有必要结果换一道题就翻车了。比如问法变成“对于每个 d 从 1 到 K输出 gcd 恰好等于 d 的序列数量”这时候你用 φ 公式只能得到一个总答案完全没法拆开。再比如题目要求做某种前缀和、或者要对 G[d] 做额外的卷积手里没有 G 数组就很难继续。所以我的建议是两种写法都要能闭着眼写出来。φ 公式是“最优解”倒序容斥是“保底解”也是理解整道题的根基。在真实比赛中先想到倒序容斥再用恒等式去优化是很自然的思维链。4.3 同类题目的延伸从 gcd 到 lcm再到多组询问这道题的思路可以迁移到很多相近的问法上。比如求“所有序列 gcd 恰好等于 1 的个数”那就可以直接用Σ_{d1}^K μ(d) * floor(K/d)^N这就是莫比乌斯反演的标准形式。再比如求 lcm 相关的问题思路会从“每个数都能被 d 整除”变成“每个数都是 d 的因数”计数方式变成 2 的幂次而不是 floor 的 N 次幂但整体框架仍然类似。还有一个很现实的拓展如果有多组询问每次给定 N 和 K但 K 变化范围很大那么可以先预处理 φ 的前缀和然后用整除分块处理floor(K/d)相同的区间。因为floor(K/d)的取值只有 O(√K) 种每组询问可以做到 O(√K log N)比每次都 O(K) 快很多。这个技巧在很多 gcd 求和类问题里都是高频考点。5. 实战踩坑记录与调试建议5.1 最容易翻车的负余数我在第一次提交时自以为取模写得很严谨结果 WA 了一版。当时的错误写法是ll g F[d]; for (int x d d; x K; x d) { g - G[x]; } G[d] g % MOD;问题出在g可能会变成很大的负数。比如连续减去好几个接近 MOD 的数最终g可能是 -4e14这时候直接% MOD得到的还是负数赋给G[d]会导致后面所有的减法都失效。正确的做法是在循环内每减一次都立刻归一化g (g - G[x] MOD) % MOD;因为g和G[x]都在[0, MOD-1]单次减法最小值是负数但一定大于 -MOD加一次 MOD 就能保证结果非负。养成“减法必取模”的习惯后这一整类问题都不会再出这种低级错误。5.2 快速幂与取模的细节std::pow(K / d, N)想都不要想它返回浮点数不取模而且大数下会丢精度。快速幂的指数b用 long long底数a也要 long long因为res * a在 MOD1e97 时大约是 1e18 级别long long 刚好装得下。如果未来遇到更大的 MOD比如 1e12 以上两个数相乘会超过 long long那时候需要换成快速乘或者用__int128兜底。5.3 越界写法的雷区容斥写法里常见的越界源头是for (int j 2; j * d K; j)j * d在 d 和 j 都很大的时候可能先溢出 int 再参与比较导致循环少算甚至死循环。我现在的习惯是直接for (int x d d; x K; x d)安全且语义清楚复制到任何题里都不需要操心溢出问题。5.4 用 N1 和全 K 序列做边界验证有两个特别好的手算边界强烈建议每次写完都测一下。当 N1 时每个序列只有一个元素gcd 就等于这个元素本身所以答案应该是12...K K(K1)/2。用 φ 公式验证Σ φ(d) * floor(K/d)也正好等于这个结果这是一个非常强的检验点。如果两边不等几乎可以确定是 φ 筛写错了或者在容斥里漏了某个倍数。另一个边界是考虑“所有元素都等于 K”的序列。这样的序列只有一个它的 gcd 是 K所以答案至少不小于 K。而在倒序容斥里G[K] F[K] 1^N 1贡献K * 1 K正好对应这个序列。看到答案里包含这个 K说明边界计算没有漏掉 dK 的位置。5.5 对拍评估脚本遇到这种计数题我强烈建议写一个暴力程序对拍。别看 N 和 K 不能太大但在 N≤4、K≤6 的规模下暴力枚举所有序列是完全可行的。比如用 Pythonimport math import itertools def brute(N, K): ans 0 for arr in itertools.product(range(1, K 1), repeatN): g 0 for x in arr: g math.gcd(g, x) ans g return ans然后用它和 C 程序跑同样的小数据逐一对比。我当时的习惯是生成十几个随机小用例保证 N 和 K 都覆盖边界情况比如 N1、K1、N2、K2、N4、K5 等等。这个对拍过程能帮你发现几乎所有逻辑错误比盯着代码看半天高效得多。5.6 小坑汇总表症状原因修法答案变成负数减法后取模时机不对每轮减法立刻(g - G[x] MOD) % MOD大样例跑得很慢用了浮点pow或重复快速幂预处理floor(K/d)^N全部用快速幂边界 K1 时答案不对没有把 G[1] 初始化或漏算手动验算 K1答案应为 1容斥顺序写反从小到大枚举 d必须从 K 到 1 倒序枚举倍数越界j * d溢出 int改成x d的写法φ 筛缺少 φ[1]欧拉筛没有给 phi[1] 赋初值初始化phi[1] 1说起来这些都是我在实际提交里真真切切踩过的坑尤其是负余数那一次WA 了一版才发现是取模时机的问题。如果你时间紧张至少记住两件事倍数枚举用x d减法每轮都取模。把这两个习惯刻进肌肉记忆这类 gcd 计数题的通过率能高一大截。