
% solution(x): P_1 (x) /\ P_2 (x) /\ P_3 (x) /\ P_4 (x)

solution(X) :- isUniverse(2, 62, U), element(X, U),
               p1(X, U), p2(X, U), p3(X, U), p4(X, U).

%solution(X) :- isUniverse(2, 62, U), element(X, U),
%               p1(X, U), p2(X, U), p3(X, U).
%| ?- solution(X).
% 
%X = pair(7,4) ? ;
% 
%X = pair(8,3) ? ;
% 
%X = pair(9,2) ? ;
% 
%X = pair(13,4) ? 
% 
%yes

%solution(X) :- isUniverse(2, 62, U), element(X, U),
%               p1(X, U), p2(X, U).
%| ?- solution(X).
% 
%X = pair(6,5) ? ;
% 
%X = pair(7,4) ? ;
% 
%X = pair(8,3) ? ;
% 
%X = pair(9,2) ? 
% 
%yes



% P_1(x): E y. ((x_1 * x_2 = y_1 * y_2) /\ x =/= y)

p1(X, U) :- element(Y, U),
	    prod(X, P), prod(Y, P),
	    notEqual(X, Y),
            !.


% P_2(x): A_2(x) /\ B_2(x)
% A_2(x): E y. ((x_1 + x_2 = y_1 + y_2) /\ x =/= y)
% B_2(x): A y. (P_1 (y) \/ x_1 + x_2 =/= y_1 + y_2)

p2(X, U) :- a2(X, U), b2(X, U, U).

a2(X, U) :- element(Y, U),
	    sum(X, S), sum(Y, S),
	    notEqual(X, Y),
            !.

b2(_, [], _).
b2(X, [Y | L], U) :- p1(Y, U), !, b2(X, L, U).
b2(X, [Y | _], _) :- sum(X, S), sum(Y, S), !, fail.
b2(X, [_ | L], U) :- b2(X, L, U).


% P_3(x): A y. (x = y \/ ~((x_1 * x_2 = y_1 * y_2) /\ P_2 (y)))

p3(X, U) :- a3(X, U, U).

a3(_, [], _).
a3(X, [X | L], U) :- !, a3(X, L, U).
a3(X, [Y | _], U) :- prod(X, S), prod(Y, S), p2(Y, U), !, fail.
a3(X, [_ | L], U) :- a3(X, L, U).


% P_4(x): A y. (x = y \/ ~((x_1 + x_2 = y_1 + y_2) /\ P_3 (y)))

p4(X, U) :- a4(X, U, U).

a4(_, [], _).
a4(X, [X | L], U) :- !, a4(X, L, U).
a4(X, [Y | _], U) :- sum(X, S), sum(Y, S), p3(Y, U), !, fail.
a4(X, [_ | L], U) :- a4(X, L, U).


% isUniverse(l, h, u): `u' is a lexicographically sorted list of all pairs
%                       of form `pair (x, y)' with l <= x, y <= h, x >= y.

isUniverse(Lo, Hi, U) :- makePairList(Lo, Hi, Lo, Lo, U).

makePairList( _, Hi, Hi, Hi, [pair(Hi, Hi)])     :- !.
makePairList(Lo, Hi,  I,  I, [pair( I,  I) | L]) :-
      !, IP is I + 1, makePairList(Lo, Hi, IP, Lo, L).
makePairList(Lo, Hi,  I,  J, [pair( I,  J) | L]) :-
      JP is J + 1, makePairList(Lo, Hi, I, JP, L).


% Sums and products of pairs.

sum(pair(X1, X2), S)  :- S is X1 + X2.

prod(pair(X1, X2), S) :- S is X1 * X2.


% notEqual(x, y): x =/= y.

notEqual(X, X) :- !, fail.
notEqual(_, _).


% Element of a list.

element(X, [X | _]).
element(X, [_ | L]) :- element(X, L).


