题目链接:哆啦A梦传送门
题意:
求: ∑ i = 1 n ∑ j = 1 n g c d ( i , j ) \begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)\end{aligned} i=1∑nj=1∑ngcd(i,j)
题解:
∑ i = 1 n ∑ j = 1 n g c d ( i , j ) = ∑ d = 1 n d ∑ i = 1 n ∑ j = 1 n g c d ( i , j ) = d = ∑ d = 1 n d ∑ i = 1 n d ∑ j = 1 n d g c d ( i , j ) = 1 = ∑ d = 1 n d ( ( 2 ∑ i = 1 n d ∑ j = 1 i [ ( i , j ) = 1 ] ) − 1 ) . . . . . . . 解 释 1 = ∑ d = 1 n d ( ( 2 ∑ i = 1 n d ϕ ( i ) ] ) − 1 ) \begin{aligned}&\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)\\ =&\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)=d\\ =&\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}gcd(i,j)=1\\ =&\sum_{d=1}^{n}d((2\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{i}[(i,j)=1])-1).......解释1\\ =&\sum_{d=1}^{n}d((2\sum_{i=1}^{\frac{n}{d}}\phi(i)])-1)\end{aligned} ====i=1∑nj=1∑ngcd(i,j)d=1∑ndi=1∑nj=1∑ngcd(i,j)=dd=1∑ndi=1∑dnj=1∑dngcd(i,j)=1d=1∑nd((2i=1∑dnj=1∑i[(i,j)=1])−1).......解释1d=1∑nd((2i=1∑dnϕ(i)])−1)
我们再设:
S
(
n
)
=
∑
i
=
1
n
ϕ
(
i
)
S(n)=\sum_{i=1}^{n}\phi(i)
S(n)=∑i=1nϕ(i)
解释1见:
故只需杜教筛下欧拉函数前n项和就可以了。
杜教筛求欧拉函数前n项和
那么最后题目就变为:
S
(
n
)
=
∑
i
=
1
n
−
∑
d
=
2
n
S
(
(
⌊
n
d
⌋
)
S(n)=\sum_{i=1}^{n}-\sum_{d=2}^{n}S((\left \lfloor \frac{n}{d} \right \rfloor)
S(n)=i=1∑n−d=2∑nS((⌊dn⌋)
a
n
s
=
∑
d
=
1
d
d
∗
(
2
∗
S
(
n
d
)
−
1
)
ans=\sum_{d=1}^{d}d*(2*S(\frac{n}{d})-1)
ans=d=1∑dd∗(2∗S(dn)−1)
分块下解决就好了。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<tr1/unordered_map>
using namespace std;
typedef long long LL;
const int N=6e6;
const LL mod=1e9+7;
bool vis[N+10];
int pri[N+10],phi[N+10],tot;
tr1::unordered_map<LL,LL> w;
LL sum[N+10];
LL inv2=500000004;
void init()
{
phi[1]=1;
for(int i=2;i<=N;i++)
{
if(!vis[i]){
pri[++tot]=i;
phi[i]=i-1;
}
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
int x=pri[j];
vis[i*x]=true;
if(i%x==0){
phi[i*x]=phi[i]*x;break;
}
phi[i*x]=phi[i]*(x-1);
}
}
///预处理前N项欧拉函数和
for(int i=1;i<=N;i++)
sum[i]=(sum[i-1]+phi[i])%mod;
}
LL s1(LL x) ///前x项和
{
x=x%mod;
return x*(x+1)%mod*inv2%mod;
}
///杜教筛求欧拉函数前n项和
LL djs(LL x)
{
if(x<=N) return sum[x];
if(w[x]) return w[x];
LL ans=s1(x);
for(LL l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
ans=(ans-(r-l+1)%mod*djs(x/l)%mod+mod)%mod;
}
return w[x]=ans;
}
int main()
{
LL n;
init();
scanf("%lld",&n);
LL ans=0;
for(LL l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+(s1(r)-s1(l-1)+mod)%mod*(2*djs(n/l)%mod-1+mod)%mod)%mod;
}
printf("%lld\n",ans%mod);
return 0;
}