首先我们转化。
[egin{aligned}
&sum_{i=1}^{n} sum_{j=i+1}^{n} gcd(i,j) \
&=frac {sum_{i=1}^{n}sum_{j=1}^{n} gcd(i,j) - sum_{i=1}^{n} gcd(i,i)}{2}\
end{aligned}
]
然后就是套路
[sum_{i=1}^{n} sum_{j=1}^{n} gcd(i,j)\
sum_{d=1}^{n} d sum_{i=1}^{frac nd}sum_{j=1}^{frac nd} [gcd(i,j) = 1]\
sum_{d=1}^{n} d sum_{i=1}^{frac nd}sum_{j=1}^{frac nd} sum_{p|i,p|j}mu(p)\
sum_{d=1}^{n} d sum_{p=1}^{frac nd} mu(p)sum_{i=1}^{frac n{dp}}sum_{j=1}^{frac n{dp}} 1\
sum_{d=1}^{n} d sum_{p=1}^{frac nd} mu(p)sum_{i=1}^{frac n{dp}}sum_{j=1}^{frac n{dp}} 1\
]
令
[g(n) = sum_{p=1}^{n} mu(p) sum_{i=1}^{frac np} sum_{j=1}^{frac np} 1\
g(n) = sum_{p=1}^{n} mu(p) frac np frac np\
]
数论分块求
原式变为
[sum_{d=1}^{n} d g(frac nd)\
]
再数论分块。
code
#include<bits/stdc++.h>
using namespace std;
const int N = 2e6;
int prime[N / 10], cnt, p[N + 5], mu[N + 5];
long long g(int n){
long long ret = 0;
for(int l = 1 ,r; l <= n; l = r + 1){
r = min(n, n/(n/l));
ret = ret + 1ll * (mu[r] - mu[l - 1]) * (n/l) * (n/l);
}
return ret;
}
int n;
int main(){
p[0] = p[1] = 1; mu[1] = 1;
for(int i = 2; i <= N; ++ i){
if(!p[i]) { prime[++ cnt] = i; mu[i] = -1; }
for(int j = 1; j <= cnt && 1ll * prime[j] * i <= N; ++ j){
p[prime[j] * i] = 1;
if(i % prime[j] == 0) break;
mu[prime[j] * i] = -mu[i];
}
}
for(int i = 1; i <= N; ++ i) mu[i] += mu[i - 1];
scanf("%d",&n);
long long ans = 0;
for(int l = 1, r; l <= n; l = r + 1){
r = min(n, n/(n/l));
ans = ans + 1ll * (r + l) * (r - l + 1) / 2 * g(n/l);
}
ans -= 1ll * (1 + n) * n / 2;
printf("%lld
",ans / 2);
return 0;
}