sodoku(M,S)

function S = sodoku(M,S)
%S = sodoku(M,[S])
%
%A recursive program that solves 'sodoku' puzzles.
%
%Inputs:  M  partially filled 9x9 matrix with zeros in 'blank' cells
%         S  list of solutions (only used during recursive calls)
%
%Outputs: S  list of solutions
%
%Example:
%
%M = [0,0,1,9,0,0,0,0,8;6,0,0,0,8,5,0,3,0;0,0,7,0,6,0,1,0,0;...
%     0,3,4,0,9,0,0,0,0;0,0,0,5,0,4,0,0,0;0,0,0,0,1,0,4,2,0;...
%     0,0,5,0,7,0,9,0,0;0,1,0,8,4,0,0,0,7;7,0,0,0,0,9,2,0,0];
%
%S = sodoku(M)
%
%Written by G.M. Boynton, 6/3/05
%Optimized by Leif Warland 8/10/05 

%If this is the first call, then zero out the solution matrix
if ~exist('S')
    S = zeros([size(M),0]);
end

%find the blank cells, or zeroes
[I J] = find(~M);
if isempty(I)  %If there aren't any zeros, then we have a solution!
    S(:,:,size(S,3)+1) = M;  %save it
else %otherwise, find the cell with the least number of possible solutions
    minsize = 9;
    for k = 1:length(I)
        newsize = 9;
        ii = (ceil(I(k)/3)-1)*3+1;
        jj = (ceil(J(k)/3)-1)*3+1;
        mm = M(ii:ii+2,jj:jj+2); %these are the indices of the 3x3 block containing that cell    
        for r = 1:9
            if any(mm(:)==r);           newsize = newsize -1;
            elseif   any(M(I(k),:)==r); newsize = newsize -1;
            elseif any(M(:,J(k))==r);   newsize = newsize -1;
            end
        end
        if newsize<minsize
            i = I(k);j = J(k);
            minsize = newsize;
            MM = mm;
        end
    end
    for k=1:9  %loop through all 9 possibilities
        if ~any(M(i,:)==k) & ~any(M(:,j)==k) & ~any(MM(:)==k)  %OK for column, row, and 3x3 block
            M(i,j) = k;  %put this number in,
            S = sodoku(M,S); %and call this function recursively!
        end
    end
end