题目大意：
对于给定的$n(n<2^{31})$，求$displaystylesum_{a=1}^nsum_{b=1}^n[(a+b)|ab]$。
思路：
设$d=gcd(a,b),a=id,b=jd$，则原式=$displaystylesum_{d=1}^nsum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{lfloorfrac n d
floor}[(i+j)|ijd]$。
因为$gcd(i,j)=1$，所以$gcd(i,i+j)=gcd(j,i+j)=1$，所以$(i+j)|ijd$只能是$(i+j)|d$，即原式=$displaystylesum_{d=1}^nsum_{i=1}^{lfloorfrac n d
floor}sum_{j=1}^{lfloorfrac n d
floor}[(i+j)|d]$。
设$k=frac d{(i+j)}$，则$b=jd=jk(i+j)leq n$，即对于每组$i,j$，符合条件的$k$有$lfloorfrac n{(i+j)j}
floor$个。
由于$igeq1,(i+j)jleq n$，所以$jleqsqrt n-1$。令$i<j$，则原式=$displaystylesum_{j=1}^{sqrt n-1}sum_{i=1}^{j-1}[gcd(i,j)=1]lfloorfrac n{(i+j)j}
floor$。
$i$可以数论分块枚举，但不是每一个$i$都对答案有贡献，用$iperp j$表示$gcd(i,j)=1$考虑计算一段$displaystylesum_{i=l}^r[iperp j]$。
显然，$displaystylesum_{i=l}^r[iperp j]=displaystylesum_{i=1}^r[iperp j]-displaystylesum_{i=1}^{l-1}[iperp j]$。而$displaystylesum_{i=1}^n[iperp j]=sum_{i=1}^nsum_{d|gcd(i,j)}mu(d)=sum_{d|j}lfloorfrac nd
floor$。
枚举$j$的约数即可。
细节：
注意判断被$0$除，一旦$lfloorfrac n{(i+j)j}
floor=0$就break掉。
1 #include<cmath> 2 #include<cstdio> 3 #include<cctype> 4 #include<algorithm> 5 typedef long long int64; 6 inline int getint() { 7 register char ch; 8 while(!isdigit(ch=getchar())); 9 register int x=ch^'0'; 10 while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0'); 11 return x; 12 } 13 const int N=46341,M=4793; 14 bool vis[N]; 15 int mu[N],p[M],d[M]; 16 inline void sieve(const int &n) { 17 mu[1]=1; 18 for(register int i=2;i<=n;i++) { 19 if(!vis[i]) { 20 p[++p[0]]=i; 21 mu[i]=-1; 22 } 23 for(register int j=1;j<=p[0]&&i*p[j]<=n;j++) { 24 vis[i*p[j]]=true; 25 if(i%p[j]==0) { 26 mu[i*p[j]]=0; 27 break; 28 } 29 mu[i*p[j]]=-mu[i]; 30 } 31 } 32 } 33 int main() { 34 const int n=getint(); 35 sieve(sqrt(n)); 36 int64 ans=0; 37 for(register int j=1;j<=n/(j+1);j++) { 38 d[0]=0; 39 for(register int i=1;i*i<j;i++) { 40 if(j%i==0) d[++d[0]]=i; 41 } 42 for(register int i=sqrt(j);i;i--) { 43 if(j%i==0) d[++d[0]]=j/i; 44 } 45 for(register int l=1,r;l<j&&n/(l+j)/j;l=r+1) { 46 int64 tmp=0; 47 r=std::min(n/(n/(l+j)/j)/j-j,j-1); 48 for(register int i=1;d[i]<=r;i++) { 49 tmp+=(int64)(r/d[i]-(l-1)/d[i])*mu[d[i]]; 50 } 51 ans+=tmp*(n/j/(l+j)); 52 } 53 } 54 printf("%lld ",ans); 55 return 0; 56 }