题解 P1390 公约数的和

题目描述

给定 \(n\) ，求：

\[\sum_{i=1}^n\sum_{j=i+1}^n \gcd(i,j) \]

\(1 \leq n\leq 2\times 10^6\)

顺便一提，本题有多倍经验： UVA11426 拿行李（极限版） GCD - Extreme (II) ， UVA11424 GCD - Extreme (I) ， SP3871 GCDEX - GCD Extreme 。

前置知识：

推式子的能力，欧拉函数，前缀和

Solution

为了方便，我们设：

\[f(n)=\sum_{i=1}^n\sum_{j=i+1}^n \gcd(i,j) \]

同时，用 \((i,j)\) 表示 \(\gcd(i,j)\) 。

我们来举几个简单的例子看看有什么规律：

\[f(1) = 0 \\ f(2) = (1,2) \\ f(3) = (1,2) + (1,3) + (2,3) \\ f(4) = (1,2) + (1,3) + (1,4) + (2,3) + (2,4) + (3,4) \]

把同样的项对其看看：

\[f(1) = 0 \\ f(2) = (1,2) \\ f(3) = (1,2) + (1,3) + (2,3) \\ f(4) = (1,2) + (1,3) + (2,3) + (1,4) + (2,4) + (3,4) \]

我们可以发现：

\[f(n)=f(n-1)+\sum_{i=1}^{n-1}\gcd(i,n) \]

证明一下：

\[f(n)=\sum_{i=1}^n\sum_{j=i+1}^n \gcd(i,j) \\ f(n-1)=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n-1} \gcd(i,j) \]

当 \(n-1\) 变成 \(n\) 的时候， \(i\) 可以循环到 \(n\) ，但是 \(j\) 却不能继续循环，所以当 \(i=n\) 时，对答案无贡献。

那么对于 \(\sum_{i=1}^{n-1}\) 里面的每一项，相当于都加上了 \(\gcd(i,n)\) （这里之比 \(f(n-1)\) 多了一个 \(j=n\) 的情况）。

于是我们就得到了：

\[f(n) = (\sum_{i=1}^{n-1}\sum_{j=i+1}^{n-1} \gcd(i,j)) + \sum_{i=1}^{n-1}\gcd(i,n) \\ \therefore f(n) = f(n-1)+\sum_{i=1}^{n-1}\gcd(i,n) \]

我们可以令 \(g(n)=\sum_{i=1}^n\gcd(i,n)\) ，这样的话：

\[f(n) = g(1)+g(2)+g(3)+\cdots+g(n-1)\\ =\sum_{i=1}^{n-1} g(i) \]

也就是对 \(g\) 做一次前缀和。

那么我们考虑如何求出 \(g(n)\) 。

方法一：枚举约数：

我们可以尝试推式子：

\[g(n)=\sum_{i=1}^{n} \gcd(i,n) \\ = \sum_{k|n} k \times \sum_{i=1}^{n} [\gcd(i,n) = k] \\ = \sum_{k|n} k \times \sum_{i=1}^{n} [\gcd(\frac{i}{k},\frac{n}{k})=1] \\ =\sum_{k|n} k \times\sum_{i=1}^{\frac{n}{k}} \gcd(i,\frac{n}{k}=1) \\ = \sum_{k|n} k \times \varphi(\frac{n}{k}) \]

注意：因为这里的 \(g(n)\) 实际计算时不需要计算 \(\gcd(n,n)\) ，所以特殊定义 \(\varphi(1)=0\) 、

这样，先预处理出所有的 \(\varphi\) ，就可以枚举 \(n\) 的所有约数 \(k\) ，就可以求出 \(g(n)\) 了。

我们需要求出所有 \(i\in[1,n]\) 之间所有 \(g(i)\) 的值，然后求一次前缀和。

每次求 \(g(i)\) 的时间复杂度是 \(\mathcal{O}(\sqrt{n})\) ，总的时间复杂度 \(\mathcal{O}(n\sqrt{n})\) ，但是 \(n\) 是 \(4\times 10^6\) ，会超时，代码就不贴了。

方法二：倍数法

看到上面的枚举约数，又看到要求出所有的 \(g\) 值，是不是跟我们判断素数的情况很像？

那么我么也可以考虑像埃氏筛法那样做，循环一次 \(i\) ，然后枚举 \(i\) 的所有倍数，累加答案，具体见代码。

时间复杂度和埃氏筛是一样的： \(\mathcal{O}(n\log\log n)\) 。

代码如下：

#include <cstdio>
#include <cstring>
#define int long long
const int N = 2e6 + 5;
int prime[N] ,cnt; bool isprime[N];
int phi[N];
inline void Prime(int n) { //线性筛求 phi ,原理就不解释了，注意这里定义 phi[1] = 0
	memset(isprime ,true ,sizeof(isprime));
	isprime[0] = isprime[1] = false;
	for (int i = 2;i <= n; i++) {
		if (isprime[i]) {
			prime[++cnt] = i;
			phi[i] = i - 1;
		}
		for (int j = 1;j <= cnt && prime[j] * i <= n; j++) {
			isprime[prime[j] * i] = false;
			if (i % prime[j] == 0) {
				phi[i * prime[j]] = phi[i] * prime[j];
				break;
			}
			phi[i * prime[j]] = phi[i] * (prime[j] - 1);
		}
	}
}
int n;
int f[N] ,s[N];
signed main() {
	scanf("%lld" ,&n);
	Prime(n);
	for (int i = 1;i <= n; i++) //注意 i 从 1 开始循环
		for (int j = i * 2;j <= n; j += i) // 枚举 i 的所有倍数 j
			f[j] += i * phi[j / i]; // 按照上面的公式进行计算
   	//因为对于 f[i] 并没有累加 i * phi[1] 所以这里的 f[i] 代表的是上面的 g[i - 1]
	for (int i = 1;i <= n; i++) //求前缀和
		s[i] = s[i - 1] + f[i];
	printf("%lld\n" ,s[n]);
	return 0;
}