【BZOJ2671】—Calc（莫比乌斯反演）

传送门

考虑实际求
$ans=sum_{i=1}^{n}sum_{j=i+1}^n[i+j|i*j]$

令$d=gcd(i,j),i=i'd,j=j'd$

$i+j|i*j ightarrow d(i'+j')|d^2i'j'$

则$i'+j'|di'j'$

由于$gcd(i',j')=1$所以$gcd(i'+j',i'j')=1$
所以也要求的就是$i'+j'|d$

$sum_{i=1}^nsum_{j=i+1}^{n}sum_{d=1}^{n/j}[u+v|d][gcd(u,v)=1]$
$=sum_{i=1}^nsum_{j=i+1}^{n}frac{n}{v(u+v)}[gcd(u,v)=1]$

反演一下就是

$sum_{d}mu(d)sum_{i,j}frac{n}{max(i,j)*(i+j)d^2}$

枚举$k=i+j$钦定一个大小顺序

$sum_{d}mu(d)sum_{k=1}sum_{i=k/2+1}^{k-1}frac{n}{ikd^2}$

然后卡一下上界枚举$d,k$，对$i$整除分块
复杂度据说是$O(n^{frac 3 4})$

#include<bits/stdc++.h>
using namespace std;
#define gc getchar
inline int read(){
    char ch=gc();
    int res=0,f=1;
    while(!isdigit(ch))f^=ch=='-',ch=gc();
    while(isdigit(ch))res=(res+(res<<2)<<1)+(ch^48),ch=gc();
    return f?res:-res;
}
#define re register
#define pb push_back
#define cs const
#define pii pair<int,int>
#define fi first
#define se second
#define ll long long
#define poly vector<int>
#define bg begin
cs int mod=1e9+7;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline void Add(int &a,int b){(a+=b)>=mod?(a-=mod):0;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline void Dec(int &a,int b){(a-=b)<0?(a+=mod):0;}
inline int mul(int a,int b){return 1ll*a*b>=mod?1ll*a*b%mod:a*b;}
inline void Mul(int &a,int b){a=mul(a,b);}
inline int ksm(int a,int b,int res=1){for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));return res;}
inline void chemx(ll &a,ll b){a<b?a=b:0;}
inline void chemn(int &a,int b){a>b?a=b:0;}
cs int N=100005;
int mu[N],pr[N],tot;
bitset<N> vis;
inline void init(){
    mu[1]=1;
    for(int i=2;i<N;i++){
        if(!vis[i])pr[++tot]=i,mu[i]=-1;
        for(int j=1;i*pr[j]<N&&j<=tot;j++){
            vis[i*pr[j]]=1;
            if(i%pr[j]==0)break;
            mu[i*pr[j]]=-mu[i];
        }
    }
}
ll res,n;
int main(){
    init();
    n=read();
    for(ll d=1;d*d<=n;d++){
        ll now=0;
        for(ll k=1,lim=min(n/d/d,(ll)sqrt(n*2/d/d));k<=lim;k++){
            ll mx=n/d/d/k;
            for(ll j=k/2+1,nxt;j<k&&j<=mx;j=nxt+1){
                nxt=min(k-1,mx/(mx/j));
                now+=(nxt-j+1)*(mx/j);
            }
        }
        res+=now*mu[d];
    }
    cout<<res;
}