【GDSOI 2016】互补约数

Description

Solution1

一题标准的莫比乌斯反演，
变式：

$$ans=\sum_{i=1}^n\sum_{j|i}gcd(i,\frac{j}{i})=\sum_{i=1}^n\sum_{j=1}^{\lfloor\frac{n}{i}\rfloor}gcd(i,j)$$
设gcd(i,j)=d的数有 $f(d)$个， $g(d)$表示 $d|gcd(i,j)$的个数,
$$g(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(id)$$ ，
根据反演：
$$f(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(id)\mu(i)$$
那么
$$ans=\sum_{d=1}^{\sqrt{n}}f(d)*d$$
（因为gcd(i,j)不可能大于 $\sqrt n$）
展开：
$$ans=\sum_{d=1}^{\sqrt{n}}d*\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}g(id)\mu(i)$$
用 $T=id$
$$ans=\sum_{T=1}^{n}g(T)\sum_{d|T}d*\mu(\frac{T}{d})$$
因为 $g(d)$表示 $d|gcd(i,j)$的个数，所有可以表示成： $gcd(d*x,d*y)$，
枚举其中的x，
$$g(T)=\sum_{x=1}^{\lfloor\frac{n}{T^2}\rfloor}\lfloor\frac{n}{T^2*x}\rfloor$$
所以当 $T>\sqrt n$时， $g(T)=0$无意义，第一个枚举的界限可以降低到 $\sqrt n$
$$ans=\sum_{T=1}^{\sqrt n}g(T)\sum_{d|T}d*\mu(\frac{T}{d})$$
$g(T)$用分块，后面那一堆暴力算，

复杂度：约$O(n^{3\over 4})$（我不会算，CZY用微积分证明了(%%%)）

Solution2

$$ans=\sum_{T=1}^{\sqrt n}\varphi(T)\sum_{x=1}^{\lfloor\frac{n}{T^2}\rfloor}\lfloor\frac{n}{T^2*x}\rfloor$$
因为 $\varphi(n)=\sum_{d|n}\frac{n}{d}*\mu(d)$;（ 详情点这里）
所以我们发现跟上面的式子一模一样。

Code

#include<cstdio>
#include<cmath>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long LL;
const int N=400000;
int pr[N/3],mu[N];
int b[N/3];
bool prz[N];
LL g[N],n,m,ans;
void dofirst(int n)
{
    mu[1]=1;
    fo(i,2,m)
    {
        if(!prz[i])pr[++pr[0]]=i,mu[i]=-1;
        for(int j=1;j<=pr[0] && i*pr[j]<=m;j++)
        {
            prz[i*pr[j]]=1;
            if(!(i%pr[j]))break;
            mu[i*pr[j]]=-mu[i];
        }
    }
}
void zh(int q,int e,int T)
{
    if(q>pr[0]){g[T]+=(LL)e*mu[T/e];return;}
    fo(i,0,b[q])zh(q+1,e,T),e*=pr[q];
}
LL fk(LL n)
{
    LL ans=0,t=1,r;
    while(t<=n)
    {
        r=n/(n/t);
        ans+=n/t*(r-t+1);
        t=r+1;
    }
    return ans;
}
int main()
{
    freopen("gcd.in","r",stdin);
    freopen("gcd.out","w",stdout);
    scanf("%lld",&n);m=sqrt(n);
    dofirst(m);
    fo(I,1,m)
    {
        pr[0]=0;int i=I;
        fo(j,2,sqrt(i))if(!(i%j))
        {
            pr[++pr[0]]=j;b[pr[0]]=0;
            while(!(i%j))i/=j,b[pr[0]]++;
        }
        if(i-1)pr[++pr[0]]=i,b[pr[0]]=1;
        zh(1,1,I);
    }
    fo(T,1,m)ans+=g[T]*fk(n/T/T);
    printf("%lld\n",ans);
    return 0;
}