由欧拉线定理可知三角形的中心在$(5/3,0)$处。
设三角形三边长分别是$a,b,c$，周长$p=a+b+c$，外接圆半径为$R$，则又由该页面的公式可得
$$9R^2=a^2+b^2+c^2+25\le p^2+25,$$
由此可得R不超过33334。综合上述两个结论可得以下一组方程
$$x_1^2+y_1^2=x_2^2+y_2^2=x_3^2+y_3^2, \quad x_1+x_2+x_3=5, \quad y_1+y_2+y_3=0.$$
枚举$x_1,y_1$，对每一组固定的$x_1,y_1$有
\begin{align*}
x_2^2+y_2^2&=x_1^2+y_1^2, \\
x_3^2+y_3^2&=x_1^2+y_1^2, \\
x_2+x_3&=5-x_1, \\
y_2+y_3&=-y_1, \\
\end{align*}
前两式相减可得
$$(5-x_1)(x_2-x_3)-y_1(y_2-y_3)=0,$$
替换$x_3,y_3$得
$$(5-x_1)(2x_2+x_1-5)-y_1(2y_2+y_1)=0,$$
即
$$2x_2=\frac{y_1(2y_2+y_1)}{5-x_1}+(5-x_1)=\frac{2y_1}{5-x_1}\cdot y_2+\frac{y_1^2+(5-x_1)^2}{5-x_1},$$
若$x_1\neq 5$，则将其代入第一式得到关于$y_2$的二次方程:
$$\left(\frac{2y_1}{5-x_1}\cdot y_2+\frac{y_1^2+(5-x_1)^2}{5-x_1}\right)^2+4y_2^2=4x_1^2+4y_1^2,$$
整理后得
$$\left(\frac{4y_1^2}{(5-x_1)^2}+4\right)\cdot y_2^2+\frac{4y_1\left(y_1^2+(5-x_1)^2\right)}{(5-x_1)^2}\cdot y_2+\frac{\left(y_1^2+(5-x_1)^2\right)^2}{(5-x_1)^2}-4x_1^2-4y_1^2=0,$$
解此二次方程，若有整数解，相应地再判断$x_2$是否为整数，若两者都是整数，则$x_3,y_3$也是整数，从而得一组满足要求的顶点，求出边长后可计算周长。
若$x_1=5$，则或者$y_1=0$或者$y_2=y_3=-y_1/2$，前者$x_2,y_2,x_3,y_3$只能取$\pm 3,\pm4$之间的组合，后者则需要$y_2,y_3$是整数，且$x_2,x_3$等于$\pm\sqrt{x_1^2+3y_1^2/4}$也是整数。
枚举$x_1,y_1$时可以利用对称性只枚举$0\le y_1\le x_1$的情况（根据实际运行的结果，$x_1$只需枚举至$16000$），再交换两者坐标，以及添加$\pm$号组合。
最后，稳妥起见，再验证一次三点各自与$(5,0)$的连线确实垂直与另外两点的连线，并对所得的三角形去重。最终结果是$2816417.1055$（共有155个这样的三角形）。
注：为便于作浮点数除法，以下为Python 3代码。另外，计算量很大，因此代码中打印了进度信息。
注2：上述关于$y_2$的二次方程，可以计算其判别式得到更简洁一些的条件，不过并不能剪枝，从代码的角度而言仍是直接用浮点数计算更容易书写。
import math
f=lambda a1,b1,a2,b2:math.sqrt((a1-a2)*(a1-a2)+(b1-b2)*(b1-b2))
bound=100000
target=16000
tick=target//100
p=0
s=set()
def update(x1,y1,x2,y2,x3,y3):
global p
global t
a,b,c=f(x1,y1,x2,y2),f(x2,y2,x3,y3),f(x3,y3,x1,y1)
k=str(sorted([(x1,y1),(x2,y2),(x3,y3)]))
if a+b>c and b+c>a and a+c>b and a+b+c<=bound and ((x1-5)*(x2-x3)+y1*(y2-y3))==((x2-5)*(x1-x3)+y2*(y1-y3))==((x3-5)*(x1-x2)+y3*(y1-y2))==0 and k not in s:
s.add(k)
p+=a+b+c
for x in range(1,target):
for y in range(x+1):
if 9*(x*x+y*y)<=10**10+25:
for x1,y1 in {(x,y),(-x,y),(x,-y),(-x,-y),(y,x),(y,-x),(-y,x),(-y,-x)}:
r2=x1*x1+y1*y1
if x1!=5:
a=4*y1*y1/((5-x1)*(5-x1))+4
b=(4*y1*y1*y1)/((5-x1)*(5-x1))+4*y1
c=(y1*y1+(5-x1)*(5-x1))*(y1*y1+(5-x1)*(5-x1))/((5-x1)*(5-x1))-4*x1*x1-4*y1*y1
delta=b*b-4*a*c
if delta>=0:
d=math.sqrt(delta)
for u in [(-b-d)/(2*a),(-b+d)/(2*a)]:
if u.is_integer():
y2=int(u)
if abs(2*y1*y2+y1*y1)%abs(5-x1)==0:
v=(2*y1*y2+y1*y1)//(5-x1)+(5-x1)
if abs(v)%2==0:
x2=v//2
x3=5-x1-x2
y3=-y1-y2
update(x1,y1,x2,y2,x3,y3)
else:
if y1==0:
for x2,y2,x3,y3 in [(3,4,-3,-4),(-3,4,3,-4),(3,-4,-3,4),(-3,-4,3,4)]:
update(x1,y1,x2,y2,x3,y3)
else:
d2=x1*x1+3*y1*y1/4
d=math.sqrt(d2)
if abs(y1)%2==0 and d.is_integer():
y2=y3=-y1//2
x2,x3=int(d),-int(d)
update(x1,y1,x2,y2,x3,y3)
if x%tick==0:print(f"{x//tick} percent completed, current value:{p:.4f},{len(s)}")
print(f"{p:.4f},{len(s)},{t}")
|
