设k=gcd(a,b),a=x*k,b=y*k。
x*k+y*k=(x+y)*k整除x*y*k*k。
因为gcd(x,y)=1,gcd(x+y,x)=gcd(x+y,y)=1,所以x+y整除k。
Ans=∑⌊n√−1⌋y=1∑y−1x=1[gcd(x,y)=1]⌊ny∗(x+y)⌋
对于 ⌊ny∗(x+y)⌋ 的值进行枚举。
被sqrt卡了。。。
#include<iostream>
#include<cstdio>
#include<algorithm>
#define ll long long
#define N 100000
using namespace std;
int n,mu[N],k[N],ss[N],b[N],l,h,t;
ll Ans;
int main()
{
scanf("%d",&n);
mu[1]=1;
for (int i=2;i<=50000;i++)
{
if (ss[i]==0) mu[i]=-1,b[++l]=i;
for (int j=1;j<=l;j++)
{
int t=i*b[j];
if (t>50000) break;
ss[t]=1;
if (i%b[j]==0){mu[t]=0;break;}
mu[t]=-mu[i];
}
}
for (int i=2;i<=n/(i+1);i++)
{
l=0;int j;
for (j=1;j*j<i;j++)
if (i%j==0) k[++l]=j,k[++l]=i/j;
if (j*j==i) k[++l]=j;
sort(k+1,k+l+1);
for (h=1;h<i&&i<=n/(h+i);h=t+1)
{
t=min(n/(n/i/(h+i))/i-i,i-1);
ll sum=0;
for (int j=1;k[j]<=t;j++)
sum+=mu[k[j]]*(t/k[j]-(h-1)/k[j]);
Ans+=n/i/(h+i)*sum;
}
}
printf("%lld\n",Ans);
return 0;
}