GCD求和难题解析：容斥原理与欧拉函数反演

“E - Sum of gcd of Tuples (Hard)”，这题名我第一次看到时心里就咯噔一下：Hard标签加上Sum of gcd，八成又是一道数论计数题。翻译成大白话就是：给定两个正整数 N 和 K，把所有长度为 N、每个元素都在 1 到 K 之间的序列全部枚举出来，对每个序列求一次最大公约数，最后把所有 gcd 加起来，输出对 1e9+7 取模的结果。

题目数据范围是 1 ≤ N, K ≤ 1e5。第一次写暴力程序试跑时，我还没意识到问题的严重性：序列总数是 K 的 N 次方，别说 1e5，哪怕 N=10、K=10 也是 100 亿个序列，这数字直接劝退。所以这道题硬做必死，必须靠数论工具把“枚举所有序列”这个维度压缩掉。这篇文章适合三类人：准备算法竞赛的选手、在笔试面试前想补数论计数的朋友、以及单纯好奇“gcd 求和”这类问题怎么优雅求解的读者。我会把两条主流做法——倒序容斥和欧拉函数反演——完整拆开，附上能直接跑的代码和踩坑记录。

1. 题目在问什么：N 维数组的全序列 gcd 求和

1.1 把题意翻译成人话

假设 N=2、K=2，那所有序列一共 4 个：(1,1)的 gcd 是 1，(1,2)的 gcd 是 1，(2,1)的 gcd 是 1，(2,2)的 gcd 是 2，总和就是 5。题目就是让你计算这个总和对 1e9+7 取模。

这个定义本身没有什么高深的地方，小学生都能理解“最大公约数”是啥。真正麻烦的是“所有序列”四个字：长度为 N，每个位置有 K 种选择，序列总数是 K^N。N 和 K 都是 1e5，K^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) 是可接受的——K=1e5 时，调和级数展开也就约 120 万次操作，完全在比赛时限内。如果 K 变成了 1e7 甚至 1e9，做法又要换血：欧拉函数配合数论分块可以做到 O(√K log N)，这个变化我放到第 5 节讲。先把基础版本吃透，再谈优化。

1.3 先给暴力判死刑

有人可能会说：N 和 K 最大 1e5，但题目又没让输出方案，只是求和，我能不能动态规划？比如设 dp[i][g] 表示长度为 i 的序列 gcd 为 g 的个数？到这一步已经接近正解了，但转移时你还是避不开枚举下一个数，复杂度至少 O(N·K²)。K=1e5 时没法做。

还有人想用容斥预处理每一个 gcd 的倍数个数，这个思路就是正解之一。但无论哪条路，都要先承认一个前提：不能枚举序列。我写了个暴力程序对拍小数据时，N=4、K=6 就开始肉眼可见地卡顿，N=5、K=8 基本跑不动。暴力只能用来验证小数据，不能用来求解。想清楚这一点，再看下面的解法就顺了。

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，可能更大。比如 d=2 时，序列(4, 4)也满足“每个元素都是 2 的倍数”，但它真正的 gcd 是 4。更准确地说，all[d]包含了所有 gcd 为 d、2d、3d、…… 的序列。所以它是“gcd 能被 d 整除”的序列数，不是“gcd 恰好等于 d”的序列数。

如果我能得到一组数组cnt[d]，表示“gcd 恰好等于 d”的序列数，那答案就是 Σ d·cnt[d]。现在问题变成了：怎么从all反推cnt？

2.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 答案求和与复杂度分析

完整流程分三步：

1. 预处理all[d] = powmod(K/d, N, MOD)。
2. 从 d=K 到 d=1 倒序遍历，对每个 d 枚举倍数 j=2d,3d,...，执行cnt[d] = (cnt[d] - cnt[j] + MOD) % MOD。
3. 累加答案ans = (ans + d * cnt[d]) % MOD。

复杂度方面，第二层枚举的循环总次数是 ⌊K/1⌋ + ⌊K/2⌋ + ... + ⌊K/K⌋ ≈ K ln K，在 K=1e5 时约为 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。于是：

答案 = Σ_{d=1}^K φ(d) · (⌊K/d⌋)^N

这个式子漂亮多了：没有数组互相减，没有倒序遍历，只需要把每个 d 的欧拉函数值和快速幂结果乘起来累加。我第一次推出这个结论时拍了一下大腿——原来 Hard 题的真正难点就在“交换求和次序”这一步。

3.3 两种做法的联系与选择建议

欧拉函数反演和容斥法本质上是同一个问题的两个侧面。容斥法直接维护“恰好 gcd=d”的数组，欧拉反演则是把 gcd 的数值本身拆成约数贡献的和。从结果看，两者时间复杂度都在 O(K log N) 附近，K=1e5 时都能跑进几十毫秒。

我个人的选择策略是：

• 如果现场推公式状态好，优先用欧拉反演，代码短，不容易在遍历方向上翻车。
• 如果一时间忘了 φ 的预处理细节，容斥法更稳，因为它只需要快速幂和倍数枚举，不需要欧拉函数的知识。
• 如果题目里的函数不是 f(g)=g，而是 f(g)=g²、f(g)=1（求 gcd=1 的序列数）等，容斥法可以改维护加权和，欧拉法则看 f 的约数展开形式——通常我会先判断 f 能不能写成约数函数的卷积，能写就反演，不能写就容斥硬扛。

