#### Perfect Square Tiling

```% Attempts to perfectly pack squares with the given sides (in decreasing order
%   of magnitude) into a rectangle with the given height and width. The
%   solution is expressed using Bouwkamp notation. In this notation, brackets
%   are used to group adjacent squares with flush tops, and then the groups are
%   sequentially placed in the highest (and leftmost) possible slots.
% e.g. solve([18,15,14,10,9,8,7,4,1], 33, 32, X)
%   gives X=[[18,14],[4,10],[15,7],[1,9],[8]]
% e.g. solve([25,24,23,22,19,17,11,6,5,3], 47, 65, X)
%   gives X=[[25,17,23],[11,6],[5,24],[22,3],[19]]
% e.g. solve([36,33,28,25,16,9,7,5,2], 61, 69, X)
%   gives X=[[36,33],[5,28],[25,9,2],[7],[16]]
% e.g. solve([211,179,167,157,149,143,135,113,100,93,88,87,67,62,50,34,33,27,
%             25,23,22,19,16,15,4], 503, 503, X)
%   gives X=[[211,135,157],[113,22],[179],[149,62],[25,88],[87],[100,167],
%            [143,93],[33,67],[50,27,16],[15,34],[23,4],[19]]
solve(Squares, Height, Width, Groups):-
compatible(Squares, 0, Height, Width),
for(2, 5, MaxLen),
solve_1(Squares, [r(Height,Width)], MaxLen, Groups), !.
% Remove the cut in the previous line if all solutions are required.

solve_1([], _, _, []).
solve_1(Squares, Rectangles, MaxLen, [Group|Groups]):-
highest(Rectangles, r(Height,Width)),
group(Width, Squares, MaxLen, Group0),
efface_all(Group0, Squares, Squares1),
permute(Group0, Group),
split(Rectangles, Height, Group, Rectangles1),
join(Rectangles1, Rectangles2),
solve_1(Squares1, Rectangles2, MaxLen, Groups).

% Check that the sum of the areas of the squares is equal to the area of the
%   rectangle.
% e.g. compatible([18,15,14,10,9,8,7,4,1], 0, 33, 32)
compatible([], Sum, Height, Width):-
Sum is Height * Width, !.
compatible([], _, _, _):-
write(`Incompatible data`), nl,
fail.
compatible([Square|Squares], Sum, Height, Width):-
Sum1 is Sum + Square * Square,
compatible(Squares, Sum1, Height, Width).

% Returns the rectangle with the greatest height.
% e.g. highest([r(15,22),r(9,30)], X) gives X=r(15,22)
highest([Rectangle0|Rectangles0], Rectangle):-
highest_1(Rectangles0, Rectangle0, Rectangle).

highest_1([], Rectangle, Rectangle).
highest_1([r(H1,W1)|Rectangles0], r(H2,_), Rectangle):-
H1 >= H2, !,
highest_1(Rectangles0, r(H1,W1), Rectangle).
highest_1([_|Rectangles0], Rectangle0, Rectangle):-
highest_1(Rectangles0, Rectangle0, Rectangle).

% group(Width, Squares, MaxLen, Group) is true if the sum of the elements of
%   the Group (chosen from the Squares) is equal to the Width, the number of
%   elements in Group is less than or equal to MaxLen, and the elements of the
%   Group are in the same order as they are in the Squares.
%   On backtracking, all solutions are found.
% e.g. findall(X, group(65, [25,24,23,22,19,17,11,6,5,3], 3, X), Xs)
%   gives Xs=[[25,23,17],[24,22,19]]
group(0, _, _, []):-!.
group(Width, [Square|Squares], MaxLen, [Square|Group]):-
MaxLen > 0,
Width >= Square,
Width1 is Width - Square,
MaxLen1 is MaxLen - 1,
group(Width1, Squares, MaxLen1, Group).
group(Width, [_|Squares], MaxLen, Group):-
group(Width, Squares, MaxLen, Group).

% Replaces the first rectangle of the given Height by one rectangle for each
%   element (a square of side Side) of Group, a new rectangle having height
%   Height-Side and width Side.
% e.g. split([r(15,22),r(9,10)], 15, [15,7], X)
%   gives X=[r(0,15),r(8,7),r(9,10)]
split([r(Height,_)|Rectangles0], Height, Group, Rectangles):-!,
split_1(Group, Height, Rectangles1),
append(Rectangles1, Rectangles0, Rectangles).
split([Rectangle|Rectangles0], Height, Group, [Rectangle|Rectangles]):-
split(Rectangles0, Height, Group, Rectangles).

split_1([], _, []).
split_1([Side|Group], Height, [r(Height1,Side)|Rectangles]):-
Height1 is Height - Side,
Height1 >= 0,
split_1(Group, Height, Rectangles).

% Replaces two or more consecutive rectangles having the same height by a
%   single rectangle of the same height, and whose width is the sum of the
%   widths of the consecutive rectangles.
% e.g. join([r(3,1),r(3,2),r(3,3),r(1,1),r(3,4),r(3,5)], X)
%   gives X=[r(3,6),r(1,1),r(3,9)]
join([Rectangle|Rectangles0], Rectangles):-
join_1(Rectangles0, Rectangle, Rectangles).

join_1([], Rectangle, [Rectangle]).
join_1([r(Height,Width)|Rectangles0], r(Height,Width1), Rectangles):-!,
Width2 is Width + Width1,
join_1(Rectangles0, r(Height,Width2), Rectangles).
join_1([Rectangle0|Rectangles0], Rectangle, [Rectangle|Rectangles]):-
join_1(Rectangles0, Rectangle0, Rectangles).

/* append(Xs, Ys, Zs) is true if Zs is the result of appending the list Xs   */
/*   to the list Ys.                                                         */
%append([], Ys, Ys).
%append([X|Xs], Ys, [X|Zs]):-append(Xs, Ys, Zs).

/* efface_all(Xs, Ys, Zs) is true if Zs is the result of deleting the first  */
/*   occurrence of all elements of the list Xs from the list Ys.             */
efface_all([], Ys, Ys).
efface_all([X|Xs], Ys, Zs):-remove(X, Ys, Ws), !, efface_all(Xs, Ws, Zs).

/* for(I, J, K) is true if K is an integer between I and J inclusive.        */
for(I, J, I):-I =< J.
for(I, J, K):-I < J, I1 is I + 1, for(I1, J, K).

/* member(X, Xs) is true if the element X is contained in the list Xs.       */
%member(X, [X|_]).
%member(X, [_|Xs]):-member(X, Xs).

/* permute(Xs, Ys) is true if Ys is a permutation of the list Xs.            */
permute([], []).
permute(Xs, [Y|Ys]):-remove(Y, Xs, Zs), permute(Zs, Ys).

/* remove(X, Ys, Zs) is true if Zs is the result of removing one occurrence  */
/*   of the element X from the list Ys.                                      */
%remove(X, [X|Ys], Ys).
%remove(X, [Y|Ys], [Y|Zs]):-remove(X, Ys, Zs).
```