Prolog: pythagorean triple

Question

I have this code that uses an upper bound variable N that is supposed to terminate for X and Y of the pythagorean triple. However it only freezes when it reaches the upper bound. Wasn't sure how to use the cut to stop the backtracking. Code is:

is_int(0).
is_int(X) :- is_int(Y), X is Y+1.
minus(S,S,0).
minus(S,D1,D2) :- S>0, S1 is S-1, minus(S1,D1,D3), D2 is D3+1.

pythag(X,Y,Z,N) :- int_triple(X,Y,Z,N), Z*Z =:= X*X + Y*Y.

int_triple(X,Y,Z,N) :- is_int(S), minus(S,X,S1), X>0, X<N,
                                  minus(S1,Y,Z), Y>0, Y<N.

Will be called, for example with,

?- pythag(X,Y,Z,20).

Answer 1

First, let us test your solution:

?- pythag(X,Y,Z,20).
X = 4,
Y = 3,
Z = 5 ;
X = 3,
Y = 4,
Z = 5 ;
X = 8,
Y = 6,
Z = 10 ;
X = 6,
Y = 8,
Z = 10 ;
X = 12,
Y = 5,
Z = 13 ;
X = 5,
Y = 12,
Z = 13 ;
X = 12,
Y = 9,
Z = 15 ;
X = 9,
Y = 12,
Z = 15 ;
X = 15,
Y = 8,
Z = 17 ;
X = 8,
Y = 15,
Z = 17 ;
X = 16,
Y = 12,
Z = 20 ;
X = 12,
Y = 16,
Z = 20 ...

Looks perfect to me! All answers are correct solutions! ... up to and including this last solution. After that, your program loops.

Before we try to identify the problem, just hold on for a moment: You must be pretty patient to go through 12 (that is: twelve) answers only to find that loop. Do you think that this method will also work for bigger cases? How many answers are you willing to look at before you give up? Isn't there a simpler way to find out about the problem?

There is one interesting observation here: The answers found have (almost) nothing to do with the looping of the program! That is: By looking at the answers, you get (frequently – as in this case) no clue about the actual cause of the loop! So why not turn off all the answers and concentrate on the relevant part! In fact, we can do this as follows:

?- pythag(X,Y,Z,20), false.
** LOOP **

Now, all answers have been removed due to the goal false. What remains is just the final outcome: either termination, or non-termination, or some error. Nothing else. This should facilitate our observations about termination a bit - no more blinding answers scrolling over the screen. Note that this does not solve the problem in general. After all, how long are we willing to wait? 1s ? 1m?

The actual reason of non-termination can be best understood by looking at a relevant failure slice. That is a fragment of the program whose non-termination implies the non-termination of the whole program. See this answer for more details. Here is the relevant failure slice of your program for query pythag(X,Y,Z,20), false:

pythag(X,Y,Z,N) :-
   int_triple(X,Y,Z,N), false,
   Z*Z =:= X*X + Y*Y.

int_triple(X,Y,Z,N) :-
   is_int(S), false,
   minus(S,X,S1), X>0, X<N,
   minus(S1,Y,Z), Y>0, Y<N.

is_int(0) :- false.
is_int(X) :-
   is_int(Y), false,
   X is Y+1.

Note that there are not many things left of your program. E.g., the actual equation is gone (that's more or less the logic part...). Still, this fragment is relevant. And as long as you do not change something within that fragment, the problem will persist! That is guaranteed for a pure monotonic program as this one...

Here is my preferred solution: It uses length/2 and between/3, two frequently supported predicates of the Prolog prologue.

pythag2(X,Y,Z,N) :-
   length(_, N),
   between(1,N,X),
   between(1,N,Y),
   between(1,N,Z),
   Z*Z =:= X*X + Y*Y.

Answer 2

I was recently as well thinking about a Prolog solution to find Pythagorean triples. I came up with a slightly different code. Assume we have a function:

isqrt(a) = floor(sqrt(a))

It is then enough to enumerate x and y, and to check whether x*x+y*y is the square of some z. Namely to check for:

h = x*x+y*y, z = isqrt(h), z*z = h ?

The function isqrt can be implemented via bisection. For symmetry breaking we can enumerate y after x. Assuming N = 99 the resulting code is:

% between(+Integer, +Integer, -Integer)
between(Lo, Hi, _) :-
   Lo > Hi, !, fail.
between(Lo, _, Lo).
between(Lo, Hi, X) :-
   Lo2 is Lo+1, between(Lo2, Hi, X).

% bisect(+Integer, +Integer, +Integer, -Integer)
bisect(Lo, Hi, X, Y) :-
    Lo+1 < Hi, !,
    M is (Lo+Hi) // 2,
    S is M*M,
    (S > X -> bisect(Lo, M, X, Y);
     S < X -> bisect(M, Hi, X, Y);
     M = Y).
bisect(Lo, _, _, Lo).

% pythago(-List)
pythago(X) :-
   X = [A,B,C],
   between(1, 99, A),
   between(A, 99, B),
   H is A*A+B*B,
   bisect(0, H, H, C),
   C =< 99, H =:= C*C.

There should be 50 such Pythagorean tripples, see also Sloan's A046083:

?- findall(-, pythago(_), L), length(L, N).
N = 52.

One might like to cross check with the following CLP(FD) solution.

:- use_module(library(clpfd)).

% pythago3(-List)
pythago3(X) :-
   X = [A,B,C],
   X ins 1..99,
   A*A+B*B #= C*C,
   A #=< B,
   label(X).

It gives the same number of solutions:

?- findall(-, pythago3(_), L), length(L, N).
N = 50.