LJJ爱数数

LJJ爱数数

求\(\sum_{i=1}^n\sum_{j=1}^n\sum_{k=1}^n\epsilon(gcd(i,j,k))(\frac{1}{i}+\frac{1}{j}==\frac{1}{k}),n\leq 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)=1\Rightarrow 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+j\Rightarrow g\leq\sqrt{2n}\Rightarrow k\leq n \]

现在考虑\(I\)

\[1\leq i=Ig,j=(g-I)g\leq n\Downarrow \]

\[I\in[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)=1\Rightarrow 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,首先显然\(m\in[1,\sqrt{N}]\)

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

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

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

---

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

但是因为\(i,j\in Z\),所以必然存在一个\(z\in 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;
        }
    }
}