【JZOJ 4496】【GDSOI 2016】第一题 互补约数 （两种解法）

Description

求∑i=1n∑d|igcd(d,id)

对于100%的数据，n≤1011。

Analysis

方法一

这是一个不需要莫比乌斯反演的方法。

∑i=1n∑d|igcd(d,id)

=∑i∗j≤ngcd(i,j)

=∑i∗j≤n∑d|gcd(i,j)φ(d)

这里用到了一个φ的性质，∑d|nφ(d)=n，具体证明可以看我的博客。

=∑i∗j≤n∑d|i⋂d|jφ(d)

调换主体，得

=∑d=1n√φ(d)∑i∗j≤⌊nd2⌋1

=∑d=1n√φ(d)∑i=1⌊nd2⌋⌊nd2i⌋

随着d的增长，d2增长迅速，因此⌊nd2⌋会迅速收敛到很小，并且有许多的⌊nd2i⌋的值是一样的，所以可以分块处理。时间复杂度是很优的。

方法二

这次我们用莫比乌斯反演。

∑i=1n∑d|igcd(d,id)

=∑i=1n∑j=1⌊ni⌋gcd(i,j)

设f(d)表示上式中满足gcd(i,j)=d的i,j有多少对;

g(d)表示上式中满足d|gcd(i,j)的i,j有多少对。

首先有

g(d)=∑i=1⌊nd⌋f(d∗i)

反演

f(d)=∑i=1⌊nd⌋g(d∗i)∗μ(i)

而我们要求的Ans=∑d=1nd∗f(d)

=∑d=1nd∑i=1⌊nd⌋g(d∗i)∗μ(i)

设T=d∗i，再调换主体

=∑T=1ng(T)∑i|TTi∗μ(i)

但其实T只用枚举到n−−√，因为n−−√是在i∗j≤n的范围内gcd(i,j)的最大取值，此时i=j=n−−√

那么g(T)怎么求呢？

若d|gcd(i,j)，则可以设dx,dy分别为任意的i,j，因为必有d|gcd(dx,dy)

所以我们枚举x，因为dx∗dy≤n，所以x只需要枚举到⌊nd2⌋

g(T)=∑x=1⌊nT2⌋⌊nT2x⌋

所以Ans=∑T=1n√∑x=1⌊nT2⌋⌊nT2x⌋∗∑i|TTi∗μ(i)

这个东西直接做，预处理一些东西，就可以AC了。

但我们的追求更远。

上面这个式子有没有一点眼熟？

我们把方法一最后得出的结论搬下来

Ans=∑d=1n√φ(d)∑i=1⌊nd2⌋⌊nd2i⌋，

两个式子比较一下，许多东西非常相似，所以我们可以提出这样一个猜想

φ(n)=∑d|nnd∗μ(d)

证明如下：

设f(n)=φ(n),g(n)=n

根据上面说的那个欧拉函数的性质，得

g(n)=n=∑d|nφ(d)=∑d|nf(d)

反演

f(n)=φ(n)=∑d|ng(nd)∗μ(d)=∑d|nnd∗μ(d)

得证！

殊途同归！

Code

#include<cstdio>
#include<cmath>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=316230;
int phi
,pri
;
bool bz
;
void pre(int n)
{
phi[1]=1;
fo(i,2,n)
{
if(!bz[i])
{
phi[i]=i-1;
pri[++pri[0]]=i;
}
fo(j,1,pri[0])
{
int t=i*pri[j];
if(t>n) break;
bz[t]=1;
if(i%pri[j]==0)
{
phi[t]=phi[i]*pri[j];
break;
}
phi[t]=phi[i]*phi[pri[j]];
}
}
}
int main()
{
freopen("gcd.in","r",stdin);
freopen("gcd.out","w",stdout);
ll n,ans=0;
scanf("%lld",&n);
int m=int(sqrt(n));
pre(m);
fo(p,1,m)
{
ll i,l=n/p/p,t=0;
for(ll j=1;j<=l;j=i+1)
{
i=l/(l/j);
t+=(i-j+1)*(l/j);
}
ans+=phi[p]*t;
}
printf("%lld",ans);
fclose(stdin);fclose(stdout);
return 0;
}