LJJ爱数数

求(sum_{i=1}^nsum_{j=1}^nsum_{k=1}^nepsilon(gcd(i,j,k))(frac{1}{i}+frac{1}{j}==frac{1}{k}),nleq 10^{12})

解

显然无法用Mobius反演，问题在于(frac{1}{i}+frac{1}{j}==frac{1}{k})，要将其转换为gcd条件。

法一：先约数拆分，再证明对应相等

分数我们无法处理,所以有

[(i+j)k=ij ]

设(g=gcd(i,j),I=i/g,J=j/g),接着有

[(I+J)k=IJg,gcd(I,J)=1 ]

[gcd(i,j,k)=gcd(g,k)=1 ]

考虑分数反证,设

[frac{I+J}{g}=frac{IJ}{k}=d ]

(if(d>1))

[d|IJ,d|(I+J) ]

而

[gcd(I,J)=1Rightarrow gcd(I+J,IJ)=1 ]

所以与原命题矛盾，故(d=1),因此(I+J=g,IJ=k,gcd(g,I)=1),满足了等式，也满足了(gcd(g,k)=1),所以自然枚举g，再枚举I,以此k也就确定了，接下来问题在于范围，显然要回到(i,j,k)

[g^2=i+jRightarrow gleqsqrt{2n}Rightarrow kleq n ]

现在考虑(I)

[1leq i=Ig,j=(g-I)gleq nDownarrow ]

[Iin[max(1,frac{g^2-n}{g}),min([n/g],g-1)] ]

于是总上，设(L(g),R(g))为在g意义下(I)的范围,有

[ans=sum_{g=1}^{sqrt{2n}}sum_{I=L(g)}^{R(g)}(gcd(g,I)==1) ]

注意到后式取值范围受到g的影响，所以对后式单独维护

设

[f(d)=sum_{I=L(g)}^{R(g)}(gcd(g,I)==d) ]

[F(d)=([R(g)/d]-[(L(g)-1)/d])(d|g) ]

由Mobius反演定理，我们有

[f(d)=sum_{d|x}(x|g)([R(g)/x]-[(L(g)-1)/x])mu(x/d) ]

因此

[f(1)=sum_{x|g}([R(g)/x]-[(L(g)-1)/x])mu(x) ]

所以

[ans=sum_{g=1}^{sqrt{2n}}sum_{x|g}([R(g)/x]-[(L(g)-1)/x])mu(x) ]

[=sum_{x=1}^{sqrt{2n}}mu(x)sum_{x|g}^{^{sqrt{2n}}}([R(g)/x]-[(L(g)-1)/x]) ]

显然后式可以暴力算，时间复杂度应为(O(sqrt{n}log(sqrt{n})))。

参考代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#define il inline
#define ri register
#define ll long long
using namespace std;
bool check[2000000];
int prime[150000],pt,mu[2000000],
    R[2000000],L[2000000];
il void prepare(int);
template<class free>il free Min(free,free);
template<class free>il free Max(free,free);
int main(){
    ll n,sn,i,j,ans(0);
    scanf("%lld",&n),sn=sqrt(n<<1),prepare(sn);
    for(i=1;i<=sn;++i)L[i]=Max(i-n/i,1ll),
                          R[i]=Min(i-1,n/i);
    for(i=1;i<=sn;++i)
        for(j=i;j<=sn;j+=i)
            ans+=mu[i]*(R[j]/i-(L[j]-1)/i);
    printf("%lld",ans);
    return 0;
}
il void prepare(int n){
    int i,j;check[1]|=mu[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mu[i]=-1;
        for(j=1;j<=pt&&i*prime[j]<=n;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mu[i*prime[j]]=~mu[i]+1;
        }
    }
}
template<class free>il free Max(free a,free b){return a>b?a:b;}
template<class free>il free Min(free a,free b){return a<b?a:b;}

法二：先证明对应相等再约数拆分

