#Spring 2016 ExpMath Class Project
#Jinyoung Park

Help:=proc()
print(`GenSudoku(n):generates a Sudoku matrix with n blanks,`):
print(` which has a unique solution`):
print(`SolveSudoku(M):gives the solution to M`):
end:

with(combinat):

#######################Procedures for the Sudoku Solver


#FindBlank: detects blank spaces
FindBlank:=proc(M) local i,j,S:

S:=[]:

for i from 1 to nops(M) do
 for j from 1 to nops(M[1]) do
  if M[i][j]=0 then
   S:=[op(S),[i,j]]:
  fi:
 od:
od:

S:

end:

#FindPoss: finds possibilities for given blank
FindPoss:=proc(M,i,j) local Msize,Poss,Rowset,Columnset,Localset,k1,k2:

Msize:=nops(M[1]):

Poss:={seq(k,k=1..Msize)}:

Rowset:={seq(M[i][k],k=1..Msize)}:

Poss:=Poss minus Rowset:

Columnset:={seq(M[k][j],k=1..Msize)}:

Poss:=Poss minus Columnset:

 if i <=3 then
  if j <=3 then
   Localset:={seq(seq(M[k1][k2],k1=1..3),k2=1..3)}:
  elif j <=6 then
   Localset:={seq(seq(M[k1][k2],k1=1..3),k2=4..6)}:
  else
   Localset:={seq(seq(M[k1][k2],k1=1..3),k2=7..9)}:
  fi:
 elif i <=6 then
   if j <=3 then
   Localset:={seq(seq(M[k1][k2],k1=4..6),k2=1..3)}:
  elif j <=6 then
   Localset:={seq(seq(M[k1][k2],k1=4..6),k2=4..6)}:
  else
   Localset:={seq(seq(M[k1][k2],k1=4..6),k2=7..9)}:
  fi:
 else
   if j <=3 then
   Localset:={seq(seq(M[k1][k2],k1=7..9),k2=1..3)}:
  elif j <=6 then
   Localset:={seq(seq(M[k1][k2],k1=7..9),k2=4..6)}:
  else
   Localset:={seq(seq(M[k1][k2],k1=7..9),k2=7..9)}:
  fi:
 fi:

Poss:=Poss minus Localset:  

Poss:

end:


#SizeOfPoss: outputs the size of possibilities for given blank
SizeOfPoss:=proc(M,i,j)
 nops(FindPoss(M,i,j)):
end:



#FillSimpleBlank : fills the blank if it has only one possibility
FillSimpleBlank:=proc(M) local Blank,k1,Solu,k,M1:

Solu:=[]:

Blank:=FindBlank(M):

M1:=M:

for k from 1 to nops(Blank) do
 if nops(FindPoss(M,op(Blank[k])))=0 then
  RETURN(FAIL):
 elif nops(FindPoss(M,op(Blank[k])))=1 then
  Solu:=[op(Solu),[Blank[k],op(FindPoss(M,op(Blank[k])))]]:
 fi: od:

if Solu=[] then
 RETURN(FAIL):
else
  for k1 from 1 to nops(Solu) do
  M1[Solu[k1][1][1]][Solu[k1][1][2]]:=Solu[k1][2]:
 od:
fi:

M1:

end:

#FillSimpleBlankIm: Improved FillSimpleBlank(one at a time)
FillSimpleBlankIm:=proc(M) local Blank,k1,Solu,k,M1:

Blank:=FindBlank(M):

M1:=M:

for k from 1 to nops(Blank) do
 if nops(FindPoss(M1,op(Blank[k])))=0 then
  RETURN(FAIL):
 elif nops(FindPoss(M1,op(Blank[k])))=1 then
  M1[Blank[k][1]][Blank[k][2]]:=op(FindPoss(M1,op(Blank[k]))):
 fi:
od:

M1:

end:

#FillSimpleBlank1:Fills simple blanks as many as possible
FillSimpleBlank1:=proc(M) local M1:

M1:=FillSimpleBlankIm(M):

 if (M1=M) or (M1=FAIL) then
  RETURN(M):
 else
  FillSimpleBlank1(M1):
 fi:

end:


#MinPoss: finds the entry with the minimum possibilities(>=2)
MinPoss:=proc(M) local B,i,min:

