#sudoku by Bobby Griffin, 2/13/06
#edited by Bobby Griffin, 4/24/06 (thanks to Andrew Baxter)

Help:=proc():  
print(`Su(M), poss(M,n,rootN,B), Su2(P,n,rootN,B)`):
print(`descendants(P,n,rootN,B), doubleCell(P,n,rootN,B), remove2(P,i,left,right,n,rootN,B,type)`):
print(`singleCell(P,n,rootN,B), singleVal(P,n,rootN,B), solve1(P,B,i,j,opValue,n,rootN)`):
print(`whichBlock(i,j,rootN), champion(P,n), children(P,n,rootN,B,y,z)`):
end:

#################################################################
#Su(M):  The main procedure called by the user
#
#Input:  M, a list of lists, size n^2 by n^2, of a Sudoku puzzle
#
#Output:  Y, a set of solutions for the given Sudoku puzzle
#         or ERROR, if the size is not a perfect square
#################################################################
Su:=proc(M) local pos,n,rootN,B,wB,i,j,k,X,Y,Xi,Xij:
option remember:

n:=nops(M):
rootN:=sqrt(n):
if rootN <> floor(rootN) then
  error(`Not a valid Sudoku problem`):
fi:

B:=[]:
for i from 1 to n do
  B:=[op(B), []]:
od:

for i from 1 to n do
  for j from 1 to n do
    wB := whichBlock(i,j,rootN):
    B[wB] := [op(B[wB]), [i,j]]:
  od:
od:

pos:=poss(M,n,rootN,B):

X:=Su2(pos,n,rootN,B):
Y:={}:

for i from 1 to nops(X) do
  Xi:=X[i]:
  for j from 1 to n do
    Xi[j] := [op(1..n,Xi[j])]:
    Xij := Xi[j]:
    Xij := [op(1..n,Xij)]:
    Xij := [seq(op(Xij[k]),k=1..n)]:
    Xi[j] := Xij:
  od:
  Y:=Y union {Xi}:
od:

return Y:
   
end:
##########################################
#END Su(M)
##########################################


####################################################################################
#poss(M,n,rootN,B):  Takes the original input of M and returns the possibility 
#  matrix.  The possibility matrix has a set of possible elements at each 
#  location (i,j) 1<=i,j<=n.
#  Also, it contains a total value of how many times a number can occur in 
#  each division
#  e.g. P[i][n+1][j] gives how many times the number 'j' can occur in the 'ith' row
#  n+1 corresponds to rows, n+2 to columns, n+3 to blocks
#
#Input:  M, the original input
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#          e.g. B[i][j] gives the coordinates of the 
#          'jth' element in the 'ith' block
#
#Output:  P, the possibility matrix
####################################################################################
poss:=proc(M,n,rootN,B) local P,i,j,k,Mij:
option remember:

P:=[]:

for i from 1 to n do
  P:= [op(P), []]:
  for j from 1 to n do
    P[i] := [op(P[i]), {seq(k, k=1..n)}]:
  od:
  P[i] := [op(P[i]), [n$n], [n$n], [n$n]]:
od:

for i from 1 to n do
  for j from 1 to n do
    Mij := M[i][j]:
    if Mij <> 0 then
      P:=solve1(P,B,i,j,M[i][j],n,rootN)
    fi:
  od:
od:

return P:

end:
###################################################
#END poss(M,n,rootN,B)
###################################################


##############################################################################
#Su2(P,n,rootN,B):  descdendants' recursive partner, this procedure calls
#singleCell,singleVal, and doubleCell until no change is detected
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#
#Output:  Recursively shrinks the possibility matrix along with descendants
##############################################################################
Su2:=proc(P,n,rootN,B) local P1,change,change1,Pnew,X:
option remember:
change:=true:
P1:=P:

while change do
  change1:=true:
  while change1 do
    X:=singleCell(P1, n, rootN, B):
    change:=X[2]:
    change1:=X[2]:
    P1 := X[1]:

    X:=singleVal(P1, n, rootN, B):
    change:=change or X[2]:
    change1:=change1 or X[2]:
    P1 := X[1]:
  od:

  X:=doubleCell(P1,n,rootN,B):
  if P1 <> X then
    P1:=X:
    change:=true:
  fi:
od:

Pnew := descendants(P1,n,rootN,B):

return Pnew:

end:
#############################################
#END Su2(P,n,rootN,B)
#############################################


