题目链接

题目大意

给出n，求出满足以下两个条件的数对(a,b)的数量
$1. 1\le a<b \le n$

$2. (a+b)|(a*b)$

做法

对于数对$(a,b)$，令$g=gcd(a,b)$，再令$p*g=a,q*g=b$
那么有$(p*g+q*g)|(p*g*q*g)$
两边同时约掉一个g，$(p+q)|(p*q+g)$
由于$p,q$互质，所以$(p+q)$与$(p*q)$也一定互质，这个结论应该还是比较显然的。
所以还要满足整除关系的话，$g$一定是$(p+q)$的倍数，我们令$g=k*(p+q)$
别忘了我们还有第一个限制条件，我们将b代入。
$b=g*q=k*(p+q)*q\le n$
显然$q\le \sqrt n$
令$m=\sqrt{n}$
枚举$q$，可以得到

$ans=\sum_{q=1}^{m}\sum_{p=q+1}^{m}[gcd(p,q)==1]*\lfloor\frac{n}{y(x+y)}\rfloor$

看到这个$[gcd(x,y)==1]$，想到莫比乌斯反演

$ans=\sum_{q=1}^{m}\sum_{p=q+1}^{m}\lfloor\frac{n}{y(x+y)}\rfloor*\sum_{d|gcd(x,y)}\mu(d)$

还是老套路，改为枚举$d$和倍数

$ans=\sum_{d=1}^{m}\mu(d)\sum_{i=1}^{\lfloor\frac{m}{d}\rfloor}\sum_{j=1}^{i-1}\lfloor\frac{n}{d*d(i+j)*i}\rfloor$

再进行一步简单的变形，改为枚举$i+j$和$i$

$ans=\sum_{d=1}^{m}\mu(d)\sum_{i=2}^{2*\lfloor\frac{m}{d}\rfloor-1}\sum_{j=\frac{i}{2}+1}^{i-1}\lfloor\frac{n}{d*d*i*j}\rfloor$

我们先枚举$d$，再枚举$i$，根据调和级数，这一步的复杂度是$O(m*logm)$的，有了$d$和$i$了之后，再对$j$整除分块，复杂度大概是$O(m^{\frac{3}{2}}*logm)$，也就是$O(n^{\frac{3}{2}}*log\sqrt n)$

代码

#include<cstdio>
#include<algorithm>
#include<cctype>
#include<cstring>
#include<iostream>
#include<cmath>
#define LL long long
#define N (100005)
using namespace std;
bool nprime[N];
int pri,r;
LL ans;
int prime[N],mu[N];
template <typename T> void read(T&t) {
    t=0;
    bool fl=true;
    char p=getchar();
    while (!isdigit(p)) {
        if (p=='-') fl=false;
        p=getchar();
    }
    do {
        (t*=10)+=p-48;p=getchar();
    }while (isdigit(p));
    if (!fl) t=-t;
}
int n,m;
int main(){
    read(n);
    m=sqrt(n);
    mu[1]=1;
    for (int i=2;i<=m;i++){
        if (!nprime[i]){
            prime[++pri]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=pri&&i*prime[j]<=m;j++){
            nprime[i*prime[j]]=1;
            if (i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
    }
    for (int d=1;d<=m;d++){
        if (mu[d]==0) continue;
        for (int i=2;i<=2*(m/d)-1;i++){
            int t=n/d/d/i;
            if (t==0) break;
            //for (int j=i/2+1;j<i;j++) ans+=mu[d]*(t/j);
            for (int j=i/2+1;j<i;j=r+1){
                if (t/j==0) break;
                r=min(i-1,t/(t/j));
                ans+=(r-j+1)*(t/j)*mu[d];
            }
        }
    }
    printf("%lld",ans);
    return 0;
}