題目描述

給定 \(n\) ，求：

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

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

前置知識：

推式子的能力，尤拉函式，字首和

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]; // 按照上面的公式進行計算
	for (int i = 1;i <= n; i++) //求字首和
		s[i] = s[i - 1] + f[i];
	printf("%lld\n" ,s[n]);
	return 0;
}