function [B,level,err]=yass(A,N,block)
% YASS Yet Another Sudoku Solver
%
% Usage: [B,level,err]=yass(A,N,block)
% Inputs:
% A: a puzzle (9 x 9 array).
% N: the maximum number of solutions to seek, default N=1.
% block: block definition for non-standard Sudoku.
% Outputs:
% B: the solution array. For multiple solutions, B is a cell array.
% level: the level of puzzle: 1: easy, 2: normal, 3: hard.
% err: true if the puzzle is not valid.
%
% Example:
%{
% Example 1: A level-1 puzzle:
A=[0 4 0 0 0 0 3 0 0
0 0 2 0 0 0 0 5 0
0 0 0 0 8 1 0 0 0
0 0 0 0 0 0 4 1 8
0 0 0 5 0 9 0 0 0
0 0 0 0 0 0 0 0 6
0 0 0 4 3 0 2 0 0
1 0 0 0 0 0 0 0 0
0 0 0 7 0 0 0 0 0];
[B,level]=yass(A)
%B =
% 6 4 1 9 5 7 3 8 2
% 9 8 2 3 6 4 1 5 7
% 3 7 5 2 8 1 9 6 4
% 7 5 9 6 2 3 4 1 8
% 8 1 6 5 4 9 7 2 3
% 4 2 3 1 7 8 5 9 6
5 5 9 8 4 3 6 2 7 1
% 1 3 7 8 9 2 6 4 5
% 2 6 4 7 1 5 8 3 9
%level =
% 1
% Example 2: find 100 solutions of an empty puzzle
tic,B=yass(zeros(9),100);toc
% Elapsed time is 0.866873 seconds.
% Example 3: a non-standard sudoku
block = [1 1 1 2 2 2 2 3 3
1 1 1 2 2 3 3 3 3
1 1 4 5 2 2 6 3 3
1 4 4 5 5 2 6 6 3
4 4 4 5 5 5 6 6 6
7 4 4 5 5 8 6 6 9
7 7 4 5 8 8 6 9 9
7 7 7 7 8 8 8 9 9
7 7 8 8 8 9 9 9 9];
[B,level]=yass(A,1,block)
%B =
% 8 4 1 0 5 0 3 0 0
% 0 0 2 0 4 0 1 5 0
% 0 0 0 3 8 1 0 0 0
% 0 0 0 3 0 0 4 1 8
% 0 1 0 5 0 9 0 0 3
% 0 0 0 0 0 0 0 9 6
% 0 0 8 4 3 0 2 7 1
% 1 0 0 0 9 0 0 0 0
% 0 0 0 7 1 0 9 0 0
% level =
% 1
%}
% By Yi Cao at Cranfield University on 25 Feb 2008.
%
% check inputs and outputs
error(nargchk(1,3,nargin));
error(nargoutchk(0,3,nargout));
% default setting
if nargin<2
N=1;
end
if nargin<3,
h = reshape(1:9,3,3);
g = ceil((1:9)/3);
block=h(g,g);
end
% Initialization
[I,J]=find(ones(9));
K=block(:);
B=A;
C=setC(B,I,J,K);
% Solver without recursion
[B,C,err,level]=nonrecsolver(B,C,I,J,K);
if ~any(~B(:)) || err, B=reshape(B,9,[]); return; end
% recursive solver
level=3;
B1=cell(N,1);
[B,C,err,N0,B1]=branch(B,C,I,J,K,1,B1);
B=B1(1:N0-1);
if N==1
B=reshape(B{1},9,[]);
end
% Basic non-recursive solver
function [B,C,err,level]=nonrecsolver(B,C,I,J,K)
% solver through finding the sole option to a grid
level=1;
[B,C,err]=findone(B,C,I,J,K);
if ~any(~B(:)) || err, return; end
% solver through finding grid pairs
level=2;
[B,C,err]=findtwo(B,C,I,J,K);
% Set grid options according to current board
function C=setC(B,I,J,K)
zB=~B;
C=false(81,9);
nc=find(~zB);
C(nc+(B(nc)-1)*81)=true;
C(zB,:)=true;
nB = nnz(B);
idxB = find(B);
for k=1:nB
[C,err]=setOneC(C,idxB(k),I,J,K);
if err
return
end
end
% update board according to grid option map C, then update C afterwards
function [B,C,err]=updateBC(B,C,I,J,K)
zB=~B;
newOne=sum(C(zB(:),:),2)==1;
err=false;
while sum(newOne)>0
zBidx=find(zB);
newOne=zBidx(newOne);
[i,j]=find(C(newOne,:));
B(newOne(i))=j;
for k=newOne(:)'
[C,err]=setOneC(C,k,I,J,K);
if err
return
end
end
zB=~B;
newOne=sum(C(zB(:),:),2)==1;
end
% set map C at grid z
function [C,err]=setOneC(C,z,I,J,K)
zz = I==I(z) | J==J(z) | K==K(z);
zz(z) = false;
cc=C(z,:);
C(zz,cc)=false;
err=any(~any(C(zz,:),2));
% Solver through finding a grid which has a sole option
function [B,C,err]=findone(B,C,I,J,K)
flag=true;
while flag
[B,C,err]=updateBC(B,C,I,J,K);
if err || ~nnz(~B)
return
end
[B,C,err]=findoneC(B,C,I,J,K);
if err || ~nnz(~B)
return
end
C0=C;
C=crosscheck(C,~B,I,J,K);
flag=~isequal(C,C0);
end
% Sub-solver of finding sole option
function [B,C,err]=findoneC(B,C,I,J,K)
err=false;
zB=~B;
flag=1;
while any(zB(:)) && flag,
flag=0;
zC=reshape(find(zB),1,[]);
for z=zC,
is=I==I(z); is(z)=false;
js=J==J(z); js(z)=false;
bs=K==K(z); bs(z)=false;
cc = C(z,:) & (~any(C(is,:)) | ~any(C(js,:)) | ~any(C(bs,:)));
if sum(cc)==1,
flag=1;
ks = (is | js | bs) & zB(:);
C(z,:)=cc;
C(ks,cc)=false;
if any(~any(C(ks,:),2))
err=true;
return
end
[B,C,err]=updateBC(B,C,I,J,K);
if err
return
end
zB(z)=false;
break;
end
end
end
% Solver find grid pair
function [B,C,err]=findtwo(B,C,I,J,K)
C0=C;
zB=~B;
[C,err]=find2sub(C,I,zB);
if err
return
end
[C,err]=find2sub(C,J,zB);
if err
return
end
[C,err]=find2sub(C,K,zB);
if err
return
end
if isequal(C0,C)
return,
end;
[B,C,err]=findone(B,C,I,J,K);
% Sub-solver of finding pair
function [C,err]=find2sub(C,I,zB)
err=false;
for c=1:9,
kc=(I==c)&zB(:);
zc=find(kc)';
nz=length(zc);
if nz<4,
continue,
end
for i=1:nz-1,
ic=zc(i);
kc(ic)=false;
for j=i+1:nz,
jc=zc(j);
kc(jc)=false;
cc0=C(ic,:) | C(jc,:);
cc=cc0 & ~any(C(kc,:));
if sum(cc)==2
C(ic,~cc)=false;
C(jc,~cc)=false;
cc0=cc;
end
if sum(cc0)==2
C(kc,cc0)=false;
if any(~sum(C(kc,:),2))
err=true;
return
end
end
kc(jc)=true;
end
kc(ic)=true;
end
end
% Recursive solver
function [B,C,err,N,B1]=branch(B,C,I,J,K,N,B1)
C0=C;
B0=B;
sumC=sum(C,2);
sumC(B>0)=Inf;
[nz,wz]=min(sumC);
xz=find(C(wz,:));
for k=1:nz
B(wz)=xz(k);
C(wz,xz([1:k-1 k+1:end]))=false;
[C,err]=setOneC(C,wz,I,J,K);
if ~err
[B,C,err]=nonrecsolver(B,C,I,J,K);
end
if ~err
if ~any(~B(:))
B1{N}=B;
N=N+1;
if N>numel(B1)
return
end
else
[B,C,err,N,B1]=branch(B,C,I,J,K,N,B1);
if N>numel(B1)
return
end
end
end
B=B0;
C=C0;
end
err=true;
function [C,err]=crosscheck(C,zB,I,J,K)
err=false;
N=1:9;
I(~zB)=0;
J(~zB)=0;
K(~zB)=0;
for k=1:9
row=false(1,9);
col=false(1,9);
zK = K == k;
Ik = I(zK);
Jk = J(zK);
row(Ik)=true;
col(Jk)=true;
for r=N(row)
rI = I == r;
zI = zK & rI;
kI = zK & ~rI;
pI = rI & ~zK;
z = any(C(zI,:),1) & ~any(C(kI,:),1);
if sum(z)>0
C(pI,z)=false;
if any(~any(C(pI,:),2))
err=true;
return
end
end
z = any(C(zI,:),1) & ~any(C(pI,:),1);
if sum(z)>0
C(kI,z)=false;
if any(~any(C(kI,:),2))
err=true;
return
end
end
end
for c=N(col)
cJ = J == c;
zJ = zK & cJ;
kJ = zK & ~cJ;
pJ = cJ & ~zK;
z = any(C(zJ,:),1) & ~any(C(kJ,:),1);
if sum(z)>0
C(pJ,z)=false;
if any(~any(C(pJ,:),2))
err=true;
return
end
end
z = any(C(zJ,:),1) & ~any(C(pJ,:),1);
if sum(z)>0
C(kJ,z)=false;
if any(~any(C(kJ,:),2))
err=true;
return
end
end
end
end