####################################################################################
#descendants(P,n,rootN,B):  Su2's recursive partner.
#  Uses champion and children to find valid solutions to a given possibility matrix
#  Then uses Su2 to reduce these possibilities
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#
#Output:  P1, a set containing the reduced (via Su2) forms of the children of 
#           the champion of P
#         or {}, if champion returns FAIL
####################################################################################
descendants:=proc(P,n,rootN,B) local i,j,P1,champ,y,z,CH:
option remember:

champ:=champion(P,n):
if champ = FAIL then
  return {}:
fi:


if champ[3] = n+1 then
  return {P}:
fi:

y:=champ[1]:
z:=champ[2]:

CH:=children(P,n,rootN,B,y,z):


P1:={}:
for i from 1 to nops(CH) do
  P1:=P1 union Su2(CH[i],n,rootN,B):
od:

return P1:

end:
###########################################################
#END descendants(P,n,rootN,B)
###########################################################


#####################################################################################
#doubleCell(P,n,rootN,B): checks if any cells in the same division with two elements 
#  are exactly the same.  If so, it uses remove2 to remove these two elements from 
#  other cells in the same division
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#
#Output:  P1, an updated possibility matrix
#####################################################################################
doubleCell:=proc(P,n,rootN,B) local i,j,left,right,P1,temp,nR,nC,nB:
option remember:

P1:=P:
for i from 1 to n do
  for left from 1 to n-1 do   
    nR:=nops(P1[i][left]):
    nC:=nops(P1[left][i]):
    nB:=nops(P1[B[i][left][1]][B[i][left][2]]):  
    for right from left+1 to n do      
      if nR = 2 then
        if P1[i][left] = P1[i][right] then
          P1:=remove2(P1,i,left,right,n,rootN,B,1):
        fi:
      fi:
      if nC = 2 then
        if P1[left][i] = P1[right][i] then
          P1:=remove2(P1,i,left,right,n,rootN,B,2):
        fi:
      fi:
      if nB = 2 then
        if (P1[B[i][left][1]][B[i][left][2]] = 
                P1[B[i][right][1]][B[i][right][2]]) then
          P1:=remove2(P1,i,left,right,n,rootN,B,3):
        fi:
      fi:
    od:
  od:
od:
return P1:
 
end:
###########################################################################
#END doubleCell(P,n,rootN,B)
###########################################################################


##############################################################################
#remove2(P,i,left,right,n,rootN,B,type):  Called by doubleCell in the event
#  that there are two cells in the same division with only the exact same 
#  two elements.  remove2 then removes these two elements from elsewhere in 
#  the division
#
#Input:  P, the possibility matrix
#        i, the index of where the match occurs (e.g. the 'ith' row)
#        left, the left-most cell that matches
#        right, the right-most cell that matches
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#        type, corresponds to the different divisions
#          1=row, 2=column, 3=block
#
#Output:  P1, the updated possibility matrix
##############################################################################
remove2:=proc(P,i,left,right,n,rootN,B,type) local j,P1,op1,op2,wB,y,z:
option remember:

P1:=P:

if type = 1 then
  op1:=P1[i][left][1]:
  op2:=P1[i][left][2]:
  for j from 1 to n do
    if j<>left and j<>right then
      if member(op1,P1[i][j]) then
        P1[i][j] := P1[i][j] minus {op1}:
        P1[i][n+1][op1] := P1[i][n+1][op1] - 1:
        P1[j][n+2][op1] := P1[j][n+2][op1] - 1:
        wB:=whichBlock(i,j,rootN):
        P1[wB][n+3][op1] := P1[wB][n+3][op1] - 1:
      fi:
      if member(op2,P1[i][j]) then
        P1[i][j] := P1[i][j] minus {op2}:
        P1[i][n+1][op2] := P1[i][n+1][op2] - 1:
        P1[j][n+2][op2] := P1[j][n+2][op2] - 1:
        wB:=whichBlock(i,j,rootN):
        P1[wB][n+3][op2] := P1[wB][n+3][op2] - 1:
      fi:
    fi:
  od:
fi:

if type = 2 then
  op1:=P1[left][i][1]:
  op2:=P1[left][i][2]:
  for j from 1 to n do
    if j<>left and j<>right then
      if member(op1,P1[j][i]) then
        P1[j][i] := P1[j][i] minus {op1}:
        P1[j][n+1][op1] := P1[j][n+1][op1] - 1:
        P1[i][n+2][op1] := P1[i][n+2][op1] - 1:
        wB:=whichBlock(j,i,rootN):
        P1[wB][n+3][op1] := P1[wB][n+3][op1] - 1:
      fi:
      if member(op2,P1[j][i]) then
        P1[j][i] := P1[j][i] minus {op2}:
        P1[j][n+1][op2] := P1[j][n+1][op2] - 1:
        P1[i][n+2][op2] := P1[i][n+2][op2] - 1:
        wB:=whichBlock(j,i,rootN):
        P1[wB][n+3][op2] := P1[wB][n+3][op2] - 1:
      fi:
    fi:
  od:
