function SOLUTIONS=sudoku(PUZZLE)
%Retrieves all solutions to the sudoku given by PUZZLE.
% PUZZLE    is a 9x9 matrix with 0's as the missing entries
% SOLUTIONS is a 1xN cell where SOLUTIONS{i} is the i-th solution.
%
%Example: P=[0 4 0 0 0 5 1 0 0
%            5 6 0 0 0 3 0 0 2
%            0 0 2 6 0 8 0 0 4
%            7 0 5 0 0 0 0 3 0
%            0 0 0 0 0 0 0 0 0
%            0 1 0 0 0 0 5 0 7
%            2 0 0 9 0 6 4 0 0
%            3 0 0 8 0 0 0 2 6
%            0 0 6 4 0 0 0 8 0];
%         S=sudoku(P);
%         then S is a 1x1 cell and
%         S{1}=[8 4 7 2 9 5 1 6 3
%               5 6 1 7 4 3 8 9 2
%               9 3 2 6 1 8 7 5 4
%               7 8 5 1 6 4 2 3 9
%               4 2 3 5 7 9 6 1 8
%               6 1 9 3 8 2 5 4 7
%               2 5 8 9 3 6 4 7 1
%               3 7 4 8 5 1 9 2 6
%               1 9 6 4 2 7 3 8 5];

global SOLUTIONS
SOLUTIONS=cell(0);

jisudoku(PUZZLE)
function jisudoku(PUZZLE)
global SOLUTIONS

A=cell(9); num=zeros(9); %---initialize---
[X,Y]=meshgrid([0 0 0 3 3 3 6 6 6],[0 0 0 1 1 1 2 2 2]);
T=X+Y+1;
I=repmat(find(T==1),[1 9])+repmat([0 3 6 27 30 33 54 57 60],[9 1]);

%---list available elements for free entries,
%   i.e. unique in row/column/field
for c1=1:numel(PUZZLE)
  x=ceil(c1/9);
  y=mod(c1-1,9)+1;
  if ~PUZZLE(y,x)
    A{c1}=1:9;
    A{c1}=setdiff(A{c1},PUZZLE(:,x));
    A{c1}=setdiff(A{c1},PUZZLE(y,:));
    A{c1}=setdiff(A{c1},PUZZLE(I(:,T(y,x))));
  end
  num(c1)=numel(A{c1});
end
  
if any(any(~num&~PUZZLE))%---free entry but no available => quit!  
  return
end
  
[y,x]=find(num==1);      %---exactly one available => put!
if length(x)
  for c1=1:length(x)
    PUZZLE(y(c1),x(c1))=A{y(c1),x(c1)};
  end
  jisudoku(PUZZLE);
  return
end

cnt=0;
for c1=1:9               %---forall 3x3 fields
  u=[];
  for c2=I(:,c1)'        %----collect availables 
    u=[u A{c2}];
  end
  [u,U]=jiunique(u);  
  if length(intersect(1:9,[u PUZZLE(I(:,c1))']))<9
    %---union in field of available & occupied != 1:9 => quit!
    return
  end
  for c3=find(U==1)      %---entry with unique available => put!
    for c2=I(:,c1)'
      if ismember(u(c3),A{c2})
        PUZZLE(c2)=u(c3);
        cnt=cnt+1;        
      end
    end
  end
end

if cnt
  jisudoku(PUZZLE);
  return
end

if all(all(PUZZLE))      %---found solution
  disp('solution')
  SOLUTIONS{numel(SOLUTIONS)+1}=PUZZLE;
else  
  [mv,mp]=min(reshape(num+10*~~PUZZLE,[1 81]));
  disp('fiendish')
  for c1=A{mp}
    PUZZLE(mp)=c1;
    jisudoku(PUZZLE);
    PUZZLE(mp)=0;
  end
end

function [u,U,I,J]=jiunique(x)
  [u,I,J]=unique(x);
  o=ones(size(x));
  U=full(sparse(o,J,o));