B:=FindBlank(M):

for i from 2 to nops(B) do
  min:=B[1]: 
  if SizeOfPoss(M,op(B[i]))<SizeOfPoss(M,op(min)) then
 min:=B[i]:
  fi:
od:

min:

end:


#MinPossIm: Improved MinPoss
MinPossIm:=proc(M) local B,i,min:

B:=FindBlank(M):

min:=B[1]:
  
for i from 2 to nops(B) do
 if SizeOfPoss(M,op(min))=2 then
  RETURN(min):
 else
  if SizeOfPoss(M,op(B[i]))<SizeOfPoss(M,op(min)) then
   min:=B[i]:
  fi: fi:
od:

min:

end:

#CheckPoss:=Checks the possible situation for given entry
CheckPoss:=proc(M,i,j,k) local M1:

M1:=M:

M1[i][j]:=k:

FillSimpleBlank1(M1):

end:

#SolveSudoku(M):gives the solution to M
SolveSudoku:=proc(M) local M0,minposs,minpossset,i0,M1,M2,S:
  if M=FAIL then
   RETURN(FAIL):
  fi:

  M0:=FillSimpleBlank1(M):  

  if FindBlank(M0)=[] then
   RETURN({M0}):
  fi:


  minposs:=MinPossIm(M0):
  minpossset:=FindPoss(M0,minposs[1],minposs[2]):
   if minpossset={} then
    RETURN(FAIL):
   fi:
     S:={}:
     for i0 from 1 to nops(minpossset) do
     M2:=CheckPoss(M0,minposs[1],minposs[2],minpossset[i0]):
     M1:=SolveSudoku(M2):
      if (M1<>FAIL and M1<>{}) then
        if type(M1,set) then
         S:=S union M1:
        else
         S:=S union {M1}:
        fi:
      fi:
     od:

if type(S,set) then
 S:
else
 {S}:
fi:

end:

#Primitive version of SolveSudoku; do not use this
SolveSudoku1:=proc(M) local M0,minposs,minpossset,i0,M1,M2:
  M0:=FillSimpleBlank1(M):  

  if FindBlank(M0)=[] then
   RETURN(M0):
   
  else
  minposs:=MinPossIm(M0):
  minpossset:=FindPoss(M0,minposs[1],minposs[2]):
   if minpossset={} then
    RETURN(FAIL):
   else
     for i0 from 1 to nops(minpossset) do
     M2:=CheckPoss(M0,minposs[1],minposs[2],minpossset[i0]):
     M1:=SolveSudoku1(M2):
      if (M1<>FAIL) then         
        print(M1):
      fi:
     od:
   fi:   
fi:

end:

###############Procedures for the Sudoku Generator

#DetectEntries(M):Detect Entries from 1-9
DetectEntries:=proc(M) local S,i,j,k:

for i from 1 to 9 do
 S[i]:={}:
od:

for k from 1 to 9 do
for i from 1 to 9 do
 for j from 1 to 9 do
  
  if M[i][j]=k then
   S[k]:=S[k] union {[i,j]}:
  fi:
  od:
 od:
od:

[seq(S[i],i=1..9)]:

end:

#PermuteEntries(M,p):Permute entries of M with respect to the given permutation
PermuteEntries:=proc(M,p) local S,i,j,M1:
 S:=DetectEntries(M):
 M1:=M:

for i from 1 to 9 do
 for j in S[i] do
  M1[j[1]][j[2]]:=p[i]:
 od:
od:

M1:

end:

#PermuteRows(M,p):Permute rows1-3,4-6,7-9 of M with respect to the given permutation
PermuteRows:=proc(M,p) local M1,i:
 M1:=M:

for i from 1 to 3 do
 M1[i]:=M[p[i]]:
od:

for i from 4 to 6 do
 M1[i]:=M[p[i-3]+3]:
od:

for i from 7 to 9 do
 M1[i]:=M[p[i-6]+6]:
od:

M1:

end:

#PermuteColumns(M,p):Permute columns1-3,4-6,7-9 of M with respect to the given permutation
PermuteColumns:=proc(M,p) local M1,i,j:
 M1:=M:

for j from 1 to 9 do
for i from 1 to 3 do
 M1[j][i]:=M[j][p[i]]:
od:

