Infinite loop in prolog? Or just very slow?

Question

I'm trying to figure out if I have an infinite loop in my Prolog program, or if I just did a bad job of writing it, so its slow. I'm trying to solve the square sum chains problem from the dailyprogrammer subreddit. Given a number N, find an ordering of the numbers 1-N (inclusive) such that the sum of each pair of adjacent numbers in the ordering is a perfect square. The smallest N that this holds for is 15, with the ordering [8, 1, 15, 10, 6, 3, 13, 12, 4, 5, 11, 14, 2, 7, 9]. This is the code that I'm trying to use to solve the problem:

is_square(Num):- is_square_help(Num, 0).

is_square_help(Num, S):- Num =:= S * S.
is_square_help(Num, S):- 
    Num > S * S,
    T is S+1,
    is_square_help(Num, T).
is_square_help(Num, S):- Num < S * S, fail.

contains(_, []):- fail.
contains(Needle, [Needle|_]).
contains(Needle, [_|Tail]):- contains(Needle, Tail).

nums(0, []).
nums(Num, List) :- length(List, Num), nums_help(Num, List).

nums_help(0, _).
nums_help(Num, List) :- 
    contains(Num, List),
    X is Num - 1,
    nums_help(X, List).

square_sum(Num, List) :- 
    nums(Num, List),
    square_sum_help(List).

square_sum_help([X, Y|T]) :- 
    Z is X + Y,
    is_square(Z),
    square_sum_help(T).

Currently, when I run square_sum(15, List)., the program does not terminate. I've left it alone for about 10 minutes, and it just keeps running. I know that there are problems that take a long time to solve, but others are reportedly generating answers in the order of milliseconds. What am I doing wrong here?

Answer 1

SWI-Prolog allows this compact implementation

square_sum(N,L) :-
    numlist(1,N,T),
    select(D,T,R),
    adj_squares(R,[D],L).

adj_squares([],L,R) :- reverse(L,R).
adj_squares(T,[S|Ss],L) :-
    select(D,T,R),
    float_fractional_part(sqrt(S+D))=:=0,
    adj_squares(R,[D,S|Ss],L).

that completes really fast for N=15

edit as suggested, building the list in order yields better code:

square_sum(N,L) :-
    numlist(1,N,T),
    select(D,T,R),
    adj_squares(R,D,L).

adj_squares([],L,[L]).
adj_squares(T,S,[S|L]) :-
    select(D,T,R),
    float_fractional_part(sqrt(S+D))=:=0,
    adj_squares(R,D,L).

edit

the code above becomes too slow when N grows. I've changed strategy, and attempt now to find an Hamiltonian path into the graph induced by the binary relation. For N=15 it looks like

(here is the code to generate the Graphviz script:

square_pairs(N,I,J) :-
    between(1,N,I),
    I1 is I+1,
    between(I1,N,J),
    float_fractional_part(sqrt(I+J))=:=0.
square_pairs_graph(N) :-
    format('graph square_pairs_N_~d {~n', [N]),
    forall(square_pairs(N,I,J), format(' ~d -- ~d;~n', [I,J])),
    writeln('}').

)

and here the code for lookup a path

hamiltonian_path(N,P) :-
    square_pairs_struct(N,G),
    between(1,N,S),
    extend_front(1,N,G,[S],P).

extend_front(N,N,_,P,P) :- !.
extend_front(Len,Tot,G,[Node|Ins],P) :-
    arg(Node,G,Arcs),
    member(T,Arcs),
    \+memberchk(T,Ins),
    Len1 is Len+1,
    extend_front(Len1,Tot,G,[T,Node|Ins],P).


struct_N_of_E(N,E,S) :-
    findall(E,between(1,N,_),As),
    S=..[graph|As].

square_pairs_struct(N,G) :-
    struct_N_of_E(N,[],G),
    forall(square_pairs(N,I,J), (edge(G,I,J),edge(G,J,I))).
