[JZOJ4496] 【GDSOI 2016】互补约数

Description

求

∑d=1N∑i|dgcd(i,di)

Solution

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

原式显然可以化。

∑i=1N∑j=1⌊Ni⌋gcd(i,j)

然后就是套路
设

f(d)=∑i=1N∑j=1⌊Ni⌋[gcd(i,j)=d]

g(d)=∑i=1⌊Nd⌋f(id)

根据g的定义（f那个式子中d|i,d|j的组数）
可得

g(T)=∑i=1⌊NT2⌋⌊NT2i⌋

因为d|i,d|j,i∗j<=N
所以d<=N−−√

Ans=∑d=1N√df(d)=∑d=1N√d∑d|Tμ(Td)g(T)

交换主体

=∑T=1N√g(T)∑d|Tμ(Td)d=∑T=1N√g(T)∑d|Tμ(d)Td

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

∑d|Tμ(d)Td=φ(T)

证明懒得打了，具体参见大牛LYD729的BLOG
然后

=∑T=1N√g(T)φ(T)

g带进去

=∑T=1N√φ(T)∑i=1⌊NT2⌋⌊NT2i⌋

然后分块可以在很优的复杂度跑出来

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;}