4. 可直接运行的代码与实测数据

4.1 容斥版实现（C++17）

#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; vector<ll> 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]，计算完后它变成“恰好 gcd=d 的序列数”。因为 d 是倒序遍历的，访问cnt[j]时 j>d，此刻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; vector<bool> is_composite(K + 1, false); vector<int> primes; vector<ll> 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; vector<int> 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 进制数加一，从末位开始进位，比递归写起来更不容易爆栈。我对拍时遇到过一个很有意思的情况：两种高级算法在小数据下都能对上，但调到 K=20、N=3 时容斥版突然不对。后来发现是我的对拍程序本身把序列枚举漏了一种情况——末位进位到尽头时直接 break 的逻辑写错了。教训是：对拍前先单独测暴力函数本身，比如 N=1、K=5 时暴力结果必须是 15，这样能确认参照物可信。

4.4 性能实录与边界测试

在我的机器上，K=1e5、N=1e5 时，容斥版耗时约 30ms，欧拉反演版约 50ms，差别很小。瓶颈都在模乘上，不在枚举上。所以选手可以根据自己对哪个方法熟选哪个，不需要纠结性能。

边界测试建议跑这样几组：

这些特殊边界非常值得测试：N=1 可以验证你对 gcd 定义的理解，K=1 可以暴露数组下标越界和空循环问题。如果这些测试都过了，再上大数据基本不会出幺蛾子。

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)) = Σ_{d=1}^K h(d) · (⌊K/d⌋)^N

两个常用特例：

• f(n)=n时，h(d)=φ(d)，就是原题。
• f(n)=[n==1]时，h(d)=μ(d)（莫比乌斯函数），于是“gcd 恰好等于 1 的序列数” =Σ_{d=1}^K μ(d)·(⌊K/d⌋)^N。

至于f(n)=n²这类，理论上可以求出对应的 h，但现场推导往往比容斥慢。我的建议是：优先判断能不能套现成的 φ、μ、1 函数；不行就回到容斥，给每个 d 的贡献乘上 f(d)，思路同样清晰。

5.2 当 K 变大到 1e7 甚至 1e9：数论分块优化

基础做法是 O(K log N)，如果 K=1e7，算下来大概 1.7 亿次模乘，勉强有点悬。K=1e9 时 O(K) 完全不可行。这时要用到数论分块。

观察公式Σ_{d=1}^K φ(d)·(⌊K/d⌋)^N，注意到⌊K/d⌋的取值最多只有 O(√K) 种。对于取值相同的连续区间 [l, r]，贡献可以写成(⌊K/l⌋)^N × Σ_{d=l}^{r} φ(d)。只要预处理 φ 的前缀和S[x] = Σ_{i=1}^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 实战套路清单

几次比赛复盘之后，我把这类题总结成了一张清单：

1. 看到“所有序列的 gcd 求和/求计数”，立刻往约数维度想。
2. 先写出暴力对拍程序，哪怕只能跑 N≤4 的小数据，后面验证公式全靠它。
3. 想用欧拉/莫比乌斯反演，先确认 f(n) 能不能写成约数卷积。
4. 想用容斥，记住倒序遍历和减法取模。
5. 最后扫一眼 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 最大 1e5，cnt 最大 1e9 左右，乘积 1e14 已经超出 int 范围，必须用 long long。更隐蔽的是快速幂里的a * a % mod，a 是 long long 时乘法结果 1e18 恰好卡在 long long 上界附近，不会溢出，但如果你把 mod 换成更大的数，这里就要用__int128过渡。竞赛里养成一个习惯：凡是有两个 1e9 级变量相乘，先转 long long 再说。

6.3 线性筛 phi 的初始化位置

欧拉反演版本里，phi[1]必须手动设为 1。我见过有人只初始化数组为 0，然后从 i=2 开始筛，最后发现φ(1)=0导致答案偏小。同理，is_composite初始化为 false，素数从 i=2 开始识别，别把 0 和 1 标成质数。这些小细节在代码 review 时几乎没人提醒，只有对拍才会把你拉回现实。

6.4 对拍程序要够“脏”

最后一个坑来自我的对拍习惯：不要只测几组精心挑选的小数据，要随机。我的做法是写一个循环，N 随机取 1 到 4，K 随机取 1 到 7，跑 1000 组，全部一致才算放心。随机对拍能覆盖很多手算时根本想不到的边界，比如 N=1、K=6，或者 K 是质数、N 是奇数这种组合。曾经有一次，我的高级算法在 N=3、K=4 这组随机数据上暴露了 bug，而那组数据恰好是手算样例没覆盖的。

我在实际解题时的最后一步是：用大样例实测一下耗时，然后再提交。这道题的两个做法我前前后后写了不下五遍，每一遍都会在某些细节上栽跟头。如果你正要备战比赛，我的建议很简单：先别看题解，自己把暴力对拍程序写出来，然后尝试从“统计 gcd 恰好为 d 的序列数”这个角度推一遍容斥，再推一遍欧拉反演。两道都推通之后，你收获的不只是这道题的解法，而是整个“约数反演”类题目的通用思考方式，以后遇到再诡异的 gcd 计数题，心里都会有个底。