注意到(ik+jk=ij),类似因式分解形式，于是添项因式分解，

[k^2+ij-ik-jk=k^2 ]

[k(k-i)+j(i-k)=k^2 ]

[(j-k)(i-k)=k^2 ]

再进行约数拆分

[j-k=ln^2,i-k=lm^2,k=lmn ]

[gcd(i,j,k)=1Rightarrow l=1,j-k=n^2,i-k=m^2,k=mn ]

[j=mn+n^2,i=mn+m^2,k=mn ]

[gcd(mn+n^2,mn+m^2,mn)=gcd(gcd(m,n)(m+n),mn)=1 ]

[gcd(m,n)=1 ]

故充分条件为(j-k=n^2,i-k=m^2,gcd(m,n)=1),故考虑枚举n，m，为了与题目数据范围的n区分，把数据范围的n改为N,现在关键在于找数据范围，回到i，j,首先显然(min[1,sqrt{N}])

[j=mn+n^2,i=mn+m^2,k=mnin[1,N]Downarrow ]

[nleq min(frac{-m+sqrt{m^2+4N}}{2},frac{N-m^2}{m}) ]

但事实上，我们只能枚举整数，但是m,n有可能为无理数，但是我们可以证明这样没有影响

设(m=sqrt{x},n=sqrt{y}(x,yin Z)),于是有(i=sqrt{xy}+x,j=sqrt{xy}+y,k=sqrt{xy})

但是因为(i,jin Z),所以必然存在一个(zin Z)，满足(xy=z^2)，容易知道x,y,z两两不互质，所以有

(i=x+z,j=y+z,k=z),容易知道(i,j,k)不互质，那么与题设矛盾，于是不可能。

设其范围为(r(m)),于是有

[ans=sum_{i=1}^{sqrt{N}}sum_{j=1}^{r(i)}epsilon(gcd(i,j)) ]

同法一，后式是变的，维护后式，设

[f(d)=sum_{j=1}^{r(i)}(gcd(i,j)==d) ]

[F(d)=(d|i)[r(i)/d] ]

由Mobius反演定理，我们有

[f(d)=sum_{d|x,x|i}[r(i)/x]mu(x/d) ]

[f(1)=sum_{x|i}[r(i)/x]mu(x) ]

代入原式即

[ans=sum_{i=1}^{sqrt{N}}sum_{x|i}[r(i)/x]mu(x)=sum_{x=1}^{sqrt{N}}mu(x)sum_{x|i}[r(i)/x] ]

显然后式处理平均只要(log(sqrt(N))),所以时间复杂度应为(sqrt{N}log(sqrt{N}))。

参考代码：

#include <iostream>
#include <cstdio>
#include <cmath>
#define il inline
#define ri register
#define ll long long
using namespace std;
bool check[1000001];
int prime[100001],pt,mu[1000001],
    li[1000001];
il void prepare(int);
template<class free>
il free Min(free,free);
int main(){
    ll n,ans(0),sn,i,j;
    scanf("%lld",&n),sn=sqrt(n),prepare(sn);
    for(i=1;i<=sn;++i)li[i]=Min((n-i*i)/i,(ll)(sqrt(i*i+4*n)-i)/2);
    for(i=1;i<=sn;++i)
        for(j=i;j<=sn;j+=i)
            ans+=mu[i]*(li[j]/i);
    printf("%lld",ans);
    return 0;
}
template<class free>
il free Min(free x,free y){
    return x<y?x:y;
}
il void prepare(int n){
    int i,j;check[1]|=mu[1]|=true;
    for(i=2;i<=n;++i){
        if(!check[i])prime[++pt]=i,mu[i]=-1;
        for(j=1;j<=pt&&i*prime[j]<=n;++j){
            check[i*prime[j]]|=true;
            if(!(i%prime[j]))break;
            mu[i*prime[j]]=~mu[i]+1;
        }
    }
}