It does if you reduce the domain to positive integers, taken advantage of the problem being quadratic on all vars: s = Solve[a^2 + b^2 - c^2 == 0 && 0 <= a < 100 && 0 <= b < 100 && 0 <= c < 100, {a, b, c}, Integers] ListPointPlot3D[{a, b, c} /. s] ![Mathematica graphics](https://i.sstatic.net/vhseD.png) **Edit** And you can verify that 2 Sum[SquaresR[2, i^2], {i, 0, 99}] -1 == (2 Length@# -Count[#, 0]) &@Flatten[{a, b, c} /. s] (* True *)