Bisschen Brute-Force und ohne Pythagoras, rechtwinklige Dreiecke und Co.
\((x-y)(x+y)=2^2\cdot (2\cdot 5^3) \)
\((x-y)(x+y)=2\cdot (2^2\cdot 5^3)\)
\((x-y)(x+y)=(2\cdot 5^2)\cdot (2^2\cdot 5)\)
... usw. alle Kombinationen.
Man hat dann jeweils \(x-y=a\) und \(x+y=b\) und berechnet dann:$$\frac{1}{2}\begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}\begin{pmatrix} a\\b \end{pmatrix}$$ So gilt z. B. für die erste Gleichung \((x-y)(x+y)=2^2\cdot (2\cdot 5^3)\Rightarrow a=2^2=4 \, \land \, b=2\cdot 5^3=250\):$$\frac{1}{2}\begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}\begin{pmatrix} 4\\250 \end{pmatrix}=\begin{pmatrix} 127\\123 \end{pmatrix}$$ Natürlich gehört zu jedem Lösungspaar \((u,v)\) auch \((-u,v),(u,-v),(-u,-v)\)