【LOJ6482】LJJ 愛數數 數論

題目大意

　　給你 \(n\)，求
\[ \sum_{a=1}^n\sum_{b=1}^n\sum_{c=1}^n[\frac{1}{a}+\frac{1}{b}=\frac{1}{c}][\gcd(a,b,c)=1]\\ \]
　　\(n\leq {10}^{12}\)

題解

\[ \begin{align} &\sum_{a=1}^n\sum_{b=1}^n\sum_{c=1}^n[\frac{1}{a}+\frac{1}{b}=\frac{1}{c}][\gcd(a,b,c)=1]\\ =&\sum_{a=1}^n\sum_{b=1}^n\sum_{c=1}^n[c(a+b)=ab][\gcd(a,b,c)=1]\\ =&\sum_{a=1}^n\sum_{b=1}^n[(a+b)\mid ab][\gcd(a,b,c)=1]\\ \end{align} \]

　　經過打表能夠發現，一對數 \(a,b(a\leq b)\) 知足條件的充要條件是 \(b\leq n\) 且 \(\frac{a}{\gcd(a,b)}+\frac{b}{\gcd(a,b)}=\gcd(a,b)\)

　　證實：

　　若 \(\gcd(a,b)=1\)，則 \(\gcd(a+b,ab)=1\)

　　記 \(g=\gcd(a,b),a=ga',b=gb'\)，則
\[ a+b=g(a'+b')\\ ab=g^2a'b'\\ \frac{ab}{a+b}=\frac{g^2a'b'}{g(a'+b')}=\frac{ga'b'}{a'+b'}\\ \]
　　因此 \((a'+b')\mid g\)

　　若 \(a'+b'\neq g\)，則 \(\gcd(a,b,c)=\frac{g}{a'+b'}\)，因此 \(a'+b'=g\) 因此 \(a=a'(a'+b'),b=b'(a'+b'),c=a'b'\)

解法一

　　記 \(x=a',y=b'\)
\[ \begin{align} ans&=\sum_{x=1}^n\sum_{y=x}^n[xy+y^2\leq n]\gcd(x,y)=1\\ &=\sum_{d=1}^\sqrt{n}\mu(d)\sum_{x=1}\sum_{y=x}[xy+y^2\leq \frac{n}{d^2}]\\ &=\sum_{d=1}^\sqrt{n}\mu(d)\sum_{y=1}\min\left(\left\lfloor\frac{\left\lfloor\frac{n}{d^2}\right\rfloor-y^2}{y}\right\rfloor,y\right) \end{align} \]
　　時間複雜度：\(O(\sqrt n\log n)\)

解法二

　　\(a+b=g^2\)

　　枚舉 \(g\)，那麼 \(\gcd(g,a')=1\)

　　當 \(g\) 比較小的時候 \((1\leq g\leq \sqrt n)\)，\(a'\) 有 \(\varphi(g)\) 種取值，能夠直接篩

　　當 \(g\) 比較大的時候 \((\sqrt n<g<\sqrt{2n})\)，
\[ \sum_{i=1}^\frac{n}{g}[\gcd(g,i)=1]\\ =\sum_{d\mid g}\mu(d)\lfloor\frac{n}{gd}\rfloor \]
　　能夠暴力枚舉因子。

　　時間複雜度：\(O(\sqrt n\log n)\)

代碼

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cstdlib>
#include<ctime>
#include<utility>
#include<functional>
#include<cmath>
#include<vector>
//using namespace std;
using std::min;
using std::max;
using std::swap;
using std::sort;
using std::reverse;
using std::random_shuffle;
using std::lower_bound;
using std::upper_bound;
using std::unique;
using std::vector;
typedef long long ll;
typedef unsigned long long ull;
typedef double db;
typedef std::pair<int,int> pii;
typedef std::pair<ll,ll> pll;
void open(const char *s){
#ifndef ONLINE_JUDGE
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
void open2(const char *s){
#ifdef DEBUG
    char str[100];sprintf(str,"%s.in",s);freopen(str,"r",stdin);sprintf(str,"%s.out",s);freopen(str,"w",stdout);
#endif
}
int rd(){int s=0,c,b=0;while(((c=getchar())<'0'||c>'9')&&c!='-');if(c=='-'){c=getchar();b=1;}do{s=s*10+c-'0';}while((c=getchar())>='0'&&c<='9');return b?-s:s;}
void put(int x){if(!x){putchar('0');return;}static int c[20];int t=0;while(x){c[++t]=x%10;x/=10;}while(t)putchar(c[t--]+'0');}
int upmin(int &a,int b){if(b<a){a=b;return 1;}return 0;}
int upmax(int &a,int b){if(b>a){a=b;return 1;}return 0;}
int gcd(int a,int b)
{
    return b?gcd(b,a%b):a;
}
int check(int x,int y)
{
    return x*y%(x+y)==0;
}
const int N=1000010;
int b[N];
int pri[N];
int miu[N];
int cnt;
int main()
{
//  open("loj6482");
    ll n;
    scanf("%lld",&n);
    miu[1]=1;
    for(int i=2;i<=1000000;i++)
    {
        if(!b[i])
        {
            pri[++cnt]=i;
            miu[i]=-1;
        }
        for(int j=1;j<=cnt&&i*pri[j]<=1000000;j++)
        {
            b[i*pri[j]]=1;
            if(i%pri[j]==0)
                break;
            miu[i*pri[j]]=-miu[i];
        }
    }
    ll ans=0;
    for(int i=1;(ll)i*i<=n;i++)
        if(miu[i])
        {
            ll s=0;
            for(int j=1;;j++)
            {
                int z=min((n/i/i-(ll)j*j)/j,(ll)j);
                if(z<=0)
                    break;
                s+=z;
            }
            ans+=miu[i]*s;
        }
    ans=ans*2-1;
    printf("%lld\n",ans);
//  int ans=0;
//  for(int i=1;i<=n;i++)
//      for(int j=1;j<=n;j++)
//          if(gcd(i,j)==1&&i*(i+j)<=n&&j*(i+j)<=n)
//              ans++;
//  printf("%d\n",ans);
//  return 0;
//  int n=1000;
//  for(int i=1;i<=n;i++)
//      for(int j=i;j<=n;j++)
//          if(gcd(i,j)==1&&check(i*(i+j),j*(i+j))&&gcd(gcd(i*(i+j),j*(i+j)),i*j)==1)
//              if(j*(i+j)<=n)
//                  printf("%d %d %d\n",i*(i+j),j*(i+j),i*j);
//  return 0;
//  for(int i=1;i<=n;i++)
//      for(int j=i;j<=n;j++)
//          if(i*j%(i+j)==0)
//          {
//              int z=i*j/(i+j);
//              if(gcd(gcd(i,j),z)==1)
//                  printf("%d %d %d %d %d %d\n",i,j,z,gcd(i,j),i/gcd(i,j),j/gcd(i,j));
//          }
//  return 0;
}