题目链接:哆啦A梦传送门
题意:
求 ∑ i = 1 n ∑ j = 1 n l c m ( i , j ) \begin{aligned}\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j) \end{aligned} i=1∑nj=1∑nlcm(i,j)
我们设:
a
n
s
=
∑
i
=
1
n
∑
j
=
1
n
l
c
m
(
i
,
j
)
=
∑
d
=
1
n
d
∑
i
=
1
n
d
∑
j
=
1
n
d
[
(
i
,
j
)
=
1
]
i
∗
j
\begin{aligned}ans&=\sum_{i=1}^{n}\sum_{j=1}^{n}lcm(i,j)\\ &=\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}[(i,j)=1]i*j \end{aligned}
ans=i=1∑nj=1∑nlcm(i,j)=d=1∑ndi=1∑dnj=1∑dn[(i,j)=1]i∗j
这步公式推导见:详情
参考博客
我们知道:
∑ i = 1 n [ ( n , i ) = 1 ] ∗ i = n ∗ φ ( n ) + [ n = = 1 ] 2 \begin{aligned}\sum_{i=1}^{n}[(n,i)=1]*i=\frac{n*\varphi(n)+[n==1]}{2}\end{aligned} i=1∑n[(n,i)=1]∗i=2n∗φ(n)+[n==1]
那么上式就可化为:
∑
d
=
1
n
d
∑
i
=
1
n
d
i
∑
j
=
1
n
d
[
(
i
,
j
)
=
1
]
j
=
∑
d
=
1
n
d
(
(
2
∑
i
n
d
i
∑
j
=
1
i
[
(
j
,
i
)
=
1
]
∗
j
)
−
1
)
.
.
.
.
.
.
解
释
1
=
∑
d
=
1
n
d
(
(
2
∑
i
n
d
i
i
∗
φ
(
i
)
+
[
i
=
=
1
]
2
)
−
1
)
=
∑
d
=
1
n
d
∑
i
n
d
i
2
φ
(
i
)
\begin{aligned}&\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}i\sum_{j=1}^{\frac{n}{d}}[(i,j)=1]j\\ =&\sum_{d=1}^{n}d((2\sum_{i}^{\frac{n}{d}}i\sum_{j=1}^{i}[(j,i)=1]*j)-1)......解释1\\ =&\sum_{d=1}^{n}d((2\sum_{i}^{\frac{n}{d}}i\frac{i*\varphi(i)+[i==1]}{2})-1)\\ =&\sum_{d=1}^{n}d\sum_{i}^{\frac{n}{d}}i^2\varphi(i)\end{aligned}
===d=1∑ndi=1∑dnij=1∑dn[(i,j)=1]jd=1∑nd((2i∑dnij=1∑i[(j,i)=1]∗j)−1)......解释1d=1∑nd((2i∑dni2i∗φ(i)+[i==1])−1)d=1∑ndi∑dni2φ(i)
解释1:因为两倍的话,我们多算了(1,1),故要减去1。两倍的时候,其实(x,x)是算了两次,但因为它们不互质,其实最后也没算进去的。
此时我们只需分块就能算出上式的结果,但n太大,没法预处理
i
2
ϕ
(
i
)
i^2\phi(i)
i2ϕ(i)的前n项和。
故我们知道杜教筛就是求这个的了,我们只需把
i
2
ϕ
(
i
)
i^2\phi(i)
i2ϕ(i)的前n项和给筛出来,题目就迎刃而解了。
我们设
h
(
x
)
=
x
2
φ
(
x
)
h(x)=x^2\varphi(x)
h(x)=x2φ(x),
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
S(n)=\sum_{i=1}^{n}h(i)
S(n)=∑i=1nh(i)
我们知道
n
=
∑
d
∣
n
φ
(
d
)
\begin{aligned}n=\sum_{d|n}\varphi(d)\end{aligned}
n=d∣n∑φ(d)
跟上面式子联想一下,有个
x
3
x^3
x3的关系在里面。
∑
d
=
1
n
d
3
=
n
2
∗
(
n
+
1
)
2
4
\begin{aligned}\sum_{d=1}^{n}d^3=\frac{n^2*(n+1)^2}{4}\end{aligned}
d=1∑nd3=4n2∗(n+1)2
那么我们就有文章可做了。
∑
d
=
1
n
d
3
=
∑
d
=
1
n
∑
i
∣
d
φ
(
i
)
∗
d
2
=
∑
d
=
1
n
∑
i
∣
d
h
(
i
)
(
d
i
)
2
=
∑
d
=
1
n
∑
i
=
1
n
d
h
(
i
)
(
d
∗
i
i
)
2
.
.
.
.
.
.
.
解
释
2
=
∑
d
=
1
n
d
2
∑
i
=
1
n
d
h
(
i
)
\begin{aligned}\sum_{d=1}^{n}d^3&=\sum_{d=1}^{n}\sum_{i|d}\varphi(i)*d^2\\ &=\sum_{d=1}^{n}\sum_{i|d}h(i)(\frac{d}{i})^2\\ &=\sum_{d=1}^{n}\sum_{i=1}^{\frac{n}{d}}h(i)(\frac{d*i}{i})^2.......解释2\\ &=\sum_{d=1}^{n}d^2\sum_{i=1}^{\frac{n}{d}}h(i) \end{aligned}
d=1∑nd3=d=1∑ni∣d∑φ(i)∗d2=d=1∑ni∣d∑h(i)(id)2=d=1∑ni=1∑dnh(i)(id∗i)2.......解释2=d=1∑nd2i=1∑dnh(i)
解释2:我们之前是枚举每个数d的因子,也就是h(i)中i表示的是d的因子,
d
i
\frac{d}{i}
id表示这个数d的另一个因子。
那么我们就可以等价于:枚举因子t为1到n,那么说明有
n
t
\frac{n}{t}
tn个数a,它们的因子都为t,即根据之前的表达式,h(i)中i还是为因子,此时的另外一个因子就为
i
∗
t
i
\frac{i*t}{i}
ii∗t 。也就是枚举每个数,它对它倍数的贡献。
那么上式就很明显了:
n 2 ∗ ( n + 1 ) 2 4 = ∑ i = 1 n h ( i ) + ∑ d = 2 n d 2 ∑ i = 1 n d h ( i ) 恒 等 于 S ( n ) = n 2 ∗ ( n + 1 ) 2 4 − ∑ d = 2 n d 2 S ( n d ) \begin{aligned}&\frac{n^2*(n+1)^2}{4}=\sum_{i=1}^{n}h(i)+\sum_{d=2}^{n}d^2\sum_{i=1}^{\frac{n}{d}}h(i)\\ 恒等于&S(n)=\frac{n^2*(n+1)^2}{4}-\sum_{d=2}^{n}d^2S(\frac{n}{d})\end{aligned} 恒等于4n2∗(n+1)2=i=1∑nh(i)+d=2∑nd2i=1∑dnh(i)S(n)=4n2∗(n+1)2−d=2∑nd2S(dn)
故最终结果就为:
a n s = ∑ d = 1 n d ∗ S ( n d ) \begin{aligned}ans=\sum_{d=1}^{n}d*S(\frac{n}{d})\end{aligned} ans=d=1∑nd∗S(dn)
我们只需把S()的前6000000预处理出来,再杜教筛一下就求出来了。
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<tr1/unordered_map>
using namespace std;
typedef long long LL;
const LL mod=1000000007;
const int N=6e6;
LL sum[N];///存放S的前N项
bool vis[N+10];
LL phi[N+10];
int pri[N+10],tot;
tr1::unordered_map<LL,LL> w;
LL inv2,inv4,inv6;
void init()
{
phi[1]=1;
for(int i=2;i<=N;i++)
{
if(!vis[i]){
phi[i]=i-1;
pri[++tot]=i;
}
for(int j=1;j<=tot&&i*pri[j]<=N;j++)
{
int t=pri[j];
vis[i*t]=1;
if(i%t==0){
phi[i*t]=phi[i]*t;break;
}
phi[i*t]=phi[i]*(t-1);
}
}
///预处理 \sum_{i=1}^{N}i*i*phi[i]
for(int i=1;i<=N;i++)
{
sum[i]=(sum[i-1]+1LL*phi[i]%mod*i%mod*i%mod)%mod;
}
}
LL fast_pow(LL a1,LL n)
{
LL ans=1;
while(n)
{
if(n&1) ans=ans*a1%mod;
a1=a1*a1%mod;
n>>=1;
}
return ans;
}
LL s2(LL x) ///\sum_{i=1}^{x}i*i
{
LL t=x%mod;
return t*(t+1)%mod*(2*t+1)%mod*inv6%mod;
}
LL s1(LL x) ///\sum_{i=1}^{x}i
{
x=x%mod;
return x*(x+1)%mod*inv2%mod;
}
LL djs(LL x)
{
if(x<=N) return sum[x];
if(w[x]) return w[x];
LL t=x%mod;
LL ans=t*t%mod*(t+1)%mod*(t+1)%mod*inv4%mod;
for(LL l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
ans=(ans-(s2(r)-s2(l-1)+mod)%mod*djs(x/l)+mod)%mod;
}
return w[x]=ans;
}
int main()
{
LL n;
inv2=fast_pow((LL)2,mod-2);
inv4=fast_pow((LL)4,mod-2);
inv6=fast_pow((LL)6,mod-2);
init();
scanf("%lld",&n);
LL ans=0;///分块求ans
for(LL l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
ans=(ans+(s1(r)-s1(l-1)+mod)%mod*djs(n/l)%mod)%mod;
}
printf("%lld\n",(ans));
}