#### Closest pair of a set of points in the plane

```% closest(Points, Point1, Point2, DistSqu) is true if Point1 and Point2 are two
%   points which are nearest together in the list of Points, and DistSqu is the
%   square of the distance between the two points.
% e.g. closest([p(2,3),p(5,4),p(9,6),p(4,7),p(8,1),p(7,2)], P, Q, D)
%   gives P=p(7,2), Q=p(8,1), D=2
closest(Points, Point1, Point2, DistSqu):-
Points=[_,_|_],
merge_sort(Points, 0, SortX),          % Sorted according to the x-coordinate
length(SortX, Length),
closest_1(Length, SortX, [], t(Point1,Point2,DistSqu)).

% Note the similarity with merge_sort_1/5 below.
closest_1(1, [P|Rest], Rest, t(P,P,2147483647)):-!.
closest_1(2, [P,Q|Rest], Rest, t(P,Q,DistSqu)):-!,
distance_squared(P, Q, DistSqu).
closest_1(N, Left, Rest, Closest):-
N1 is N // 2,
N2 is N - N1,
closest_1(N1, Left, Right, ClosestL),  % The closest pair in Left
closest_1(N2, Right, Rest, ClosestR),  % The closest pair in Right
closer(ClosestL, ClosestR, ClosestLR), % The closest pair in Left or Right
Right=[p(Xmedian,_)|_],                % The median of the N points
ClosestLR=t(_,_,DistSqu),
D=sqrt(DistSqu),                       % The distance between closest pair
Xmin is Xmedian - D,
left_strip(N1, Left, Xmin, StripL0),   % Points > Median-D [and < Median]
Xmax is Xmedian + D,
right_strip(N2, Right, Xmax, StripR0), % Points [>= Median and] < Median+D
merge_sort(StripL0, 1, StripL),        % Sorted according to the y-coordinate
merge_sort(StripR0, 1, StripR),        % Sorted according to the y-coordinate
join(StripL, StripR, D, ClosestLR, Closest).

closer(ClosestA, t(_,_,E), ClosestA):-
ClosestA=t(_,_,D),
D < E, !.
closer(_, ClosestB, ClosestB).

% left_strip(N, Points, X, Left) is true if Left contains the elements amongst
%   the first N elements of Points having an x-coordinate greater than X.
left_strip(0, _, _, []):-!.
left_strip(N, [Point|Points], X, [Point|Left]):-
Point=p(X1,_),
X1 > X, !,
N1 is N - 1,
left_strip(N1, Points, X, Left).
left_strip(N, [_|Points], X, Left):-
N1 is N - 1,
left_strip(N1, Points, X, Left).

% right_strip(N, Points, X, Right) is true if Right contains the elements
%   amongst the first N elements of Points having an x-coordinate less than X.
right_strip(0, _, _, []):-!.
right_strip(N, [Point|Points], X, [Point|Right]):-
Point=p(X1,_),
X1 < X, !,
N1 is N - 1,
right_strip(N1, Points, X, Right).
right_strip(N, [_|Points], X, Right):-
N1 is N - 1,
right_strip(N1, Points, X, Right).

% Finds the closest pair, if any, with one point in Ls and the other in Rs.
join([], _, _, Soln, Soln):-!.         % No more Ls; succeed
join(_, [], _, Soln, Soln):-!.         % No more Rs; succeed
join([p(_,Yl)|Ls], Rs, D, Soln0, Soln):-
Rs=[p(_,Yr)|_],
Yl =< Yr - D, !,                     % L is below R-D
join(Ls, Rs, D, Soln0, Soln).        % Try next L with same R
join(Ls, [p(_,Yr)|Rs], D, Soln0, Soln):-
Ls=[p(_,Yl)|_],
Yl >= Yr + D, !,                     % L is above R+D
join(Ls, Rs, D, Soln0, Soln).        % Try same L with next R
join([L|Ls], Rs, D, Soln0, Soln):-     % L is above R-D and below R+D
join_1(Rs, L, D, Soln0, Soln1),      % Find improvements with the current L
join(Ls, Rs, D, Soln1, Soln).        % Try next L with same R

% Finds the closest pair, if any, with one point being L and the other in Rs.
join_1([], _, _, Soln, Soln).                 % No more Rs; succeed
join_1([p(_,Yr)|_], p(_,Yl), D, Soln, Soln):-
Yl >= Yr + D, !.                            % L is above R+D; succeed
join_1([R|Rs], L, D, t(_,_,DistSqu0), Soln):- % L is below R+D
distance_squared(L, R, DistSqu1),
DistSqu1 < DistSqu0, !,                     % L and R form a closer pair
join_1(Rs, L, D, t(L,R,DistSqu1), Soln).
join_1([_|Rs], L, D, Soln0, Soln):-           % L is below R+D
join_1(Rs, L, D, Soln0, Soln).              % L and R are not a closer pair

% merge_sort(List, Axis, SortedList) is true if SortedList is the result of
%   sorting the List of 2-dimensional points along the given Axis.
% e.g. merge_sort([p(3,1),p(4,1),p(5,9),p(2,6),p(5,3)], 0, X) gives
%        X=[p(2,6),p(3,1),p(4,1),p(5,9),p(5,3)]
%      merge_sort([p(3,1),p(4,1),p(5,9),p(2,6),p(5,3)], 1, X) gives
%        X=[p(3,1),p(4,1),p(5,3),p(2,6),p(5,9)]
merge_sort(List, Axis, SortedList):-
length(List, Length),
merge_sort_1(Length, Axis, List, [], SortedList).

merge_sort_1(0, _, Rest, Rest, []):-!.
merge_sort_1(1, _, [A|Rest], Rest, [A]):-!.
merge_sort_1(N, Axis, List, Rest, Sorted):-
N1 is N // 2,
N2 is N - N1,
merge_sort_1(N1, Axis, List, TempList, SortedLeft),
merge_sort_1(N2, Axis, TempList, Rest, SortedRight),
ordered_merge(SortedLeft, SortedRight, Axis, Sorted).

ordered_merge([], Ys, _, Ys).
ordered_merge([X|Xs], [Y|Ys], Axis, [Y|Zs]):-
lt(Axis, Y, X), !,
ordered_merge([X|Xs], Ys, Axis, Zs).
ordered_merge([X|Xs], Ys, Axis, [X|Zs]):-
ordered_merge(Xs, Ys, Axis, Zs).

%  Comparisons on a given Axis of two 2-dimensional points.
lt(0, p(P,_), p(R,_)):-P < R, !.
lt(1, p(_,Q), p(_,S)):-Q < S.

distance_squared(p(X1,Y1), p(X2,Y2), DistanceSquared):-
X1_X2 is X1 - X2,
Y1_Y2 is Y1 - Y2,
DistanceSquared is X1_X2*X1_X2 + Y1_Y2*Y1_Y2.

% naive_closest(Points, Point1, Point2, DistSqu) is true if Point1 and Point2
%   are two points which are nearest together in the list of Points, and
%   DistSqu is the square of the distance between the two points.
% e.g. naive_closest([p(2,3),p(5,4),p(9,6),p(4,7),p(8,1),p(7,2)], P, Q, D)
%   gives P=p(8,1), Q=p(7,2), D=2
naive_closest(Points, Point1, Point2, DistSqu):-
naive_closest_1(Points, t(void,void,2147483647), t(Point1,Point2,DistSqu)).

naive_closest_1([], Closest, Closest).
naive_closest_1([Point|Points], Closest0, Closest):-
naive_closest_2(Points, Point, Closest0, Closest1),
naive_closest_1(Points, Closest1, Closest).

naive_closest_2([], _, Closest, Closest).
naive_closest_2([Point2|Points], Point1, t(_,_,DistSqu0), Closest):-
distance_squared(Point1, Point2, DistSqu1),
DistSqu1 =< DistSqu0, !,
naive_closest_2(Points, Point1, t(Point1,Point2,DistSqu1), Closest).
naive_closest_2([_|Points], Point1, Closest0, Closest):-
naive_closest_2(Points, Point1, Closest0, Closest).

% Generates N random points with X and Y values between Min and Max.
% e.g. rand_points(100, -999, 999, Ps), closest(Ps, P, Q, D)
% e.g. rand_points(100, -999, 999, Ps), naive_closest(Ps, P, Q, D)
rand_points(0, _, _, []):-!.
rand_points(N, Min, Max, [P|Ps]):-
N > 0,
rand_point(Min, Max, P),
N1 is N - 1,
rand_points(N1, Min, Max, Ps).

% Generates a random point with X and Y values between Min and Max.
rand_point(Min, Max, p(X,Y)):-
Min < Max,
rand_int(Min, Max, X),
rand_int(Min, Max, Y).

% rand_int(I, J, K) is true if K is a pseudo-random integer in the range I to J
%   inclusive.
rand_int(I, J, K):-K is int(rand(J - I + 1)) + I.

% length(Xs, L) is true if L is the number of elements in the list Xs.
%length(Xs, L):-length_1(Xs, 0, L).

%length_1([], L, L).
%length_1([_|Xs], L0, L):-L1 is L0 + 1, length_1(Xs, L1, L).
```