edge(G,I,J) :-
    arg(I,G,A), B=[J|A], nb_setarg(I,G,B).

Answer 2

Here is a solution using Constraint Logic Programming:

squares_chain(N, Cs) :-
        numlist(1, N, Ns),
        phrase(nums_partners(Ns, []), NPs),
        group_pairs_by_key(NPs, Pairs),
        same_length(Ns, Pairs),
        pairs_values(Pairs, Partners),
        maplist(domain, Is0, Partners),
        circuit([D|Is0]),
        labeling([ff], Is0),
        phrase(chain_(D, [_|Is0]), Cs).

chain_(1, _) --> [].
chain_(Pos0, Ls0) --> [Pos],
        { Pos0 #> 1, Pos #= Pos0 - 1,
          element(Pos0, Ls0, E) },
        chain_(E, Ls0).

plus_one(A, B) :- B #= A + 1.

domain(V, Ls0) :-
        maplist(plus_one, Ls0, Ls),
        foldl(union_, Ls, 1, Domain),
        V in Domain.

union_(N, Dom0, Dom0\/N).

nums_partners([], _) --> [].
nums_partners([N|Rs], Ls) -->
        partners(Ls, N), partners(Rs, N),
        nums_partners(Rs, [N|Ls]).

partners([], _) --> [].
partners([L|Ls], N) -->
        (   { L + N #= _^2 } -> [N-L]
        ;   []
        ),
        partners(Ls, N).

Sample query and answers:

?- squares_chain(15, Cs).
Cs = [9, 7, 2, 14, 11, 5, 4, 12, 13|...] ;
Cs = [8, 1, 15, 10, 6, 3, 13, 12, 4|...] ;
false.

A longer sequence:

?- time(squares_chain(100, Cs)).
15,050,570 inferences, 1.576 CPU in 1.584 seconds (99% CPU, 9549812 Lips)
Cs = [82, 87, 57, 24, 97, 72, 28, 21, 60|...] .

Answer 3

What you are doing wrong is mainly that you generate the whole list before you start testing.

The two clauses that call fail are pointless. Removing them will not change the program. The only reason for doing that is if you do something side-effect-y, like printing output.

Your code for generating the list, and all permutations, seems to work, but it can be done much simpler by using select/3.

You don't seem to have a base case in square_sum_help/1, and you also seem to only check every other pair, which would have lead to problems in some years or whatever when your program had gotten around to checking the correct ordering.

So, by interleaving the generation and testing, like this

square_sum(Num,List) :-
  upto(Num,[],List0),
  select(X,List0,List1),
  square_sum_helper(X,List1,[],List).

square_sum_helper(X1,Rest0,List0,List) :-
  select(X2,Rest0,Rest),
  Z is X1 + X2,
  is_square(Z,0),
  square_sum_helper(X2,Rest,[X1|List0],List).
square_sum_helper(_,[],List0,List) :- reverse(List0,List).

is_square(Num,S) :-
  Sqr is S * S,
  ( Num =:= Sqr ->
    true
  ; Num > Sqr,
    T is S + 1,
    is_square(Num,T) ).

upto(N,List0,List) :-
  ( N > 0 ->
    M is N - 1,
    upto(M,[N|List0],List)
  ; List = List0 ).

the correct result is produced in around 9 msec (SWI Prolog).

?- ( square_sum(15,List), write(List), nl, fail ; true ).
[8,1,15,10,6,3,13,12,4,5,11,14,2,7,9]
[9,7,2,14,11,5,4,12,13,3,6,10,15,1,8]

?- time(square_sum(15,_)).
% 37,449 inferences, 0.009 CPU in 0.009 seconds (100% CPU, 4276412 Lips)

Edit: fixed some typos.

Answer 4

contains/2: clause contains(_, []):- fail. is buggy and redundant at best.

you should type in the body !, fail.

But it's not needed because that what is unprovable shouldn't be mentioned (closed world assumption).

btw contains/2 is in fact member/2 (built-in)