luogu 4844 LJJ爱数数 (莫比乌斯反演+数学推导)

题目大意：求满足gcd(a,b,c)==1，1/a+1/b=1/c，a,b,c<=n的{a,b,c}有序三元组个数

因为题目里有LJJ我才做的这道题

出题人官方题解https://www.cnblogs.com/Blog-of-Eden/p/9367521.html对我帮助很大

思维很巧妙的一道题，佩服出题人Orzzz

由原式可得，$c=frac{ab}{a+b}$

令g=gcd(a,b)，A=a/g,B=b/g，显然gcd(g,c)==1，gcd(A,B)==1

带入可得$frac{ABg^{2}}{(A+B)g}=c$ <=> $frac{ABg}{A+B}=c$

因为A,B互质，所以$A,B,A+B$两两互质

由$frac{ABg}{A+B}=c$可得

因为c是整数，A与A+B互质，B与$A+B互质，所以当且仅当(A+B)|g

令G=g/(A+B)，可得ABG=c，所以G|c，而g与c互质，所以G作为g的因子，与c也一定互质，即gcd(G,c)==1，所以G只能等于1

综上，可得g=A+B,c=AB,a+b=$g^{2}$

现在大问题转化成了一般性问题，每次枚举一个g，求在一定范围内，g=A+B且gcd(A,B)==1的数对数量

显然这样的数对数量可以用莫比乌斯反演求得

由于g<=$sqrt(2n)$，g<$2*10^{6}$，通过预处理，分解质因数是时间可以优化成logn，再筛出它所有的因子，利用莫比乌斯函数的容斥性质，即可求得合法数对数量

筛出每个数的因子的时间是均摊logn(调和级数)，但由于空间限制，不能预处理出每个数的所有因子..

而A的合法范围则是保证A<=g,A*g<=n&&B*g<=n

细节需要多思考

#include <cmath>
#include <vector>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 2001000
#define maxn 2000000
#define ll long long 
#define uint unsigned int
using namespace std;

ll n;
int cnt;

int pr[N],use[N],mu[N],pmu[N],nxt[N];
void get_pr()
{
    mu[1]=pmu[1]=1;
    for(int i=2;i<=maxn;i++)
    {
        if(!use[i]) pr[++cnt]=i,nxt[i]=i,mu[i]=-1;
        pmu[i]=pmu[i-1]+mu[i];
        for(int j=1;j<=cnt&&i*pr[j]<=maxn;j++){
            use[i*pr[j]]=1,nxt[i*pr[j]]=pr[j];
            if(i%pr[j]==0){
                mu[i*pr[j]]=0;
                break;
            }else{
                mu[i*pr[j]]=-mu[i];
            }
        }
    }
}
int son[N],d[N],ps[N],num,nson;
void dfs_son(int s,int dep)
{
    if(dep>num) {son[++nson]=s;return;}
    for(int j=0;j<=d[dep];j++)
        dfs_son(s,dep+1),s*=ps[dep];
}
void get_son(int x)
{
    num=0;
    while(x!=1){
        int p=nxt[x];ps[++num]=p;d[num]=0;
        while(x%p==0) x/=p,d[num]++; 
    }
    if(x!=1) ps[++num]=x,d[num]=1;
    for(int i=1;i<=nson;i++) son[i]=0;
    nson=0;
    dfs_son(1,1);
}
ll solve(int l,int r,int g)
{
    ll ans=0;l--;
    get_son(g);
    for(int i=1;i<=nson;i++)
        ans+=1ll*mu[son[i]]*(r/son[i]-l/son[i]);
    return ans;
}

int main()
{
    //freopen("t1.in","r",stdin);
    scanf("%lld",&n);
    ll sq=sqrt(n*2);
    get_pr();
    ll ans=0;
    for(ll g=1;g<=sq;g++)
    {
        int l=max(1ll,((g*g-n)/g+( ((g*g-n)%g==0)?0:1 )));
        int r=min(g-1,n/g);
        if(l>r) continue;
        ans+=solve(l,r,g);
    }
    printf("%lld
",ans);
    return 0;
}