[BZOJ2671]Calc

题目大意：
　　对于给定的$n(n<2^{31})$，求$\displaystyle\sum_{a=1}^n\sum_{b=1}^n[(a+b)|ab]$。

思路：
　　设$d=\gcd(a,b),a=id,b=jd$，则原式=$\displaystyle\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor}[(i+j)|ijd]$。
　　因为$\gcd(i,j)=1$，所以$\gcd(i,i+j)=\gcd(j,i+j)=1$，所以$(i+j)|ijd$只能是$(i+j)|d$，即原式=$\displaystyle\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor}[(i+j)|d]$。
　　设$k=\frac d{(i+j)}$，则$b=jd=jk(i+j)\leq n$，即对于每组$i,j$，符合条件的$k$有$\lfloor\frac n{(i+j)j}\rfloor$个。
　　由于$i\geq1,(i+j)j\leq n$，所以$j\leq\sqrt n-1$。令$i<j$，则原式=$\displaystyle\sum_{j=1}^{\sqrt n-1}\sum_{i=1}^{j-1}[\gcd(i,j)=1]\lfloor\frac n{(i+j)j}\rfloor$。
​ $i$可以数论分块枚举，但不是每一个$i$都对答案有贡献，用$i\perp j$表示$\gcd(i,j)=1$考虑计算一段$\displaystyle\sum_{i=l}^r[i\perp j]$。
　　显然，$\displaystyle\sum_{i=l}^r[i\perp j]=\displaystyle\sum_{i=1}^r[i\perp j]-\displaystyle\sum_{i=1}^{l-1}[i\perp j]$。而$\displaystyle\sum_{i=1}^n[i\perp j]=\sum_{i=1}^n\sum_{d|\gcd(i,j)}\mu(d)=\sum_{d|j}\lfloor\frac nd\rfloor$。
　　枚举$j$的约数即可。

细节：
　　注意判断被$0$除，一旦$\lfloor\frac n{(i+j)j}\rfloor=0$就break掉。

#include<cmath>
#include<cstdio>
#include<cctype>
#include<algorithm>
typedef long long int64;
inline int getint() {
    register char ch;
    while(!isdigit(ch=getchar()));
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return x;
}
const int N=46341,M=4793;
bool vis[N];
int mu[N],p[M],d[M];
inline void sieve(const int &n) {
    mu[1]=1;
    for(register int i=2;i<=n;i++) {
        if(!vis[i]) {
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(register int j=1;j<=p[0]&&i*p[j]<=n;j++) {
            vis[i*p[j]]=true;
            if(i%p[j]==0) {
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
}
int main() {
    const int n=getint();
    sieve(sqrt(n));
    int64 ans=0;
    for(register int j=1;j<=n/(j+1);j++) {
        d[0]=0;
        for(register int i=1;i*i<j;i++) {
            if(j%i==0) d[++d[0]]=i;
        }
        for(register int i=sqrt(j);i;i--) {
            if(j%i==0) d[++d[0]]=j/i;
        }
        for(register int l=1,r;l<j&&n/(l+j)/j;l=r+1) {
            int64 tmp=0;
            r=std::min(n/(n/(l+j)/j)/j-j,j-1);
            for(register int i=1;d[i]<=r;i++) {
                tmp+=(int64)(r/d[i]-(l-1)/d[i])*mu[d[i]];
            }
            ans+=tmp*(n/j/(l+j));
        }
    }
    printf("%lld\n",ans);
    return 0;
}

 

​