题目分析
瞎推公式
让我们来乱推一波公式……首先我们要求:
∑i=1n∑j=1nijgcd(i,j)∑i=1n∑j=1nijgcd(i,j)
看到gcd,想到莫比乌斯反演,根据我们做莫比乌斯反演题的经验,枚举gcd的值,得到:
∑d=1nd∑i=1n∑j=1nij[gcd(i,j)==d]∑d=1nd∑i=1n∑j=1nij[gcd(i,j)==d]
然后将ii和都提取一个dd得到:
根据我们做莫比乌斯反演题的经验,我们这个时候应该用莫比乌斯函数的一个这样的性质:
∑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=1nd3∑i=1⌊nd⌋∑j=1⌊nd⌋ij∑t|gcd(i,j)μ(t)
先枚举tt可以得到:
然后设sum(x)=∑xi=1isum(x)=∑i=1xi,原式就是:
∑d=1nd3∑t=1⌊nd⌋μ(t)t2sum(⌊nT⌋)2∑d=1nd3∑t=1⌊nd⌋μ(t)t2sum(⌊nT⌋)2
设T=tdT=td可得:
∑T=1nT2sum(⌊nT⌋)2∑d|Tμ(Td)d∑T=1nT2sum(⌊nT⌋)2∑d|Tμ(Td)d
狄利克雷卷积
嗯,推到这一步,推不下去了,怎么办?
看一波题解看一看后面那个∑∑,是不是让你想到莫比乌斯反演的另一个性质呢?不是
∑d|nμ(d)d=ϕ(n)n∑d|nμ(d)d=ϕ(n)n
所以如果设id(x)=xid(x)=x,那么id∗μ=ϕid∗μ=ϕ(∗∗指狄利克雷卷积)
所以原始又华丽丽地变身成为:
由于我们知道,形如⌊nT⌋⌊nT⌋的东西是可以利用数论分块快速的搞出来的,所以我们的主要矛盾就转化为了
杜教筛
首先把杜教筛的套路摆在桌子上(SS表示的前缀和,gg是一个能让求前缀和变快的迷之函数):
欧拉函数……狄利克雷卷积……这两样东西是不是又让你想起了欧拉函数的一条神秘性质?
∑d|nϕ(d)=n∑d|nϕ(d)=n
所以你在思索gg到底是什么的时候,会想到一种可能性
然后搞一搞:
(g∗f)(i)=∑d|id2ϕ(d)i2d2=∑d|iϕ(d)i2=i3(g∗f)(i)=∑d|id2ϕ(d)i2d2=∑d|iϕ(d)i2=i3
生活如此美好!再看看我们杜教筛的时候要搞些什么:
S(n)=∑i=1ni3−∑i=2ni2S(ni)S(n)=∑i=1ni3−∑i=2ni2S(ni)
可是你发现,∑i3∑i3和∑i2∑i2都不是可以快速求的东西啊……这时候你就要想起一样神奇的玩意儿——
“小学奥数”
或者说,打表找规律,因此发现:
∑i=1ni3=(∑i=1ni)2∑i=1ni3=(∑i=1ni)2
这个可以用数学归纳法证明一下,首先对于n=1n=1的情况,显然成立。
然后若对于某个nn成立,那么
然后再打表找规律一波,发现:
∑i=1ni2=n(n+1)(2n+1)6∑i=1ni2=n(n+1)(2n+1)6
对于n=1n=1的情况,显然成立,若对于某个nn成立,那么:
好吧,
于是我们就可以用公式迅速求求求,然后杜教筛一波,在数论分块一波,搞出这道题了。
代码
#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;
}