【GDSOI 2016】互补约数

Description

Solution

看到对于gcd的求和就知道是关于莫比乌斯繁衍的题目。
先来看一看要求什么：所有数分解成两个相乘的数的gcd，那么我们可以直接枚举那些相乘的数。
所以

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

和一般的莫比乌斯繁衍的长得很像，我们试着繁衍一下。
每次繁衍的最原始函数：f[i]表示有多少个点对的gcd等于i，g[i]表示有多少个点对的gcd是i的倍数。
很显然的初始关系：

g[d]=∑i=1⌊nd⌋f[i∗d]

那么繁衍之后的关系就是：

f[d]=∑i=1⌊nd⌋g[i∗d]∗μ[i]

我们来观察一下ans与f数组的关系，因为是求和就用出现次数(f[i])*数字(i)就可以了：

ans=∑i=1nf[i]∗i

把繁衍式子代入：

ans=∑i=1ni∑d=1⌊ni⌋g[i∗d]∗μ[d]

一般式子里不允许一个数组里面有两个变量组成，我们用T=i∗d来枚举T：

ans=∑i=1ni∑d=1⌊ni⌋g[T]∗μ[d]

这样考虑直接枚举T才方便做,明显发现μ和i的出现，都是在是T的因数时且i∗d=T：

ans=∑T=1ng[T]∑i|Tμ[i]∗Ti

但是这样枚举的复杂度还是没有太大的变化，我们想原来繁衍那样观察一下g数组，g[i]是所有表示为(i*x,i*y)点对的gcd的和构成的，那么直接枚举x就可以算出y的个数：
因为i∗x∗i∗y≤n，所以xy≤⌊ni2⌋，那么x最大为⌊ni2⌋
所以

g[i]=∑x=1⌊ni2⌋⌊ni2x⌋

带入ans里面：

ans=∑T=1n∑x=1⌊nT2⌋⌊nT2x⌋∑i|Tμ[i]∗Ti

很显然可以发现，当T大于n√时，g(T)=0，所以T只用枚举到n√就可以了：

ans=∑T=1n√∑x=1⌊nT2⌋⌊nT2x⌋∑i|Tμ[i]∗Ti

发现后面关于μ[i]的地方只与i的因数有关，就是说只用枚举两个数相乘等于i然后算答案，预处理出a[i]=∑i|Tμ[i]∗Ti

fo(i,1,gen){        fo(j,1,gen/i)a[i*j]=a[i*j]+miu[i]*j;    }

那么答案就变成：

ans=∑T=1n√a[T]∑x=1⌊nT2⌋⌊nT2x⌋

我们设t=⌊nT2⌋，答案就变成：

ans=∑T=1n√a[T]∑x=1t⌊tx⌋

那么只用分块来做后面就好了。

Code

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