1. 问题引入从矩阵求和到数论本质最近在复盘一些经典的算法竞赛题目特别是蓝桥杯国赛的压轴题总能发现一些将编程技巧与数学深度结合的好例子。2018年蓝桥杯B组国赛的第六题“矩阵求和”就是这样一个典型。乍一看题目名字你可能会想这不就是一个二维数组累加吗能有多难但当你真正读题并开始思考时才会发现它巧妙地绕开了暴力计算的陷阱将问题引向了数论中一个既优美又强大的工具——欧拉函数。这道题的核心场景是这样的给定一个 n x n 的矩阵矩阵中每个位置 (i, j) 的元素值定义为 gcd(i, j)。也就是说第 i 行第 j 列的元素是 i 和 j 的最大公约数。题目要求计算这个矩阵所有元素之和并对结果取模通常是一个大质数比如 1e97。当 n 的规模达到 10^7 甚至更大时直接二重循环计算 gcd 再求和时间复杂度是 O(n² log n)这显然是无法接受的。这就迫使我们必须去寻找一个不依赖于遍历每个元素的公式。这恰恰是算法竞赛的魅力所在它考察的不是你会不会写循环而是你能不能透过问题的表象看到其背后的数学结构。矩阵中每个元素是 gcd(i, j)而我们需要求所有 gcd(i, j) 的和。这个“所有数对的最大公约数之和”的问题在数论中有一个非常经典的转化思路——利用欧拉函数。理解这个转化不仅是解这道题的关键更是深入理解数论在算法中应用的一个绝佳切入点。接下来我们就一步步拆解这个问题看看如何从暴力思路走到优雅的数论公式并最终实现高效的算法。2. 思路演进从暴力枚举到数学洞察面对“求所有 gcd(i, j) 之和”这个问题最直接的想法就是暴力枚举。我们可以写两层循环i 从 1 到 nj 从 1 到 n计算每一对 (i, j) 的 gcd然后累加。代码写起来非常简单但正如前面所说它的时间复杂度是 O(n² log n)。当 n1000 时大概需要计算百万次 gcd尚可接受但当 n10^5 时计算量就达到了百亿级别完全不可行。所以暴力法只能作为我们验证小规模数据正确性的工具绝非正解。那么优化的方向在哪里我们必须找到一种方法能够不显式地计算每一对 gcd而是通过某种聚合或计数的技巧来快速求和。这里就需要引入一个关键的数学思维计数贡献。我们不再问“每一对的 gcd 是多少”而是反过来问“有多少个数对它们的最大公约数恰好等于 d” 其中 d 是一个从 1 到 n 的可能值。如果我们能对于每一个可能的 d快速求出满足 gcd(i, j) d 的数对 (i, j) 的数量记为 cnt(d)那么最终的答案就可以表示为Sum Σ_{d1}^{n} [ d * cnt(d) ]因为对于每一对最大公约数为 d 的数它对总和的贡献就是 d。现在问题的核心就变成了如何高效计算 cnt(d)。gcd(i, j) d 意味着什么这意味着 d 能同时整除 i 和 j并且 i/d 和 j/d 互质即它们的最大公约数是 1。因为如果 i/d 和 j/d 还有大于 1 的公因子 k那么 i 和 j 的最大公约数就会是 d*k而不是 d 了。因此我们可以进行变量替换令i i / d,j j / d。那么条件gcd(i, j) d就等价于i和j是正整数。1 ≤ i ≤ floor(n/d)1 ≤ j ≤ floor(n/d)gcd(i, j) 1所以cnt(d) 就等于在1 ≤ i, j ≤ floor(n/d)的范围内满足gcd(i, j) 1的数对 (i, j) 的数量。于是一个求 gcd 之和的问题经过两步转化变成了一个求“一定范围内互质数对数量”的问题。而求解互质数对的数量正是欧拉函数大显身手的地方。3. 核心武器欧拉函数与互质数对计数为了计算 cnt(d)我们需要一个高效的方法来计算对于一个给定的上界 m这里 m floor(n/d)有多少对 (i, j) 满足 1 ≤ i, j ≤ m 且 gcd(i, j) 1。这里有一个非常巧妙且重要的公式。我们不妨先固定 i考虑有多少个 j (1 ≤ j ≤ i) 满足 gcd(i, j) 1。根据定义这正好就是欧拉函数 φ(i) 的值——它表示小于等于 i 的正整数中与 i 互质的数的个数。但是我们需要的是所有 i 和 j 的组合i 和 j 是对称的。如果我们简单地将所有 φ(i) 相加i 从 1 到 m得到的是满足j ≤ i且gcd(i, j)1的数对数量。为了得到全部数对我们可以利用对称性。考虑所有满足1 ≤ i, j ≤ m且gcd(i, j) 1的数对。它们可以分为三类i ji ji j由于对称性情况1和情况3的数量是相等的。情况2即 i j 且 gcd(i, i)1这意味着 i 必须等于 1因为只有 gcd(1,1)1当 i1 时gcd(i,i)i1。所以情况2只有 (1,1) 这一对。那么如何计算情况1或情况3的数量呢我们可以遍历 i 从 1 到 m对于每个 i与它互质且比它大的 j 有多少个这似乎不好直接算。更通用的方法是计算所有互质对 (i, j)其中 i 和 j 没有大小关系限制。这里直接给出一个经典结论在 1 到 m 的范围内互质数对 (i, j) 的总数包括 ij1 的情况等于1 2 * Σ_{k1}^{m} φ(k)。这个公式的推导可以这样理解我们先计算所有满足1 ≤ i ≤ m, 1 ≤ j ≤ i且gcd(i, j)1的数对。对于每个 i这样的 j 有 φ(i) 个。所以总数是Σ_{i1}^{m} φ(i)。这个集合包含了所有 i j 的互质对因为 j ≤ i。所有 i j 的互质对即 (1,1)。它缺少的是 i j 的互质对。由于对称性i j 的互质对数量等于 i j 的互质对数量。而 i j 的互质对数量是多少呢正是我们上面计算的Σ_{i1}^{m} φ(i)减去那个特殊的 (1,1)因为当 i1, j1 时是相等情况不是大于。所以 i j 的互质对数量是Σ_{i1}^{m} φ(i) - 1。因此总的互质对数量 (i j 的对) (i j 的对) (i j 的对) [Σφ(i) - 1] 1 [Σφ(i) - 1] 2 * Σφ(i) - 1。这和我们上面的公式1 2 * Σφ(i)是等价的吗注意上面的公式是从1开始求和即Σ_{k1}^{m} φ(k)。我们推导出的2 * Σφ(i) - 1中的Σφ(i)也是从1到m。所以两个公式是一致的1 2 * Σ_{k1}^{m} φ(k) 2 * Σ_{k1}^{m} φ(k) 1这里似乎有个笔误。让我们重新严谨推导一遍。设S(m) Σ_{i1}^{m} φ(i)。 考虑有序对 (i, j)其中 1 ≤ i, j ≤ m。 我们想计算满足 gcd(i, j) 1 的有序对数量。方法一更清晰固定 i计算有多少个 j 满足 gcd(i, j)1。对于每个固定的 ij 可以从 1 取到 m其中与 i 互质的 j 的数量并不是简单的 φ(i)因为 φ(i) 定义是 j ≤ i 且互质。当 j i 时也可能与 i 互质。所以不能直接用 φ(i)。我们需要一个更基础的公式。事实上有一个数论中常用的恒等式Σ_{d|n} φ(d) n这个公式的意思是对于任意正整数 n它的所有正因子 d 的欧拉函数值之和等于 n 本身。利用这个公式我们可以解决我们的计数问题。我们要求的是Σ_{i1}^{m} Σ_{j1}^{m} [gcd(i, j) 1]其中[ ]是艾弗森括号条件为真时值为1否则为0。这里有一个技巧利用上面提到的恒等式。因为Σ_{d|gcd(i,j)} φ(d) gcd(i, j)。特别地当 gcd(i, j) 1 时右边为1。而左边当 gcd(i,j)1 时d 能整除 gcd(i,j) 意味着 d 必须为1。所以Σ_{d|1} φ(d) φ(1) 1。这看起来是平凡的。但反过来我们可以用这个恒等式来“筛选”出 gcd(i,j)1 的数对。实际上更常用的技巧是[gcd(i, j) 1] Σ_{d|gcd(i,j)} μ(d)其中 μ(d) 是莫比乌斯函数。 或者利用欧拉函数恒等式的另一种形式。但针对本题最直接推导出的经典公式是在 1 到 m 的范围内互质有序对 (i, j) 的数量等于Σ_{d1}^{m} φ(d) * floor(m/d)²。这个公式的推导如下我们不是直接数互质对而是数所有对然后按它们的最大公约数分类。 所有数对 (i, j) 的总数是 m²。 这些数对可以根据它们的最大公约数 g 来分类。对于每个 g满足 gcd(i, j) g 的数对其数量等于满足 gcd(i, j) 1 且 1 ≤ i, j ≤ floor(m/g) 的数对数量。我们记这个数量为cnt_pair(floor(m/g))。 所以有m² Σ_{g1}^{m} cnt_pair(floor(m/g))。 根据莫比乌斯反演或二次差分可以解出cnt_pair(m) Σ_{d1}^{m} μ(d) * floor(m/d)²。 而我们知道欧拉函数 φ(n) 和莫比乌斯函数 μ(n) 有一个关系φ(n) Σ_{d|n} μ(d) * (n/d)。利用这个经过一系列变换具体过程略属于数论推导可以得到一个更易于计算且适合本题的表达式cnt_pair(m) 2 * Σ_{i1}^{m} φ(i) - 1。这个公式就是之前提到的。让我们验证一下小数据 m1: φ(1)1。公式计算21 -1 1。互质对只有(1,1)正确。 m2: φ(1)1, φ(2)1。公式计算2(11)-13。互质对有(1,1), (1,2), (2,1)。正确。 m3: φ(1)1, φ(2)1, φ(3)2。公式计算2*(112)-17。互质对有(1,1),(1,2),(1,3),(2,1),(2,3),(3,1),(3,2)。共7个正确。所以我们得到了计算 cnt_pair(m) 的可靠公式cnt_pair(m) 2 * Σ_{i1}^{m} φ(i) - 1。因此回到我们的原始问题cnt(d)即满足1 ≤ i, j ≤ n且gcd(i, j) d的数对数量就等于cnt_pair( floor(n/d) )。 即cnt(d) 2 * Σ_{k1}^{floor(n/d)} φ(k) - 1。注意这里有一个边界细节。当 d n 时floor(n/d)0求和为空公式结果为 -1实际上当 m0 时我们认为互质对数量为0。所以在应用公式时如果 floor(n/d) 0则直接令 cnt(d)0。在后续计算中d 的循环范围是 1 到 n当 d 较大时floor(n/d) 会变小公式依然适用但需要注意 m0 或 m1 时的边界情况。公式2*S(m)-1在 m0 时不适用结果为 -1需要单独处理为0。4. 算法实现欧拉筛法与前缀和优化有了公式Sum Σ_{d1}^{n} [ d * cnt(d) ]和cnt(d) cnt_pair( floor(n/d) )以及cnt_pair(m) 2 * S(m) - 1 (当 m1)其中S(m) Σ_{k1}^{m} φ(k)我们的算法框架就清晰了预处理出 1 到 n 的所有欧拉函数值 φ(i)并同时计算出它们的前缀和 S(i)。遍历 d 从 1 到 n计算m floor(n/d)。如果 m 1则pair_cnt 2 * S(m) - 1否则pair_cnt 0。将d * pair_cnt累加到最终答案中。由于答案可能很大需要在每次累加后进行取模操作。这里最大的挑战在于第一步当 n 很大比如 10^7时如何快速求出 1 到 n 所有数的欧拉函数值如果对每个数单独用公式φ(n) n * Π_{p|n} (1 - 1/p)计算时间复杂度约为 O(n√n)对于 10^7 仍然太慢。这就需要用到算法竞赛中一个非常经典且重要的技巧欧拉筛法线性筛求欧拉函数。它可以在 O(n) 的时间复杂度内一次性求出 1 到 n 所有数的欧拉函数值。其原理基于欧拉函数的几个性质如果 p 是质数则 φ(p) p - 1。如果 p 是质数且 p 能整除 n那么 φ(n*p) φ(n) * p。如果 p 是质数且 p 不能整除 n那么 φ(n*p) φ(n) * (p - 1)。欧拉筛法在筛选质数的过程中可以同时根据这些性质递推求出每个合数的 φ 值。下面给出具体的代码实现逻辑和注释。首先我们需要初始化几个数组phi[i]: 存储数字 i 的欧拉函数值。prime[]: 存储筛选出来的质数。is_prime[i]或vis[i]: 标记 i 是否为质数或是否被访问过。const int MAXN 1e7 10; // 根据题目n的最大范围设定 int phi[MAXN]; // 欧拉函数值 int prime[MAXN], cnt; // 质数表 bool vis[MAXN]; // 标记是否被筛掉true表示非质数 long long sum_phi[MAXN]; // 欧拉函数前缀和 void get_eulers(int n) { phi[1] 1; // 定义 for (int i 2; i n; i) { if (!vis[i]) { // i是质数 prime[cnt] i; phi[i] i - 1; // 质数的欧拉函数值为 i-1 } // 遍历当前已找到的所有质数 for (int j 0; j cnt i * prime[j] n; j) { vis[i * prime[j]] true; // 标记合数 if (i % prime[j] 0) { // 关键prime[j] 是 i 的最小质因子 phi[i * prime[j]] phi[i] * prime[j]; break; // 保证每个数只被其最小质因子筛一次 } else { phi[i * prime[j]] phi[i] * (prime[j] - 1); } } } // 计算前缀和 sum_phi[0] 0; for (int i 1; i n; i) { sum_phi[i] sum_phi[i-1] phi[i]; } }这段代码是线性时间复杂度的关键。if (i % prime[j] 0) break;这一行确保了每个合数只会被它的最小质因子筛掉一次从而将复杂度降为 O(n)。在筛的过程中根据 i 和 prime[j] 的关系应用欧拉函数的性质进行递推计算。预处理完 phi 和 sum_phi 之后主计算过程就很简单了const int MOD 1e9 7; long long solve(int n) { get_eulers(n); // 预处理 long long ans 0; for (int d 1; d n; d) { int m n / d; if (m 0) continue; long long pair_cnt (2 * sum_phi[m] - 1) % MOD; ans (ans d % MOD * pair_cnt % MOD) % MOD; } return ans; }4.1 复杂度分析与进一步优化上述算法的时间复杂度分为两部分预处理欧拉函数及前缀和O(n)。主循环遍历 d 从 1 到 nO(n)。总时间复杂度为 O(n)。对于 n10^7 的情况O(n) 的算法在合理的实现下使用整型变量、注意内存布局是可以在规定时间内完成的。内存方面需要存储长度为 n1 的 phi、sum_phi、vis 数组大约需要 (481)*10^7 ≈ 130MB假设 int 4字节long long 8字节bool 1字节这在大多数竞赛环境内存限制通常256MB或512MB中是允许的。然而我们还可以进行一个非常重要的优化将主循环的复杂度从 O(n) 降低到 O(√n)。这个优化基于一个观察floor(n/d)的值在 d 变化时并不是每个值都不同而是会成段地相等。例如n10。 d1, m10 d2, m5 d3, m3 d4, m2 d5, m2 d6, m1 d7, m1 d8, m1 d9, m1 d10, m1可以看到m 的值重复出现了很多次。实际上对于任意正整数 nfloor(n/d)的不同取值大约只有 2√n 个。我们可以通过数论分块也叫除法分块的技巧一次性处理所有使得floor(n/d)相同的连续 d 区间。具体来说对于一个给定的 m floor(n/d)满足 floor(n/d) m 的最大的 d 是 floor(n/m)。所以区间 [d, floor(n / floor(n/d))] 内的所有 d对应的 m 值都相同。优化后的主循环如下long long solve_optimized(int n) { get_eulers(n); long long ans 0; for (int l 1, r; l n; l r 1) { int m n / l; // 当前块内统一的 m 值 r n / m; // 当前块的右边界 // 区间 [l, r] 内的 d其 floor(n/d) 都等于 m if (m 0) continue; long long pair_cnt (2 * sum_phi[m] - 1) % MOD; // 计算 d 在 [l, r] 区间内的和等差数列求和公式 // sum_d (l r) * (r - l 1) / 2 long long cnt r - l 1; // 区间内d的个数 long long sum_d ((l r) % MOD) * (cnt % MOD) % MOD * inv2 % MOD; // inv2是2的模逆元 long long contribution sum_d * pair_cnt % MOD; ans (ans contribution) % MOD; } return ans; }这里用到了等差数列求和公式以及乘法逆元来处理除法取模因为 MOD 是质数 1e972 的逆元 inv2 (MOD1)/2。经过数论分块优化后主循环的迭代次数大约是 2√n 次对于 n10^7迭代次数约 6324 次几乎可以忽略不计。整个算法的时间复杂度主要取决于预处理的 O(n)。实操心得在竞赛中遇到floor(n/i)形式的求和或遍历一定要第一时间想到数论分块优化的可能性。这常常是能否通过大数据测试的关键。同时注意处理乘法取模和除法逆元避免溢出和精度问题。5. 代码实现细节与测试验证将上述所有步骤整合我们得到完整的解决方案。下面给出一个详细的 C 实现包含欧拉筛、前缀和、数论分块以及模运算。#include iostream #include vector using namespace std; const int MAXN 10000000; // 假设 n 最大为 1e7 const int MOD 1000000007; const int INV2 500000004; // 2 在模 MOD 下的逆元 (MOD1)/2 int phi[MAXN 5]; long long sum_phi[MAXN 5]; int prime[MAXN 5]; bool vis[MAXN 5]; int cnt; // 线性筛法求 1~n 的欧拉函数值并计算前缀和 void init_euler(int n) { phi[1] 1; for (int i 2; i n; i) { if (!vis[i]) { prime[cnt] i; phi[i] i - 1; } for (int j 0; j cnt i * prime[j] n; j) { vis[i * prime[j]] true; if (i % prime[j] 0) { phi[i * prime[j]] (long long)phi[i] * prime[j] % MOD; // 实际上phi值不会取模这里先不模以防后续求和出错。更常见的做法是用long long存。 // 更正phi数组存储精确值不取模 phi[i * prime[j]] phi[i] * prime[j]; break; } else { phi[i * prime[j]] phi[i] * (prime[j] - 1); } } } // 计算前缀和前缀和可能需要取模因为后续计算会用到 sum_phi[0] 0; for (int i 1; i n; i) { sum_phi[i] sum_phi[i-1] phi[i]; // 如果担心溢出可以在每一步对MOD取模但注意后续计算 2*sum_phi[m]-1 可能需要精确值。 // 由于 n 最大 1e7 phi[i] 最大接近 1e7前缀和最大约 5e13在 long long 范围内约9e18所以可以暂不取模。 } } long long solve(int n) { init_euler(n); long long ans 0; for (int l 1, r; l n; l r 1) { int m n / l; if (m 0) break; r n / m; // 计算 cnt_pair(m) 2 * S(m) - 1 long long pair_cnt 2 * sum_phi[m] - 1; // 计算 d 在 [l, r] 区间内的和 long long cnt r - l 1; long long sum_d (l r) * cnt / 2; // 等差数列求和结果在 long long 范围内 // 计算贡献并累加注意取模 long long contribution (sum_d % MOD) * (pair_cnt % MOD) % MOD; ans (ans contribution) % MOD; } return ans; } int main() { int n; // 假设输入 n // cin n; n 100; // 示例测试 cout solve(n) endl; return 0; }5.1 测试与验证为了确保代码正确性我们需要用暴力算法对小数据进行验证。暴力算法 O(n² log n) 实现long long brute_force(int n) { long long sum 0; for (int i 1; i n; i) { for (int j 1; j n; j) { sum __gcd(i, j); // 使用内置gcd函数或者自己实现 } } return sum; }我们可以对 n 从 1 到 100或更大只要暴力法能在短时间内算完进行测试比较solve(n)和brute_force(n)的结果是否一致。例如n1: 矩阵只有 [[1]]和为1。n2: 矩阵为 [[1,1], [1,2]]和为 11125。n3: 矩阵为 [[1,1,1], [1,2,1], [1,1,3]]和为 11112111312。运行我们的优化算法应该得到相同的结果。这是验证数论推导正确性的最直接方法。5.2 边界情况与注意事项n1 的情况公式cnt_pair(m) 2*S(m)-1中当 m1 时S(1)φ(1)1cnt_pair1正确。主循环中d1时m1贡献为 1*11。大数取模在最终累加答案时(sum_d % MOD) * (pair_cnt % MOD) % MOD可以防止中间结果溢出。但需要注意的是在计算pair_cnt 2 * sum_phi[m] - 1时sum_phi[m]可能很大最大约 5e13乘以2仍在 long long 范围内9e18所以直接计算没问题。如果 n 更大比如 10^9我们无法预处理出所有 phi就需要用更高级的方法杜教筛来求 S(m)并且要时刻注意取模。内存与初始化数组phi,sum_phi,vis的大小是 n1当 n1e7 时内存占用在可接受范围。确保全局数组初始化或局部vector初始化为0。时间复杂度预处理 O(n)主循环 O(√n)。对于 n1e7预处理是主要开销在一般的 OJ 上大约需要 0.2-0.5 秒取决于实现优化可以接受。踩坑记录在最初实现时最容易出错的地方是在欧拉筛中计算 phi 的递推式特别是if (i % prime[j] 0)分支中的phi[i * prime[j]] phi[i] * prime[j]。一定要理解其推导因为 prime[j] 是 i 的质因子所以 i * prime[j] 和 i 的质因子集合相同只是 prime[j] 的指数增加了1。根据欧拉函数公式 φ(n) n * Π (1 - 1/p)n 变成了原来的 prime[j] 倍而连乘积部分不变所以 φ 也变为原来的 prime[j] 倍。如果记错了这个性质结果就会出错。务必用小数据测试验证。6. 总结与举一反三回顾整个解题过程我们从一道看似是编程题的“矩阵求和”出发深入到了数论的核心领域。关键步骤在于问题转化将求所有gcd(i,j)的和转化为求每个公约数 d 的贡献次数cnt(d)进而转化为求一定范围内互质数对的个数cnt_pair(m)最后利用欧拉函数的前缀和来高效计算。这体现了算法竞赛中“化数为形”和“贡献法”的典型思想。欧拉函数 φ(n) 在此扮演了核心角色。它不仅是数论中的一个基础函数更是连接“最大公约数”与“互质”概念的桥梁。通过线性筛法 O(n) 预处理欧拉函数我们解决了算法的效率瓶颈。而数论分块的优化则是处理floor(n/d)形式求和的利器将复杂度从 O(n) 降为 O(√n)这对于 n 很大的情况至关重要。这道题的价值不仅仅在于解决一个具体问题更在于提供了一套解决类似“gcd求和”问题的模板。例如如果题目变形为求Σ_{i1}^{n} Σ_{j1}^{n} gcd(i, j)^kk次方和或者求Σ_{i1}^{n} Σ_{j1}^{m} gcd(i, j)矩形矩阵思路是完全一致的枚举公约数 d计算cnt(d)然后求和。区别在于cnt_pair(m)的计算公式和最终 d 的贡献形式d 还是 d^k。对于矩形矩阵n≠mcnt_pair(n, m)表示 1≤i≤n, 1≤j≤m 的互质对数量的公式会稍有不同但依然可以用欧拉函数前缀和结合二维数论分块来解决。其公式为cnt_pair(n, m) Σ_{i1}^{min(n,m)} φ(i) * floor(n/i) * floor(m/i)这可以通过莫比乌斯反演得到。在实际编码中除了掌握公式更要注重细节线性筛的写法、数组大小的设置、long long 防溢出、取模运算、数论分块边界的处理。多写多练才能将这些知识内化在比赛中快速识别并应用。最后个人在多次实现此类问题后的一点体会是数论题的代码往往不长但思维密度极高。在比赛时如果遇到类似“求和”、“计数”且与 gcd 相关的问题应优先考虑是否可以通过枚举公约数并利用欧拉函数、莫比乌斯函数等数论函数来化简。预处理这些函数的前缀和往往是优化到 O(n) 或 O(√n) 的关键。把这道“矩阵求和”吃透相当于掌握了解决一类问题的通用钥匙。