Pythagorean triples were already mentioned in Euclid's Elements.
They are numbers x, y and z that form the sides of a right-angled
triangle and that are integers. By the Pythagorean theorem the
numbers x, y and z will satisfy the following Pythagorean
equation, whereby we assume that z is the hypotenuse:
x^{2} + y^{2} = z^{2}
There is a pure Prolog and a CLP(FD) solution. The CLP(FD) solution is straight forward. The pure Prolog solution explicitly enumerates the two legs. For symmetry breaking we consider the second leg greater than the first leg. The code then uses a bisection method to see whether the sum of the squares is itself a square.
Both solutions find the same number of Pythagorean triples in the range 1 to 99:
?- findall(-, pythago(_), L), length(L, N).
N = 50.
?- findall(-, pythago3(_), L), length(L, N).
N = 50.