fi:


if type = 3 then
  y:=B[i][left][1]:
  z:=B[i][left][2]:
  op1:=P1[y][z][1]:
  op2:=P1[y][z][2]:
  for j from 1 to n do
    if j<>left and j<>right then
      y:=B[i][j][1]:
      z:=B[i][j][2]:
      if member(op1,P1[y][z]) then
        P1[y][z] := P1[y][z] minus {op1}:
        P1[y][n+1][op1] := P1[y][n+1][op1] - 1:
        P1[z][n+2][op1] := P1[z][n+2][op1] - 1:
        P1[i][n+3][op1] := P1[i][n+3][op1] - 1:
      fi:
      if member(op2,P1[y][z]) then
        P1[y][z] := P1[y][z] minus {op2}:
        P1[y][n+1][op2] := P1[y][n+1][op2] - 1:
        P1[z][n+2][op2] := P1[z][n+2][op2] - 1:
        P1[i][n+3][op2] := P1[i][n+3][op2] - 1:
      fi:
    fi:
  od:
fi:

return P1:

end:
##################################################################
#END remove2(P,i,left,right,n,rootN,B,type)
##################################################################


################################################################################
#singleCell(P,n,rootN,B):  checks for any cells with only one possible entry
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#
#output:  a list containing:
#           P1, the updated possibility matrix
#           flag, a boolean indicating whether or not a change has occurred
################################################################################
singleCell:=proc(P,n,rootN,B) local i,j,P1,temp,flag:
option remember:

P1:=P:
flag:=false:

for i from 1 to n do
  for j from 1 to n do
    temp := P1[i][j]:
    
    if nops(temp)=1 and P1[i][n+1][op(temp)]>-1 then
      P1:=solve1(P1,B,i,j,op(temp),n,rootN):
      flag:=true:
    fi:
  od:
od:

return [P1,flag]:

end:
###################################################################
#END singleCell(P,n,rootN,B)
###################################################################


###############################################################################
#singleVal(P,n,rootN,B):  checks within every division to see if a number 
#  can occur only once
#  i.e. checks P[i][n+k][j], i=1..n, k=1..3, j=1..n to see if the value is '1'
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#
#Output:  a list containing:
#           P1, the updated possibility matrix
#           flag, a boolean indicating whether or not a change has occurred
###############################################################################
singleVal:=proc(P,n,rootN,B) local i,j,k,l, P1, temp, flag,y,z:
option remember:

P1:=P:
flag:=false:

for i from 1 to n do
  for j from 1 to n do
    for k from 1 to 3 do
      if P1[i][n+k][j] = 1 then
        for l from 1 to n do
          if k = 1 then
            if member(j,P1[i][l]) then
	      P1:=solve1(P1,B,i,l,j,n,rootN):
              flag:=true:
            fi:
          else if k = 2 then
            if member(j,P1[l][i]) then
              P1:=solve1(P1,B,l,i,j,n,rootN):
              flag:=true:
            fi:
          else
            y:=B[i][l][1]:
            z:=B[i][l][2]:
            if member(j,P1[y][z]) then
              P1:=solve1(P1,B,y,z,j,n,rootN):
              flag:=true:
            fi:
          fi:
          fi:
        od:
      fi:
    od:
  od:
od: 
     
return[P1,flag]:    
    
end:
#######################################################
#END singleVal(P,n,rootN,B)
#######################################################


##############################################################################
#solve1(P,B,i,j,opValue,n,rootN): solve1 makes it so opValue is the solution
#  to the (i,j)th location in the possibility matrix
#  A dummy value of -1 is given to P[i][n+1][opValue] (similarly for column 
#  and block) in this case
#
#Input:  P, the possibility matrix
#        B, the configuation of the n, rootN-by-rootN blocks
#        i, the row in which we are solving
#        j, the column in which we are solving
#        opValue, the value for which we are solving
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#
#Output:  P1, the updated possibility matrix
#         or ERROR, if opValue does not appear in the (i,j)th location of P
##############################################################################
solve1:=proc(P,B,i,j,opValue,n,rootN) local P1,k,value,wB,wB2,temp,y,z:
option remember:

P1:=P:

value:={opValue}:

wB:=whichBlock(i,j,rootN):


temp:=P1[i][j]:
if not member(opValue,temp) then
  error(`Not a valid Sudoku problem`)
