#### A Sudoku Solver

It is slightly slower, but requires much less memory than Faster Sudoku Solver.

```/*
A Sudoku puzzle consists of a 9x9 matrix divided into nine 3x3 submatrices.
Some cells have already been allocated numbers, and empty cells are
represented by 0.  These empty cells have to be replaced by numbers between
1 and 9 in such a way that each row, column and submatrix contains each of
the numbers 1 to 9 exactly once.

e.g. go([0,0,0, 0,0,0, 0,0,1,
4,0,0, 5,0,0, 0,2,0,
2,0,0, 4,0,0, 6,0,0,

0,0,5, 6,4,0, 0,0,0,
9,0,0, 0,0,0, 0,0,8,
0,0,0, 0,7,1, 3,0,0,

0,0,1, 0,0,6, 0,0,9,
0,6,0, 0,0,3, 0,0,5,
8,0,0, 0,0,0, 0,0,0], X) gives

X = [3,7,9, 8,6,2, 5,4,1,
4,1,6, 5,3,9, 8,2,7,
2,5,8, 4,1,7, 6,9,3,

1,3,5, 6,4,8, 9,7,2,
9,4,7, 3,2,5, 1,6,8,
6,8,2, 9,7,1, 3,5,4,

5,2,1, 7,8,6, 4,3,9,
7,6,4, 1,9,3, 2,8,5,
8,9,3, 2,5,4, 7,1,6]
*/

% Executing this goal will display all solutions to the given problem.
sudoku:-
go([0,0,0,0,0,0,0,0,1,4,0,0,5,0,0,0,2,0,2,0,0,4,0,0,6,0,0,
0,0,5,6,4,0,0,0,0,9,0,0,0,0,0,0,0,8,0,0,0,0,7,1,3,0,0,
0,0,1,0,0,6,0,0,9,0,6,0,0,0,3,0,0,5,8,0,0,0,0,0,0,0,0]).

% All solutions are displayed.
go(Problem):-
init_values(Problem, 0, N, InitVVs),
N =< 3,
empty(N, EmptyVDs),
init_domains(InitVVs, EmptyVDs, InitVDs),
forward_check(InitVDs, InitVVs, FinalVVs),
sort(FinalVVs, SortedVVs),
extract_solution(SortedVVs, Solution),
show(Solution, 1, N),
fail.
go(_).

% All solutions are found on backtracking.
% e.g. go([0,0,2,4,0,0,0,0,0,0,0,0,3,2,0,0], X)
%     gives X=[1,3,2,4,2,4,1,3,4,1,3,2,3,2,4,1]
% e.g. go([0,0,0,0,0,0,0,2,4,0,0,0,0,1,0,3], X)
%     gives X=[3,2,1,4,1,4,3,2,4,3,2,1,2,1,4,3]
%       and X=[1,2,3,4,3,4,1,2,4,3,2,1,2,1,4,3]
go(Problem, Solution):-
init_values(Problem, 0, N, InitVVs),
empty(N, EmptyVDs),
init_domains(InitVVs, EmptyVDs, InitVDs),
forward_check(InitVDs, InitVVs, FinalVVs),
sort(FinalVVs, SortedVVs),
extract_solution(SortedVVs, Solution).

% Gets each given value V and its associated cell X.
% e.g. With Problem=[0,0,2,4,0,0,0,0,0,0,0,0,3,2,0,0]
%   gives [vv(2,2),vv(3,4),vv(12,3),vv(13,2)]
init_values([], X, _, _):-
N is ip(sqrt(sqrt(X))),
X > N * N * N * N, !,            % X is not a fourth power
write(`Incorrect input data length`), nl,
fail.
init_values([], X, N, []):-!,
N is sqrt(sqrt(X)).
init_values([0|Vs], X, N, VVs):-!, % Empty cell
X1 is X + 1,
init_values(Vs, X1, N, VVs).
init_values([V|Vs], X, N, [vv(X,V)|VVs]):-
X1 is X + 1,
init_values(Vs, X1, N, VVs).

% Initialises domains and dependencies considering all cells to be empty.
% The domain of a cell is the list [1,...,N^2] and the dependencies is a list
%   of all cells in the same row, column or submatrix as the cell.
% So, for N=2, the entry in VDs for cell 3 is vd(3,[1,2,3,4],[0,1,2,7,11,15,6])
empty(N, VDs):-
N2 is N * N,
findall(I, for(1, N2, I), Ds),
N4_1 is N2 * N2 - 1,
empty_1(0, N4_1, N2, N, Ds, VDs).

empty_1(I, N4_1, _, _, _, []):-
I > N4_1, !.
empty_1(X, N4_1, N2, N, Ds, [vd(X,Ds,Ys)|VDs]):-
findall(Y, dependency(X, N, N2, N4_1, Y), Ys),
X1 is X + 1,
empty_1(X1, N4_1, N2, N, Ds, VDs).

% Finds a dependency Y for a given cell X, that is, a cell in the same row,
%   same column or same submatrix as X. Finds all solutions on backtracking.
% e.g. findall(Y, dependency(3, 2, 4, 15, Y), Ys) gives Ys=[0,1,2,7,11,15,6]
dependency(X, _, N2, N4_1, Y):-
for(0, N4_1, Y),
Y =\= X,
X // N2 =:= Y // N2.                 % Y is in the same row as X
dependency(X, _, N2, N4_1, Y):-
for(0, N4_1, Y),
Y =\= X,
X mod N2 =:= Y mod N2.               % Y is in the same column as X
dependency(X, N, N2, N4_1, Y):-
for(0, N4_1, Y),
(X // N2) // N =:= (Y // N2) // N,
(X mod N2) // N =:= (Y mod N2) // N, % Y is in the same submatrix as X
X // N2 =\= Y // N2,                 %   but not in the same row as X
X mod N2 =\= Y mod N2.               %   and not in the same column as X

% Initialises the domains according to the given values. There will be an
%   entry vd(X,Ds,Ys) for each as yet unallocated cell X, where Ds is a list of
%   possible values for this cell, and Ys is a list of unallocated
%   dependencies.
% So, for the problem [0,0,2,4,0,0,0,0,0,0,0,0,3,2,0,0], the entry for cell 8
%   would be vd(8,[1,4],[9,10,11,0,4]).
init_domains([], VDs, VDs).
init_domains([vv(X,V)|VVs], VDs0, VDs):-
remove(vd(X,Ds,Ys), VDs0, VDs2),
member(V, Ds), !,
update_domains(Ys, vv(X,V), VDs2, VDs1),
init_domains(VVs, VDs1, VDs).
init_domains([_|_], _, _):-
write(`Conflict in input data`), nl,
fail.

% Updates domains and dependencies following the assignment of value V to the
%   cell X, where Ys are the dependencies of X.
update_domains([], _, VDs, VDs).
update_domains([Y|Ys], vv(X,V), VDs0, VDs):-
remove(vd(Y,Ds,Ys1), VDs0, VDs1), !,
efface(Ys1, X, Ys2), % Remove X from the dependencies of Y
efface(Ds, V, Ds1),  % Remove V from the domain of Y
Ds1=[_|_],           % Fail if a domain becomes empty
update_domains(Ys, vv(X,V), [vd(Y,Ds1,Ys2)|VDs1], VDs).

% Allocates values to the variables, and updates the domains accordingly.
%   Finds all solutions on backtracking.
forward_check([], VVs, VVs).
forward_check([VD|VDs], VVs0, VVs):-
best_and_rest(VDs, VD, vd(X,Ds,Ys), VDs1), % X has the smallest domain Ds
member(V, Ds),                             % We may backtrack here
update_domains(Ys, vv(X,V), VDs1, VDs2),
forward_check(VDs2, [vv(X,V)|VVs0], VVs).

% Finds an element vd(X,Ds,Ys) of [VD|VDs] having the smallest domain Ds, with
%   the most dependencies Ys in the case of ties.
best_and_rest([], Best, Best, []):-!.
best_and_rest(VDs, VD, VD, VDs):-!.
best_and_rest([VD|VDs], VD0, Best, [VD0|Rest]):-
VD=vd(_,D,_),
VD0=vd(_,D0,_),
is_shorter(D, D0), !,
best_and_rest(VDs, VD, Best, Rest).
best_and_rest([VD|VDs], VD0, Best, [VD0|Rest]):-
VD=vd(_,D,Ys),
VD0=vd(_,D0,Ys0),
same_length(D, D0),     % The domains have the same size
is_shorter(Ys0, Ys), !, % Select that with the most dependencies
best_and_rest(VDs, VD, Best, Rest).
best_and_rest([VD|VDs], VD0, Best, [VD|Rest]):-
best_and_rest(VDs, VD0, Best, Rest).

% Extracts the solution.
extract_solution([], []).
extract_solution([vv(_,V)|Vs], [V|Ws]):-
extract_solution(Vs, Ws).

% Displays formatted output
show([], _, _):-!,
nl.
show([X|Xs], I, N):-
I mod (N * N * N) =:= 0, !,
write(X), nl, nl,     % Horizontal gap
I1 = I + 1,
show(Xs, I1, N).
show([X|Xs], I, N):-
I mod (N * N) =:= 0, !,
write(X), nl,         % End of line
I1 = I + 1,
show(Xs, I1, N).
show([X|Xs], I, N):-
I mod N =:= 0, !,
write(X), write(` `), % Vertical gap
I1 = I + 1,
show(Xs, I1, N).
show([X|Xs], I, N):-
write(X),
I1 = I + 1,
show(Xs, I1, N).

% efface(Xs, Y, Zs) is true if Zs is the result of deleting the first
%   occurrence of the element Y from the list Xs.
efface([], _, []). % Succeed if Y is not an element of Xs
efface([Y|Xs], Y, Xs):-!.
efface([X|Xs], Y, [X|Zs]):-efface(Xs, Y, 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).

% is_shorter(Xs, Ys) is true if the list Xs contains fewer elements than
%   the list Ys.
is_shorter([], [_|_]).
is_shorter([_|Xs], [_|Ys]):-is_shorter(Xs, Ys).

% same_length(Xs, Ys) is true if the lists Xs and Ys contain the same
%   number of elements.
same_length([], []).
same_length([_|Xs], [_|Ys]):-same_length(Xs, Ys).

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

% 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).
```