P4844 LJJ爱数数

题目

P4844 LJJ爱数数

本想找到莫比乌斯反演水题练练，结果直接用了两个多小时才做完

做法

(sumlimits_{a=1}^nsumlimits_{b=1}^nsumlimits_{c=1}^n[gcd(a,b,c)=1&&frac{1}{a}+frac{1}{b}=frac{1}{c}])

([gcd(a,b,c)=1])这个好理解，但后面(frac{1}{a}+frac{1}{b}=frac{1}{c})怎么办呢？

下意识去掉分数：((a+b)c=ab)

设(g=gcd(a,b),A=frac{a}{g},B=frac{b}{g})

原式化为：((A+B)c=ABg)

( herefore frac{(A+B)c}{g}=AB)，(g)要整除((A+B)c)

由于(gcd(a,b,c)=1)，(g)不整除(c)，所以(g)整除(A+B)

( herefore frac{A+B}{g}=frac{AB}{c})，设(p=frac{A+B}{g}=frac{AB}{c})

(c=frac{A+B}{g} Longrightarrow p)整除同时整除(A,B)

(g=gcd(a,b),A=frac{a}{g},B=frac{b}{g} Longrightarrow A,B)互质

( herefore p=1)

( herefore A+B=g Longrightarrow a+b=g^2,c=AB Longrightarrow c=frac{ab}{g^2})

接下来就好做了嘛，由于(c=frac{ab}{g^2})的限制，g的上限瞬间就特别小了，所以我们枚举(g)

[egin{aligned} Ans & =sumlimits_{g=1}^{sqrt {2n}}sumlimits_{i=1}^{frac{n}{g}}[gcd(ig,g^2-ig)=g]\ &=sumlimits_{g=1}^{sqrt {2n}}sumlimits_{i=1}^{frac{n}{g}}[gcd(i,g-i)=1]\ &=sumlimits_{g=1}^{sqrt {2n}}sumlimits_{i=1}^{frac{n}{g}}[gcd(i,g)=1]\ end{aligned}]

其实前面(i)的范围并不精确：(1<=g^2-ig<=n Longrightarrow g-frac{n}{g}<=i<=g-1)

再结合之前的范围：(max(1,g-frac{n}{g})<=i<=min(g-1,frac{n}{d}))

原式变为：$$sumlimits_{g=1}^{sqrt{2n}}sumlimits_{i=max(1,g-frac{n}{g})}^{min(g-1,frac{n}{d})}[gcd(i,g)=1]$$

后半部分就是喜闻乐见的莫比乌斯反演了：

[egin{aligned} sumlimits_{i=1}^k[gcd(i,g)=1]&=sumlimits_{i=1}^ksumlimits_{j|gcd(i,g)}mu(j)\ &=sumlimits_{i=1}^ksumlimits_{j=1}^gmu(j)[j|g][j|]\ &=sumlimits_{j|g}mu(j)frac{k}{j}\ end{aligned}]

My complete code

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
#include<vector>
#include<cmath>
using namespace std;
typedef long long LL;
const int maxn=15000000;
inline LL Read(){
	LL x(0),f(1); char c=getchar();
	while(c<'0'||c>'9'){
		if(c=='-') f=-1; c=getchar();
	}
	while(c>='0'&&c<='9')
	    x=(x<<3)+(x<<1)+c-'0',c=getchar();
	return x*f;
}
int mu[maxn],prime[maxn];
bool visit[maxn];
inline void F_phi(int N){
	mu[1]=1; int tot(0);
	for(int i=2;i<=N;++i){
		if(!visit[i])
			prime[++tot]=i,
			mu[i]=-1;
		for(int j=1;j<=tot&&i*prime[j]<=N;++j){
			visit[i*prime[j]]=true;
			if(i%prime[j]==0)
			    break;
			else
				mu[i*prime[j]]=-mu[i];
		}
	}
}
LL Up;
int n;
struct node{
	int val,next;
}dis[maxn];
int num;
int head[maxn];
inline void Add(int u,int val){
	dis[++num]=(node){val,head[u]},head[u]=num;
}
inline LL Get(int d,int k){
	LL ret(0);
	for(int i=head[d];i&&abs(dis[i].val)<=k;i=dis[i].next)
	    ret+=k/dis[i].val;
	return ret;
}
int main(){
	Up=Read();
	n=sqrt(2*Up);
	F_phi(n);
	for(int i=1;i<=n;++i)
	    for(int j=1;1ll*j*i<=1ll*n;++j)
	        if(mu[j])
	            Add(i*j,mu[j]*j);
	LL ans(0);
	for(int g=1,l,r;g<=n;++g){
		l=max(1ll*1,g-Up/g),r=min(1ll*(g-1),Up/g);
		ans+=Get(g,r)-Get(g,l-1);
	}
	printf("%lld",ans);
	return 0;
}