洛谷P3769 简单的数学题 莫比乌斯反演+杜教筛

题目分析

瞎推公式

让我们来乱推一波公式……首先我们要求：

∑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和jj都提取一个dd得到：

∑d=1nd3∑i=1⌊nd⌋∑j=1⌊nd⌋ij[gcd(i,j)==1]∑d=1nd3∑i=1⌊nd⌋∑j=1⌊nd⌋ij[gcd(i,j)==1]

根据我们做莫比乌斯反演题的经验，我们这个时候应该用莫比乌斯函数的一个这样的性质：

∑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可以得到：

∑d=1nd3∑t=1⌊nd⌋μ(t)t2∑i=1⌊ntd⌋∑i=1⌊ntd⌋ij∑d=1nd3∑t=1⌊nd⌋μ(t)t2∑i=1⌊ntd⌋∑i=1⌊ntd⌋ij

然后设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∗μ=ϕ（∗∗指狄利克雷卷积）

所以原始又华丽丽地变身成为：

∑T=1nT2ϕ(T)sum(⌊nT⌋)2∑T=1nT2ϕ(T)sum(⌊nT⌋)2

由于我们知道，形如⌊nT⌋⌊nT⌋的东西是可以利用数论分块快速的搞出来的，所以我们的主要矛盾就转化为了人民日益增长的美好生活需要和不平衡不充分的发展之间的矛盾快速求f(x)=x2ϕ(x)f(x)=x2ϕ(x)的前缀和，这个应该要用杜教筛。

杜教筛

首先把杜教筛的套路摆在桌子上(SS表示ff的前缀和，gg是一个能让求前缀和变快的迷之函数）：

g(1)S(n)=∑i=1n(g∗f)(i)−∑i=2ng(i)S(ni)g(1)S(n)=∑i=1n(g∗f)(i)−∑i=2ng(i)S(ni)

欧拉函数……狄利克雷卷积……这两样东西是不是又让你想起了欧拉函数的一条神秘性质？不是

∑d|nϕ(d)=n∑d|nϕ(d)=n

所以你在思索gg到底是什么的时候，会想到一种可能性g(x)=x2g(x)=x2

然后搞一搞：

(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=1n+1)2=((i+1)(i+2)2)2=((i+1)i2)2+(i+1)3(∑i=1n+1)2=((i+1)(i+2)2)2=((i+1)i2)2+(i+1)3

然后再打表找规律一波，发现：

∑i=1ni2=n(n+1)(2n+1)6∑i=1ni2=n(n+1)(2n+1)6

对于n=1n=1的情况，显然成立，若对于某个nn成立，那么：

(n+1)(n+2)(2n+3)6=n(n+1)(2n+1)6+(n+1)2(n+1)(n+2)(2n+3)6=n(n+1)(2n+1)6+(n+1)2

好吧，伪证证明完毕。

于是我们就可以用公式迅速求求求，然后杜教筛一波，在数论分块一波，搞出这道题了。

神TM简单的数学题

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL N=8000000;
int pri
,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;
}