[s,ef,difflev]=sudokusolve(s,ms,ns)

function [s,ef,difflev]=sudokusolve(s,ms,ns)

%s0=s;

ef=0;
as=[];
if nargin<3,
    ns=Inf;
end
if nargin<2,
    ms=1;
end

step=0;
while step<ns && any(s(:)==0),
    s0=s;
    step=step+1;
    
    %cheking rules
    
    r=cell(9);
    c=cell(9);
    f=cell(9);
    rn=zeros(9);
    cn=zeros(9);
    fn=zeros(9);
    
    ix=find(s(:)==0);
    
    p=cell(1,length(ix));
    pnr=zeros(size(ix));
    
    for i=1:length(ix),
        [row_nr,col_nr]=ind2sub(size(s),ix(i));
        field_nr=(ceil(col_nr/3)-1)*3+ceil(row_nr/3);
        
        ix_row=sub2ind(size(s),row_nr*ones(1,9),1:9);
        ix_col=sub2ind(size(s),1:9,col_nr*ones(1,9));
        xv=ceil(row_nr/3)*3-2:ceil(row_nr/3)*3;
        yv=ceil(col_nr/3)*3-2:ceil(col_nr/3)*3;
        [x,y]=meshgrid(xv,yv);
        ix_field=sub2ind(size(s),x(:)',y(:)');
        
        for nr=1:9,
            st=s;
            st(ix(i))=nr;
            if rulesok(st,row_nr,col_nr),
                r{row_nr,nr}=[r{row_nr,nr} ix(i)];
                rn(row_nr,nr)=rn(row_nr,nr)+1;
                c{col_nr,nr}=[c{col_nr,nr} ix(i)];
                cn(col_nr,nr)=cn(col_nr,nr)+1;
                f{field_nr,nr}=[f{field_nr,nr} ix(i)];
                fn(field_nr,nr)=fn(field_nr,nr)+1;
                p{i}=[p{i} nr];
                pnr(i)=pnr(i)+1;
            end
        end
    end
    
    ix1=find(fn==1);
    for k=1:length(ix1),
        [i,j]=ind2sub([9 9],ix1(k));
        s(f{i,j})=j;
    end
    ix1=find(rn==1);
    for k=1:length(ix1),
        [i,j]=ind2sub([9 9],ix1(k));
        s(r{i,j})=j;
    end
    ix1=find(cn==1);
    for k=1:length(ix1),
        [i,j]=ind2sub([9 9],ix1(k));
        s(c{i,j})=j;
    end 
    
    ix2=find(pnr>=2);
    for i=1:length(ix2),
        
        [row_nr,col_nr]=ind2sub(size(s),ix(ix2(i)));
        field_nr=(ceil(col_nr/3)-1)*3+ceil(row_nr/3);
        
        fnr_r=(find(~ismember(1:3,ceil(field_nr/3)))-1)*3+mod(field_nr-1,3)+1;
        ix_fr1= floor((fnr_r(1)-1)/3)*3*9+row_nr + [0:9:18];
        ix_fr2= floor((fnr_r(1)-1)/3)*3*9+row_nr + [0:9:18];
        
        fnr_c=[ceil(field_nr/3)*3-2:field_nr-1 field_nr+1:ceil(field_nr/3)*3];
        ix_fc1=sub2ind(size(s),mod(fnr_c(1)-1,3)*3+[1:3],ones(1,3)*col_nr);
        ix_fc2=sub2ind(size(s),mod(fnr_c(2)-1,3)*3+[1:3],ones(1,3)*col_nr);
        
        rnr=[ceil(row_nr/3)*3-2:row_nr-1 row_nr+1:ceil(row_nr/3)*3];
        ix_r1=sub2ind(size(s),rnr(1)*ones(1,3),[-2:0]+ceil(field_nr/3)*3);
        ix_r2=sub2ind(size(s),rnr(2)*ones(1,3),[-2:0]+ceil(field_nr/3)*3);
        
        cnr=[ceil(col_nr/3)*3-2:col_nr-1 col_nr+1:ceil(col_nr/3)*3];
        ix_c1=sub2ind(size(s),mod(field_nr-1,3)*3+[1:3],cnr(1)*ones(1,3));
        ix_c2=sub2ind(size(s),mod(field_nr-1,3)*3+[1:3],cnr(2)*ones(1,3));
        
        remove_flag=zeros(size(p{ix2(i)}));
        for j=1:length(p{ix2(i)}),
            nr=p{ix2(i)}(j);
            
            if (~isempty(f{fnr_r(1),nr}) && all(ismember(f{fnr_r(1),nr},ix_fr1))) || (~isempty(f{fnr_r(2),nr}) && all(ismember(f{fnr_r(2),nr},ix_fr2))) ...
                    || (~isempty(f{fnr_c(1),nr}) && all(ismember(f{fnr_c(1),nr},ix_fc1))) || (~isempty(f{fnr_c(2),nr}) && all(ismember(f{fnr_c(2),nr},ix_fc2))) ...
                    || (~isempty(r{rnr(1),nr}) && all(ismember(r{rnr(1),nr},ix_r1))) || (~isempty(r{rnr(2),nr}) && all(ismember(r{rnr(2),nr},ix_r2))) ...
                    || (~isempty(c{cnr(1),nr}) && all(ismember(c{cnr(1),nr},ix_c1))) || (~isempty(c{cnr(2),nr}) && all(ismember(c{cnr(2),nr},ix_c2))),
                remove_flag(j)=1;
            end
        end
        
        ixkeep=find(~remove_flag);
        p{ix2(i)}=p{ix2(i)}(ixkeep);
        pnr(ix2(i))=length(ixkeep);
    end
                    
    ix1=find(pnr==1);

    for i=1:length(ix1),
        s(ix(ix1(i)))=p{ix1(i)};
    end
    
    if isequal(s,s0),
        if any(find(pnr==0)),
            if ~isempty(as),
                s=as(end).s;
                s(as(end).ix)=as(end).nrs(1);
                if length(as(end).nrs)==1,
                    as=as(1:end-1);
                else
                    as(end).nrs=as(end).nrs(2:end);
                end
            else
                ef=1;
                return
            end
        elseif ms
            [m,ixm]=min(pnr);
            as(end+1).s=s;
            as(end).ix=ix(ixm);
            rix=randperm(length(p{ixm}));
            nrs=p{ixm}(rix);
            as(end).nrs=nrs(2:end);
            s(ix(ixm))=nrs(1);
        else
            ef=2;
            return
        end
    end
    
    %disp('');
    
end % while
      
                             
                
%--------------------------------------------------------------------------
function rok=rulesok(s,i,j)

ix=sub2ind(size(s),i,j);
ix1=sub2ind(size(s),i*ones(1,9),1:9);
ix1=ix1(find(ix1~=ix));
ix2=sub2ind(size(s),1:9,j*ones(1,9));
ix2=ix2(find(ix2~=ix));

if i<=3,
    xv=1:3;
elseif i<=6,
    xv=4:6;
else
    xv=7:9;
end

if j<=3,
    yv=1:3;
elseif j<=6,
    yv=4:6;
else
    yv=7:9;
end

[x,y]=meshgrid(xv,yv);

ix3=sub2ind(size(s),x(:)',y(:)');
ix3=ix3(find(ix3~=ix));

rok = ~(any(s(ix)==s(ix1)) | any(s(ix)==s(ix2)) | any(s(ix)==s(ix3)));