fi:

if nops(temp)<>1 then
  for k from 1 to nops(temp) do
    P1[i][j] := temp minus {temp[k]}:
    P1[i][n+1][temp[k]] := P1[i][n+1][temp[k]] - 1:
    P1[j][n+2][temp[k]] := P1[j][n+2][temp[k]] - 1:
    P1[wB][n+3][temp[k]] := P1[wB][n+3][temp[k]] - 1:
  od:
fi:

P1[i][j] := value:
P1[i][n+1][opValue] := -1:
P1[j][n+2][opValue] := -1:
P1[wB][n+3][opValue] := -1:


for k from 1 to n do
  temp:=P1[i][k]:
  if member(opValue,temp) then
    P1[i][k] := temp minus value:
    P1[i][n+1][opValue] := P1[i][n+1][opValue] - 1:
    P1[k][n+2][opValue] := P1[k][n+2][opValue] - 1:
    wB2:=whichBlock(i,k,rootN):
    P1[wB2][n+3][opValue] := P1[wB2][n+3][opValue] - 1:  
  fi:
  
  temp:=P1[k][j]:
  if member(opValue,temp) then
    P1[k][j] := temp minus value:
    P1[k][n+1][opValue] := P1[k][n+1][opValue] - 1:
    P1[j][n+2][opValue] := P1[j][n+2][opValue] - 1:
    wB2:=whichBlock(k,j,rootN):
    P1[wB2][n+3][opValue] := P1[wB2][n+3][opValue] - 1:
  fi:

  y:=B[wB][k][1]:
  z:=B[wB][k][2]:
  temp:=P1[y][z]:
  if member(opValue,temp) then
    P1[y][z] := temp minus value:
    P1[y][n+1][opValue] := P1[y][n+1][opValue] - 1:
    P1[z][n+2][opValue] := P1[z][n+2][opValue] - 1:
    wB2:=whichBlock(y,z,rootN):
    P1[wB2][n+3][opValue] := P1[wB2][n+3][opValue] - 1:
  fi:
od:

P1[i][j]:=value:
P1[i][n+1][opValue]:=-1:
P1[j][n+2][opValue] := -1:
P1[wB][n+3][opValue] := -1:

return P1:

end:
######################################################################
#END solve1(P,B,i,j,opValue,n,rootN)
######################################################################


###################################################################
#whichBlock(i,j,rootN): returns the ordinal of the block where the 
#  (i,j)th location falls
#
#Input:  i, the row
#        j, the column
#        rootN, the square root of the size of the Sudoku puzzle
###################################################################
whichBlock:=proc(i,j,rootN) local b,ii,jj:
option remember:

ii:=ceil(i/rootN):
jj:=ceil(j/rootN):

RETURN((ii-1)*rootN+jj):

end:
#############################################
#END whichBlock(i,j,rootN)
#############################################


##################################################################################
#champion(P,n):  finds the location in the possibility matrix with the 
#  smallest number of possibilities (>1)
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#
#Output:  a list containing:
#           minI, the row of the least possibilities
#           minJ, the column of the least possibilities
#           min, the number of possibilities in this location
#         or FAIL, if a location in the possibility matrix has ZERO possibilities
##################################################################################
champion:=proc(P,n) local i,j,min,minI,minJ,nopsPij:
option remember:

minI := 0:
minJ := 0:
min := n+1:

for j from 1 to n do
  for i from 1 to n do
    nopsPij := nops(P[i][j]):

    if nopsPij = 0 then
      return FAIL:
    fi:

    if nopsPij < min and nopsPij > 1 then
      minI := i:
      minJ := j:
      min := nopsPij:
    fi:
  od:
od:

return [minI,minJ,min]:

end:
########################################################
#END champion(P,n)
########################################################


############################################################################
#children(P,n,rootN,B,y,z): returns the children of the (y,z)th location in 
#  the possibility matrix
#
#Input:  P, the possibility matrix
#        n, the size of the Sudoku puzzle
#        rootN, the square root of n
#        B, the configuation of the n, rootN-by-rootN blocks
#        y, the row of the location to find the children
#        z, the column of the location to find the children
#
#Output:  P1, a set containing the possibility matrices of all the children
############################################################################
children:=proc(P,n,rootN,B,y,z) local i, j, P1:
option remember:

P1:={}:
for i from 1 to nops(P[y][z]) do
  P1:= P1 union {solve1(P,B,y,z,P[y][z][i],n,rootN)}:
od:

return P1:

end:
########################################################
#END(children(P,n,rootN,B,y,z)
########################################################