题意
其实就是求
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^ngcd(i,j)
\]
\]
思路
建议先看一下此题的一个弱化版
推一下式子
\[\sum\limits_{i=1}^n\sum\limits_{j=1}^ngcd(i,j)
\]
\]
\[= \sum\limits_{k=1}^nk\sum\limits_{i=1}^n\sum\limits_{j=1}^n[gcd(i,j)=k]
\]
\]
\[=\sum\limits_{k=1}^nk\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{n}{k}}[gcd(i,j)=1]
\]
\]
\[=\sum\limits_{k=1}^nk\sum\limits_{i=1}^{\frac{n}{k}}2\varphi(i)-1
\]
\]
设\(\phi(i)=\varphi(1)+\varphi(2)+...+\varphi(i)\)
则原式
\[=\sum\limits_{i=1}^ni(2\phi(\frac{n}{i})-1)
\]
\]
然后就可以数论分块啦。
至于怎么比较快的求\(\phi(i)\)，基本的杜教筛喽。。
代码
//loj6074
/*
* @Author: wxyww
* @Date: 2019-03-30 12:43:48
* @Last Modified time: 2019-03-30 19:43:10
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7,N = 1000000 + 100,inv2 = (mod + 1) / 2;
ll read() {
ll x=0,f=1;char c=getchar();
while(c<'0'||c>'9') {
if(c=='-') f=-1;
c=getchar();
}
while(c>='0'&&c<='9') {
x=x*10+c-'0';
c=getchar();
}
return x*f;
}
map<ll,ll>ma;
ll n,sum[N];
int dis[N],vis[N],js;
int dls(ll x) {
if(x <= 1000000) return sum[x];
if(ma[x]) return ma[x];
ll ret = 1ll * x % mod * ((x + 1) % mod) % mod * inv2 % mod;
for(ll l = 2,r;l <= x;l = r + 1) {
r = x / (x / l);
ret -= 1ll * (r - l + 1) % mod * dls(x / l) % mod;
ret = (ret + mod) % mod;
}
return ma[x] = ret;
}
void pre() {
sum[1] = 1;vis[1] = 1;
int NN = min(n,1000000ll);
for(int i = 2;i <= NN;++i) {
if(!vis[i]) {
dis[++js] = i;
sum[i] = i - 1;
}
for(int j = 1;j <= js && dis[j] * i <= NN;++j) {
vis[dis[j] * i] = 1;
if(i % dis[j] == 0) {
sum[dis[j] * i] = sum[i] * dis[j] % mod; break;
}
sum[dis[j] * i] = (dis[j] - 1) * sum[i] % mod;
}
sum[i] += sum[i - 1];
sum[i] %= mod;
}
}
signed main() {
n = read();
pre();
ll ans = 0;
for(ll l = 1,r;l <= n;l = r + 1) {
r = n / (n / l);
ans = (ans + (1ll * (r - l + 1) % mod * ((r + l) % mod) % mod * inv2 % mod) % mod * ((2ll * dls(n / l) % mod) - 1 + mod) % mod) % mod;
}
cout<<ans;
return 0;
}