function A=sudoku(M)
% SUDOKU - rapidly find all possible solutions to a sudoku puzzle
%
% Usage: Mout = sudoku(M)
%
% M = intial sudoku matrix, with zeros for empty entries
% Mout = solution as a 9x9 matrix if there is a unique solution, or
% as a 9x9xN matrix if there are N solutions
%
% Notes: (1) The algorithm employs recursion, but does as much as it can
% with straighforward deterministic deduction at each recursion
% level to increase overall speed.
% (2) Supplying this function with an empty or overly sparse
% input matrix can create longer calulation times as the
% function searches for all possible solutions.
% (3) A "No solution" error is generated if the input puzzle has no
% valid solution.
% (4) Tested but no warranty, use at your own risk.
% (5) Michael Kleder, December 2006
%
% Examples:
%
% % find the unique solution to this puzzle in a fraction of a second:
% M = [0 0 8 0 9 0 5 0 0;0 0 1 0 7 0 4 0 0;0 0 4 0 3 0 6 0 0;
% 0 1 0 0 0 6 0 0 7;0 9 0 0 0 3 0 0 0;0 2 0 0 5 0 0 6 0;
% 0 5 0 0 4 0 0 2 0;0 0 0 8 0 0 0 3 0;6 0 0 1 0 0 0 4 0];
% sudoku(M)
%
% % find the 100 possible solutions to this puzzle in few seconds:
% M = [0 0 8 0 9 0 5 0 0;0 0 1 0 7 0 0 0 0;0 0 4 0 3 0 6 0 0;
% 0 1 0 0 0 0 0 0 7;0 9 0 0 0 3 0 0 0;0 2 0 0 5 0 0 6 0;
% 0 5 0 0 4 0 0 2 0;0 0 0 8 0 0 0 3 0;0 0 0 1 0 0 0 4 0];
% tic;M=sudoku(M);toc;size(M,3)
% input checking:
if ndims(M)~=2
error('Input matrix must be two dimensional.')
end
if any((size(M)-[9 9])~=0)
error('Input matrix must have nine rows and nine columns.')
end
if any(any(M~=floor(M))) || any(abs(M(:)-4.5)>4.5)
error('Only integers from zero to nine are permitted as input.')
end
% ----------
% main program:
A=0*M; % list of solutions to puzzle so far
[M,imp,A]=recurse(M,A); %#ok (need this syntax for recursion)
if imp
error('No solution.')
end
A=A(:,:,2:end);
return
% ----------
% recursive guess algorithm:
function [M,imp,A]=recurse(M,A)
%clc;disp(M);pause(.1)
[M,imp]=deduce(M); % perform deterministic deductions
if imp % if impossible, quit
return
end
z=find(~M(:)); % indices of unsolved entries
if isempty(z) % if solved
A(:,:,end+1)=M;
return
end
impall=zeros(1,9);
for v=1:9
Q=M;
Q(z(1))=v; % guess
[Q,impall(v),A]=recurse(Q,A);
end
imp=all(impall);
M=Q;
return
% ----------
% deterministic logic algorithm:
function [M,imp]=deduce(M)
%clc;disp(M);pause(.1)
imp=0; % not impossible yet
% solve what you can by deterministic deduction (no guessing)
Mprev = 10*M; % solution at previous stage
while any(M(:)-Mprev(:)) % iterate until no longer changing
Mprev=M;
N=ones(9,9,9);
% zero out untrue entries for input
[r,c]=find(M);
for n=1:length(r)
N(r(n),c(n),:)=0;
N(r(n),c(n),M(r(n),c(n)))=1;
end
if any(any(sum(N,3)<1))
imp=1; % impossible flag (no solution)
return
end
[r,c]=find(sum(N,3)<2); % solved entries
for n=1:length(r)
if any(any(sum(N,3)<1))
imp=1; %impossible flag (no solution)
return
end
v = find(N(r(n),c(n),:)); % value of the entry
M(r(n),c(n)) = v; % make sure entry is recorded
% clear value out of row
N(:,c(n),v)=0;
% clear value out of column
N(r(n),:,v)=0;
% top row & left column of box:
br = floor((r(n)-.5)/3)*3+1;
bc = floor((c(n)-.5)/3)*3+1;
% clear value out of box:
N(br:br+2,bc:bc+2,v)=0;
% reset value into proper place:
N(r(n),c(n),v)=1;
end
% for each entry, find lone possibilites and record
for r=1:9
for c=1:9
v=find(N(r,c,:));
if length(v)==1
M(r,c)=v;
end
end
end
% for each row, find lone possibilities and record
for r=1:9
for v=1:9
c=find(N(r,:,v));
if length(c)==1
M(r,c)=v;
end
end
end
% for each column, find lone possibilities and record
for c=1:9
for v=1:9
r=find(N(:,c,v));
if length(r)==1
M(r,c)=v;
end
end
end
% for each box, find lone possibilities and record
for r=[1 4 7]
for c=[1 4 7]
for v=1:9
Q=N(r:r+2,c:c+2,v);
[pr,pc]=find(Q);
if length(pr)==1
M(r+pr-1,c+pc-1)=v;
end
end
end
end
end
return
