题目分析
瞎推公式
让我们来乱推一波公式……首先我们要求:
∑i=1n∑j=1nijgcd(i,j)
∑
i
=
1
n
∑
j
=
1
n
i
j
g
c
d
(
i
,
j
)
看到gcd,想到莫比乌斯反演,根据我们做莫比乌斯反演题的经验,枚举gcd的值,得到:
∑d=1nd∑i=1n∑j=1nij[gcd(i,j)==d]
∑
d
=
1
n
d
∑
i
=
1
n
∑
j
=
1
n
i
j
[
g
c
d
(
i
,
j
)
==
d
]
然后将 i i 和都提取一个 d d 得到:
根据我们做莫比乌斯反演题的经验,我们这个时候应该用莫比乌斯函数的一个这样的性质:
∑d|nμ(d)=0(n>1)
∑
d
|
n
μ
(
d
)
=
0
(
n
>
1
)
∑d|nμ(d)=1(n=1)
∑
d
|
n
μ
(
d
)
=
1
(
n
=
1
)
所以原式就可变形为:
∑d=1nd3∑i=1⌊nd⌋∑j=1⌊nd⌋ij∑t|gcd(i,j)μ(t)
∑
d
=
1
n
d
3
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
n
d
⌋
i
j
∑
t
|
g
c
d
(
i
,
j
)
μ
(
t
)
先枚举 t t 可以得到:
然后设 sum(x)=∑xi=1i s u m ( x ) = ∑ i = 1 x i ,原式就是:
∑d=1nd3∑t=1⌊nd⌋μ(t)t2sum(⌊nT⌋)2
∑
d
=
1
n
d
3
∑
t
=
1
⌊
n
d
⌋
μ
(
t
)
t
2
s
u
m
(
⌊
n
T
⌋
)
2
设 T=td T = t d 可得:
∑T=1nT2sum(⌊nT⌋)2∑d|Tμ(Td)d
∑
T
=
1
n
T
2
s
u
m
(
⌊
n
T
⌋
)
2
∑
d
|
T
μ
(
T
d
)
d
狄利克雷卷积
嗯,推到这一步,推不下去了,怎么办?
看一波题解看一看后面那个
∑
∑
,是不是让你想到莫比乌斯反演的另一个性质呢?不是
∑d|nμ(d)d=ϕ(n)n
∑
d
|
n
μ
(
d
)
d
=
ϕ
(
n
)
n
所以如果设 id(x)=x i d ( x ) = x ,那么 id∗μ=ϕ i d ∗ μ = ϕ ( ∗ ∗ 指狄利克雷卷积)
所以原始又华丽丽地变身成为:
由于我们知道,形如 ⌊nT⌋ ⌊ n T ⌋ 的东西是可以利用数论分块快速的搞出来的,所以我们的主要矛盾就转化为了
杜教筛
首先把杜教筛的套路摆在桌子上(
S
S
表示的前缀和,
g
g
是一个能让求前缀和变快的迷之函数):
欧拉函数……狄利克雷卷积……这两样东西是不是又让你想起了欧拉函数的一条神秘性质?
∑d|nϕ(d)=n
∑
d
|
n
ϕ
(
d
)
=
n
所以你在思索 g g 到底是什么的时候,会想到一种可能性
然后搞一搞:
(g∗f)(i)=∑d|id2ϕ(d)i2d2=∑d|iϕ(d)i2=i3
(
g
∗
f
)
(
i
)
=
∑
d
|
i
d
2
ϕ
(
d
)
i
2
d
2
=
∑
d
|
i
ϕ
(
d
)
i
2
=
i
3
生活如此美好!再看看我们杜教筛的时候要搞些什么:
S(n)=∑i=1ni3−∑i=2ni2S(ni)
S
(
n
)
=
∑
i
=
1
n
i
3
−
∑
i
=
2
n
i
2
S
(
n
i
)
可是你发现, ∑i3 ∑ i 3 和 ∑i2 ∑ i 2 都不是可以快速求的东西啊……这时候你就要想起一样神奇的玩意儿——
“小学奥数”
或者说,打表找规律,因此发现:
∑i=1ni3=(∑i=1ni)2
∑
i
=
1
n
i
3
=
(
∑
i
=
1
n
i
)
2
这个可以用数学归纳法证明一下,首先对于 n=1 n = 1 的情况,显然成立。
然后若对于某个 n n 成立,那么
然后再打表找规律一波,发现:
∑i=1ni2=n(n+1)(2n+1)6
∑
i
=
1
n
i
2
=
n
(
n
+
1
)
(
2
n
+
1
)
6
对于 n=1 n = 1 的情况,显然成立,若对于某个 n n 成立,那么:
好吧,
于是我们就可以用公式迅速求求求,然后杜教筛一波,在数论分块一波,搞出这道题了。
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=8000000;
int pri[N],is[N+5],tot;
LL mod,n,ans,div2,div6,phi[N+5];
void init() {//预处理phi
phi[1]=1;
for(LL i=2;i<=N;++i) {
if(!is[i]) pri[++tot]=i,phi[i]=i-1;
for(LL j=1;j<=tot&&pri[j]*i<=N;++j) {
is[pri[j]*i]=1;
if(i%pri[j]) phi[pri[j]*i]=1LL*(pri[j]-1)*phi[i]%mod;
else {phi[pri[j]*i]=1LL*pri[j]*phi[i]%mod;break;}
}
}
for(LL i=1;i<=N;++i) phi[i]=(phi[i-1]+1LL*i*i%mod*phi[i]%mod)%mod;
}
//注意x%mod,不这么做会爆炸的
LL sum(LL x) {x%=mod;return (x+1)*x%mod*div2%mod;}
LL sqrsum(LL x) {x%=mod;return x*(x+1)%mod*(x+x+1)%mod*div6%mod;}
LL ksm(LL x,LL y) {
LL re=1;
for(;y;y>>=1,x=x*x%mod) if(y&1) re=re*x%mod;
return re;
}
map<LL,LL> mp;
LL getS(LL x) {//杜教筛求前缀和
if(x<=N) return phi[x];
if(mp[x]) return mp[x];
LL re=sum(x);re=re*re%mod;
for(LL i=2,j;i<=x;i=j+1) {
j=x/(x/i);
re=(re+mod-(sqrsum(j)+mod-sqrsum(i-1))%mod*getS(x/i)%mod)%mod;
}
mp[x]=re;return re;
}
int main()
{
scanf("%lld%lld",&mod,&n);
init();
div2=ksm(2,mod-2),div6=ksm(6,mod-2);
for(LL i=1,j;i<=n;i=j+1) {//数论分块求答案
j=n/(n/i);
LL k1=sum(n/i);k1=k1*k1%mod;
ans=(ans+(getS(j)-getS(i-1)+mod)%mod*k1%mod)%mod;
}
printf("%lld\n",ans);
return 0;
}