不难看出解的基本形式为
$$a=x^2, \quad b=xy, \quad c=y^2$$
随后$(ad,bd,cd)$又构成一组解。此外三边能构成三角形需要$x^2+xy>y^2$，即$y<(1+\sqrt 5)x/2$。
因此本题的主要策略先用上式找出满足$\gcd(a,b,c)=1$的本原解，再枚举无平方因子数d（否则d的平方因子可以吸收到x和y中）直到超过$2.5\times 10^{13}$为止。
无平方因子数的枚举类似第362题等题，最终结果是$41791929448408$。
注：因数据量较大，代码中打印了进度信息。
import math
bound=25*(10**12)
target=5000000
tick=50000
phi=(1+math.sqrt(5))/2
d=[0]*target
n=2
while n < target:
for i in range(n*n,target,n):
d[i]=d[i]+1
n=n+1
while n < target and d[n]>0:
n=n+1
mob=[1]*(target+1)
mob[0]=0
for x in range(2,target):
if d[x]==0:
for i in range(x,target+1,x):
mob[i]*=-1
for i in range(x*x,target+1,x**2):
mob[i] = 0
d={}
for x in range(2,target):
y=x+1
p=int(phi*x)
q=int((math.sqrt(4*bound-3*x*x)-x)/2)
while y < p+1 and x*x+y*y+x*y<=bound:
a=bound/(x*x+y*y+x*y)
b=min(p,q,int((math.sqrt(4*a*bound-3*a*a*x*x)-a*x)/(2*a)))
d[a]=d.get(a,0)+b-y+1
y=b+1
if x%tick==0:
print x/tick, "percent completed"
print sum(sum(mob[i]*(x/(i*i)) for i in range(1,int(math.sqrt(x))+1))*d[x] for x in d)+bound/3
|
