设$d=\gcd(x,y)$，且$x=ad, y=bd$, 则$n=abd/(a+b)$。$a,b$都与$a+b$互素，因此$a+b$需要能整除d。
令$d=(a+b)c$，那么$x\le y\le L$就等价于$a\le b$以及$b(a+b)c\le L$，从而解的总数是
$$\sum_{b=1}^{\sqrt{L}}\sum_{\substack{1\le a\le b-1 \\ \gcd(a,b)=1}}\left\lfloor\frac{L}{b(a+b)}\right\rfloor=\sum_{b=1}^{\sqrt{L}}\sum_{\substack{b+1\le a\le 2b-1 \\ \gcd(a,b)=1}}\left\lfloor\frac{L}{ab}\right\rfloor,$$
用莫比乌斯反演处理$\gcd(a,b)=1$的条件，再分段计算求和项($\left\lfloor\frac{L}{x}\right\rfloor$是分段常值函数)即可得结果$5435004633092$。
注：因需使用sympy计算莫比乌斯函数，以下代码为Python 3，代码中还打印了进度信息。
from sympy import sieve
L=10**12
N=10**6
mu=list(sieve.mobiusrange(1,N))
print("initialization completed, start computing...")
def cal(n,start,end):
s=0
end=min(n,end)
x=start
while x <= end and x*x <= n:
s+=n//x
x+=1
if x>end:
return s
x-=1
xPrev=x
x=n//x
while True:
xNext=min(n//x,end)
s+=x*(xNext-xPrev)
xPrev=xNext
if xNext==end:
break
x-=1
return s
tick=10000
res=0
for d in range(1,N):
if mu[d-1]:
a=L//d//d
b,c=1,0
while b*b < a:
c+=cal(a//b,b+1,2*b-1)
b+=1
res+=c*mu[d-1]
if d%tick==0:
print(d//tick,"percent completed, current sum:",res)
print(res)
|