for i from 4 to 6 do
 M1[j][i]:=M[j][p[i-3]+3]:
od:

for i from 7 to 9 do
 M1[j][i]:=M[j][p[i-6]+6]:
od:

od:

M1:

end: 

#Permute3Rows(M,p):Permute 3 chunks of rows of M with respect to the given permutation
Permute3Rows:=proc(M,p) local M1,i,j:
 M1:=M:

for i from 0 to 2 do
 for j from 1 to 3 do
 M1[3*i+j]:=M[3*(p[i+1]-1)+j]:
 od:
od:

M1:

end:

#Permute3Columns(M,p):Permute 3 chunks of rows of M with respect to the given permutation
Permute3Columns:=proc(M,p) local M1,i,j,k:
 M1:=M:

for i from 1 to 9 do
 for j from 0 to 2 do
  for k from 1 to 3 do
   M1[i][3*j+k]:=M[i][3*(p[j+1]-1)+k]:
  od:
 od:
od:

M1:

end:

#IsSudoku(M):Check if M satisfies the restrictions for Sudoku Matrix
IsSudoku:=proc(M) local i,S,T,j,k,l:

T:={1,4,7}:

if nops(M)<>9 then
 RETURN(FAIL):
fi:

for i from 1 to nops(M) do
 if nops(M[i])<>9 then
 RETURN(FAIL):
 fi:
od:

for i from 1 to 9 do
 S:={op(M[i])}:
 if S<>{seq(i,i=1..9)} then
 RETURN(FAIL):
 fi:
od:

for i in T do
 for j in T do
  if {seq(seq(M[i+k][j+l],l=0..2),k=0..2)}<>{seq(i,i=1..9)} then
   RETURN(FAIL):
  fi:
 od:
od:

true:

end:


#SudokuMatGen():generates a complete random Sudoku Matrix
SudokuMatGen:=proc() local M,p1,p2,p3,p4,p5,M1:
with(combinat):

M:=[[3, 9, 4, 1, 7, 2, 5, 8, 6], [1, 5, 7, 3, 8, 6, 2, 4, 9], [2, 8, 6, 9, 4, 5, 7, 1, 3], [5, 3, 8, 7, 9, 4, 6, 2, 1], [9, 4, 1, 2, 6, 3, 8, 7, 5], [7, 6, 2, 8, 5, 1, 3, 9, 4], [4, 1, 3, 5, 2, 8, 9, 6, 7], [6, 2, 9, 4, 3, 7, 1, 5, 8], [8, 7, 5, 6, 1, 9, 4, 3, 2]]:


p1:=randperm(9):
 if not type(p1,list) then
  RETURN(FAIL):
 fi:
p2:=randperm(3):
p3:=randperm(3):
p4:=randperm(3):
p5:=randperm(3):

M1:=PermuteEntries(M,p1):
M1:=PermuteRows(M1,p2):
M1:=PermuteColumns(M1,p3):
M1:=Permute3Rows(M1,p4):
M1:=Permute3Columns(M1,p5):

if IsSudoku(M)<>true then
 RETURN(FAIL):
fi:

M1:

end:

#BlankGen(n):generates n blanks
BlankGen:=proc(n) local blank,r1,r2,i:
with(combinat):

blank:={}:

for i from 1 while nops(blank)<n do
 r1:=rand(1..9)():
 r2:=rand(1..9)():
 blank:=blank union {[r1,r2]}:
od:

blank:

end:

#BlankMatGen(M,n):=generates n blanks in the given matrix M
BlankMatGen:=proc(M,n) local M1,B,i:
M1:=M:

B:=BlankGen(n):

for i in B do

M1[i[1]][i[2]]:=0:

od:

M1:

end:

#GenSudoku(n):=generates a Sudoku matrix with n blanks, which has a unique solution
GenSudoku:=proc(n) local M,M1:

if n>55 then
 RETURN(print(`Too many blanks, choose smaller n`)):
fi:

if n>45 then
 print(`n is big. It may take some time.`):
 print(`If it takes too long, stop it and try again.`):
fi:

M:=SudokuMatGen():

if M=FAIL then
 GenSudoku(n):
fi:

M1:=BlankMatGen(M,n):

if nops(SolveSudoku(M1))=1 then
 RETURN(M1):
else
 GenSudoku(n):
fi:

end: