[JZOJ4496] 【GDSOI 2016】互补约数

Description

求

$$\sum\limits_{d=1}^{N}\sum\limits_{i|d} gcd(i,{d\over i})$$

Solution

LF同志说过一句话。
“看到GCD，我们就要坚信它能够反演！”

原式显然可以化。

$$\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{\lfloor{N\over i}\rfloor} gcd(i,j)$$

然后就是套路
设

$$f(d)=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{\lfloor{N\over i}\rfloor} [gcd(i,j)=d]$$
$$g(d)=\sum\limits_{i=1}^{\lfloor {N\over d }\rfloor}f(id)$$
根据g的定义（f那个式子中 $d|i,d|j$的组数）
可得 $$g(T)=\sum\limits_{i=1}^{\lfloor{N\over T^2}\rfloor}\lfloor{N\over T^2i}\rfloor$$

因为$d|i,d|j,i*j<=N$
所以$d<=\sqrt N$

$$Ans=\sum\limits_{d=1}^{\sqrt N} df(d)=\sum\limits_{d=1}^{\sqrt N} d\sum\limits_{d|T} \mu({T\over d})g(T)$$
交换主体
$$=\sum\limits_{T=1}^{\sqrt N} g(T)\sum\limits_{d|T} \mu({T\over d}) d=\sum\limits_{T=1}^{\sqrt N} g(T)\sum\limits_{d|T} \mu(d){T\over d}$$

显然，你可以分块再预处理一些东西，把g代进去即可
然而，此处有一个性质

$$\sum\limits_{d|T} \mu(d){T\over d}=φ(T)$$
证明懒得打了，具体参见 大牛LYD729的BLOG
然后
$$=\sum\limits_{T=1}^{\sqrt N} g(T)φ(T)$$
g带进去
$$=\sum\limits_{T=1}^{\sqrt N} φ(T)\sum\limits_{i=1}^{\lfloor{N\over T^2}\rfloor}\lfloor{N\over T^2i}\rfloor$$
然后分块可以在很优的复杂度跑出来

Code

#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <iostream>
#include <cstring>
#include <cmath>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fod(i,a,b) for(int i=a;i>=b;i--)
#define N 320000
#define LL long long 
using namespace std;
LL pr[N],phi[N],sum[N],n;
int m,l;
bool bz[N];
void prp()
{
    l=0;
    phi[1]=sum[1]=1;
    fo(i,2,m)
    {
        if(!bz[i]) pr[++l]=i,phi[i]=i-1;
        for(int j=1;j<=l&&pr[j]*i<=m;j++)
        {
            bz[i*pr[j]]=1;
            if(i%pr[j]==0)
            {
                phi[i*pr[j]]=phi[i]*pr[j];
                break;
            }
            phi[i*pr[j]]=phi[i]*(pr[j]-1);
        }
        sum[i]=sum[i-1]+phi[i];
    }
}
int main()
{
    cin>>n;
    m=sqrt(n);
    prp();
    LL T=1,s=0;
    while(T<=m)
    {
        LL x=1,lim=n/T/T,s1=0;
        LL T1=sqrt(n/lim);
        while(x<=lim)
        {   
            LL x1=lim/(lim/x);
            s1+=(x1-x+1)*(lim/x);
            x=x1+1;
        }
        s+=(sum[T1]-sum[T-1])*s1;
        T=T1+1;
    }
    cout<<s;
}