令$x=1+b^2, y=1+c^2$，三角形面积为$A$，则由海伦公式可得$4A^2=xy-1$，从而$b,c$都需要是偶数。
令$b=2p,c=2q$，则有$A^2=4p^2q^2+p^2+q^2$，即
$$A^2-(4p^2+1)q^2=p^2,$$
此方程是Pell方程的推广，其解具有
$$\begin{pmatrix}q_{n+1} \\ A_{n+1}\end{pmatrix}=\begin{pmatrix}8p^2+1 & 4p \\ 4p(4p^2+1) & 8p^2+! \end{pmatrix}\begin{pmatrix}q_n \\ A_n\end{pmatrix},$$
的形式，另外若$(p,q,A)$是一组解，那么$(q,\pm p,A)$也都可以生成解，还需要找出其中满足$0 < A\le n$的解。
最终结果是$2919133642971$。
N=10**10
bound=1078
def sum_areas(p,q,s):
a,b,c=8*p*p+1,4*p,4*p*(4*p*p+1)
q,s=a*q+b*s,c*q+a*s
return 0 if s>N else s+sum_areas(p,q,s)+sum_areas(q,p,s)+sum_areas(q,-p,s)
print sum(sum_areas(p,0,p) for p in range(1,bound))
|
