Monday, October 6, 2008

Let's dance

There is this sudoku craze going around the world. I got to admit, for a while I were writing numbers in squares as well. But I (quickly, I like to think.) realized that it was basically just different ways of counting which numbers are already used and filling in the rest, and repeat the process over and over and over again.

Computers are good at counting, I thought, so I wrote a program to do the work for me. Using the same logic I did myself to eliminate choices and filling in squares when only one choice were left for a square. It worked well for a while, but then I found a puzzle it didn't solve, and I had to add some more rules. And it worked for a while, but then I found another puzzle. And so on.

I'm not going to show you that program. I finally gave up on the logic approach, you never know if there is some pattern you have missed until you stumble upon it. Instead I going to show the brute force approach. It's not obvious how well that would work, sudoku puzzles with as few as 17 starting hints are known. Without doing any kind of narrowing of the search, that means there are (expt 9 (- 81 17)) combinations to search through. That is a big number.

Luckily, smarter people than me have played around with sudoku as well. And one of them, Donald E. Knuth, has come up with an algorithm that's pretty good. He calls it "Dancing links". Actually, he didn't really make it to solve sudokus, Knuth wanted to pack shapes into a container.

A small example. Suppose you have the following pieces and want to pack them onto a 3x3 square, without rotating any of them. To solve this problem you could set up a matrix so that each row represents a placement of a piece. The first 9 columns represent the squares in the 3x3 square. The last three columns represent the piece used for that row (marked L, B and S for L-shape, Big square and Small square.) We get something like this (I'm using _ to make it easier to read):

                           L B S
                      
1)  1 1 1  1 _ _  _ _ _    1 _ _    ; L in top/middle row.
2)  _ _ _  1 1 1  1 _ _    1 _ _    ; L in middle/bottom row

3)  1 1 _  1 1 _  _ _ _    _ 1 _    ; Big square in top left corner
4)  _ 1 1  _ 1 1  _ _ _    _ 1 _    ; Big square in top right corner
5)  _ _ _  1 1 _  1 1 _    _ 1 _    ; Big square in bottom left corner
6)  _ _ _  _ 1 1  _ 1 1    _ 1 _    ; Big square in bottom right corner

7)  1 _ _  _ _ _  _ _ _    _ _ 1    ; Small square top left
8)  _ 1 _  _ _ _  _ _ _    _ _ 1    ; Small square top middle 
9)  _ _ 1  _ _ _  _ _ _    _ _ 1    ; etc.
10) _ _ _  1 _ _  _ _ _    _ _ 1
11) _ _ _  _ 1 _  _ _ _    _ _ 1
12) _ _ _  _ _ 1  _ _ _    _ _ 1
13) _ _ _  _ _ _  1 _ _    _ _ 1
14) _ _ _  _ _ _  _ 1 _    _ _ 1
15) _ _ _  _ _ _  _ _ 1    _ _ 1

Now we can restate the problem. We will select rows in such a way that if we add up the selected rows we get one and only one one (1) in each column. In this example row 1, 6 and 13 solves the problem. Note also that the last three columns while not representing a square that needs to be filled on the board, is treated no differently.

We can look at each column as a constraint. The last three columns forces us to use each piece only once. And the first nine columns forces us to use each square of the board only once.

Now if you're anything like me you would solve this by first picking the first row, then mentally remove all rows that have a 1 in any of the columns where the first row has a 1. Then mentally remove the columns where the row has a 1 and then looking for the first row that has a 1 in the remaining columns, then repeat that process until there are no more columns to fill. If you had started out with the third row instead of the first row, you would at one point have realized that you couldn't solve the puzzle, and you would put the columns rows back in and chosen another row. You could just start over, but that's inefficient for big puzzles.

Let's make another small puzzle, fill a 3x3 square, and use the numbers 1, 2 and 3 in it so that each row and column have all three numbers in them, but only once in each row (Obviously). I'll make a matrix based on the row, column and value of the number in the square like this:


r,c=v         r,c                r,v                c,v
             
1,1=1   1 _ _ _ _ _ _ _ _  1 _ _ _ _ _ _ _ _  1 _ _ _ _ _ _ _ _
1,1=2   1 _ _ _ _ _ _ _ _  _ 1 _ _ _ _ _ _ _  _ 1 _ _ _ _ _ _ _
1,1=3   1 _ _ _ _ _ _ _ _  _ _ 1 _ _ _ _ _ _  _ _ 1 _ _ _ _ _ _
1,2=1   _ 1 _ _ _ _ _ _ _  1 _ _ _ _ _ _ _ _  _ _ _ 1 _ _ _ _ _
1,2=2   _ 1 _ _ _ _ _ _ _  _ 1 _ _ _ _ _ _ _  _ _ _ _ 1 _ _ _ _
1,2=3   _ 1 _ _ _ _ _ _ _  _ _ 1 _ _ _ _ _ _  _ _ _ _ _ 1 _ _ _
1,3=1   _ _ 1 _ _ _ _ _ _  1 _ _ _ _ _ _ _ _  _ _ _ _ _ _ 1 _ _
1,3=2   _ _ 1 _ _ _ _ _ _  _ 1 _ _ _ _ _ _ _  _ _ _ _ _ _ _ 1 _
1,3=3   _ _ 1 _ _ _ _ _ _  _ _ 1 _ _ _ _ _ _  _ _ _ _ _ _ _ _ 1

             
2,1=1   _ _ _ 1 _ _ _ _ _  _ _ _ 1 _ _ _ _ _  1 _ _ _ _ _ _ _ _
2,1=2   _ _ _ 1 _ _ _ _ _  _ _ _ _ 1 _ _ _ _  _ 1 _ _ _ _ _ _ _
2,1=3   _ _ _ 1 _ _ _ _ _  _ _ _ _ _ 1 _ _ _  _ _ 1 _ _ _ _ _ _
2,2=1   _ _ _ _ 1 _ _ _ _  _ _ _ 1 _ _ _ _ _  _ _ _ 1 _ _ _ _ _
2,2=2   _ _ _ _ 1 _ _ _ _  _ _ _ _ 1 _ _ _ _  _ _ _ _ 1 _ _ _ _
2,2=3   _ _ _ _ 1 _ _ _ _  _ _ _ _ _ 1 _ _ _  _ _ _ _ _ 1 _ _ _
2,3=1   _ _ _ _ _ 1 _ _ _  _ _ _ 1 _ _ _ _ _  _ _ _ _ _ _ 1 _ _
2,3=2   _ _ _ _ _ 1 _ _ _  _ _ _ _ 1 _ _ _ _  _ _ _ _ _ _ _ 1 _
2,3=3   _ _ _ _ _ 1 _ _ _  _ _ _ _ _ 1 _ _ _  _ _ _ _ _ _ _ _ 1

             
3,1=1   _ _ _ _ _ _ 1 _ _  _ _ _ _ _ _ 1 _ _  1 _ _ _ _ _ _ _ _
3,1=2   _ _ _ _ _ _ 1 _ _  _ _ _ _ _ _ _ 1 _  _ 1 _ _ _ _ _ _ _
3,1=3   _ _ _ _ _ _ 1 _ _  _ _ _ _ _ _ _ _ 1  _ _ 1 _ _ _ _ _ _
3,2=1   _ _ _ _ _ _ _ 1 _  _ _ _ _ _ _ 1 _ _  _ _ _ 1 _ _ _ _ _
3,2=2   _ _ _ _ _ _ _ 1 _  _ _ _ _ _ _ _ 1 _  _ _ _ _ 1 _ _ _ _
3,2=3   _ _ _ _ _ _ _ 1 _  _ _ _ _ _ _ _ _ 1  _ _ _ _ _ 1 _ _ _
3,3=1   _ _ _ _ _ _ _ _ 1  _ _ _ _ _ _ 1 _ _  _ _ _ _ _ _ 1 _ _
3,3=2   _ _ _ _ _ _ _ _ 1  _ _ _ _ _ _ _ 1 _  _ _ _ _ _ _ _ 1 _
3,3=3   _ _ _ _ _ _ _ _ 1  _ _ _ _ _ _ _ _ 1  _ _ _ _ _ _ _ _ 1

If we now do as above and choose rows out of this matrix such that we get one and only one one (1) for each column, you can see that the first 9 columns of the matrix forces us to choose only one value for each pair r,c. The next 9 columns forces us to use each value only once for each row (If I chose the first line, (put "1" in location "1,1", I can't choose line 4 (put "1" in "1,2")). And the last columns forces us to only use each value once for each column.

This is almost how I'm going to represent the sudoku puzzle, soduku have one more constraint, the 9x9 square is divided into 9 3x3 blocs, and each block also have to contain the numbers 1-9. This is easy enough to accomplish. We just add another block to our matrix, so in addtion to (r,c) (r,v) (c,v) we have (b,v). This means our matrix is big enough that I'm not going to bother writing it out. Each row has 81 columns for the (r,c) constraint, 81 columns for (r,v), 81 columns for (c,v) and 81 columns for (b,v) => 364 columns. And the matrix has one row for each combination of r,c,v => 9*9*9 => 729 rows.

The point I wanted to illustrate with this detour before I start describing the algorithm is that you don't have to come up with a completely new algorithm for a problem if you know a good algorithm that can be used if you look at the problem in a different way.

When representing a matrix, it could be tempting to represent it as an array. It is probaly the most common way of representing matrixes. However, Knuth used double-linked circular lists for this algorithm, and for a good reason. When you remove an element from a double-linked list, the element remembers it's position in the list if you don't change the elements next/previous pointers. This makes it very easy to reinsert the element again. So when we guess wrong in our algorithm it is easy to backtrack a few steps.

The only thing to watch out for is to reinsert the elements in the opposite order they were taken out of the list, which is what you would do in a depth first search anyway.

So what is the data structure going to look like?

  • Each "1" will be represented by a matrix element, we will not represent the gaps in the matrix in any way. They are not interesting.
  • Each element get a left and right pointer which points to the previous/next element in the row. It's a circular list, so the first and the last element in each row points to each other.
  • Each element also get a up/down pointer with points to the previous/next element in the column. The top and bottom element points to a special element that is the header for that column.
  • Header elements have right/left pointers to each other in a doubly linked list as well. The first and the last header element points to a node that is our entry point to the whole matrix.
  • The header elements also has a name slot, which is not really neccessary but useful during debugging, and a size slot, which is the number of elements left in that column. We'll choose columns to search based on how many elements are left in the columns. If we search the column with the smallest number of elements first, we will reduce the search tree and the whole thing will go a lot faster.
  • Each node in the matrix also has extra slots, one triplet representing (row column value). This is used to print out the solution. And they also have a pointer to the header element of the column they are in.
  • Finally, I'm going to make a array with one element for each row in the matrix, this will point to the first element on each row, and it will be used when setting up the starting position.
The whole thing should look someting like, oh, buggr'it, go read Knuths paper. It's got pretty pictures and stuff. Maybe I'll be bored enough one day to add a diagram here.

With that out of the way, on to the implementation:


(defpackage :sudoku
    (:use :common-lisp))
  
(in-package :sudoku)

(defclass matrix-element ()
  ((up :accessor up :initarg :up)
   (down :accessor down :initarg :down)
   (left :accessor left :initarg :left)
   (right :accessor right :initarg :right)))

Now, if I were very concerned about performance I might use something else than classes for this, but Knuth says premature optimization is evil, so I wont. Besides, the best optimization is to use a good algorithm, which I think I'm doing. It's not like I made it myself.

(defmethod initialize-instance :after ((me matrix-element) &rest initargs &key &allow-other-keys)
  (setf (up me) me)
  (setf (down me) me)
  (setf (left me) me)
  (setf (right me) me))

Once nice thing about classes is that you can make them do stuff when you initialize them by adding an :after method on the initialize-instance method. What I do above is to initialize all the pointers of the matrix element to point back at the newly created element instance. This means that a newly created instance is a tiny circular list with itself as the only object and I don't have to worry about wheter a single element is a list or not, I just treat everything as lists.

(defclass column-header (matrix-element)
  ((size :accessor size :initform 0)
   (name :accessor name :initarg :name)))

(defclass column-element (matrix-element)
  ((column-header :accessor column-header :initarg :column-header)
   (triplet :accessor triplet :initarg :triplet)))

I'm creating two subclasses, to add slots for the column headers and the matrix elements as I mentioned above.

(defun insert-above (element new)
  "Insert a new element above the existing element. And increase size of column by 1"
  (setf (up new) (up element)
        (down new) element
        (down (up element)) new
        (up element) new)
  (incf (size (column-header new)))
  element)

(defun insert-rightof (element new)
  (setf (right new) (right element)
        (left new) element
        (left (right element)) new
        (right element) new)
  element)

Our insertion functions. I'm building up the matrix right to left and top to bottom. Since the columns are circular lists the easiest way to insert an element at the bottom of the list is to insert it above the column header.

(defmethod initialize-instance :after ((me column-element) &rest initargs &key &allow-other-keys)
  (if (column-header me)
    (insert-above (column-header me) me)))

When I initialize the matrix elements I'm going to pass the column header as an argument to the make-instance call. This method puts the element into the circular list for that column.

(defun remove-horizontal (element)
  (setf (left (right element)) (left element)
        (right (left element)) (right element))
  element)

(defun replace-horizontal (element)
  (setf (left (right element)) element
        (right (left element)) element)
  element)

The core of the algorithm, removing and replacing an element from a list. These elements assumes that when you reinsert the element, the list is in the same state as it was when the element was removed, in other words, that (left element) and (right element) points to each other.
(defun remove-vertical (element)
  (setf (up (down element)) (up element)
        (down (up element)) (down element))
  (decf (size (column-header element)))
  element)

(defun replace-vertical (element)
  (setf (up (down element)) element
        (down (up element)) element)
  (incf (size (column-header element)))  
  element)

Same thing vertically, this time I update the size of the column as well.

(defmacro rcv-loop (&rest body)
  `(loop for r from 0 to 8
         do (loop for c from 0 to 8
                  do (loop for v from 0 to 8
                           do ,@body))))


A macro that I'm going to use to build up the matrix. This is probably bad style since I rely on body using the variables r c and v in the body.

(defun add-constraint-headers (matrix elem1 elem2)
  (loop for i from 0 to 8
        do (loop for j from 0 to 8
            do (insert-rightof (left matrix) 
                               (make-instance 'column-header
                                              :name (format nil "~a: ~a/~a: ~a" elem1 (+ 1 i) elem2 (+ 1 j)))))))

(defun sudoinit ()
  "returns an initial matrix representing an empty board."
  (let ((matrix (make-instance 'matrix-element))
        (col-array (make-array (* 4 9 9)))
        (row-array (make-array (* 9 9 9))))
    (add-constraint-headers matrix "Row" "Col")
    (add-constraint-headers matrix "Row" "Number")      
    (add-constraint-headers matrix "Col" "Number")      
    (add-constraint-headers matrix "Block" "Number")
    (loop with foo = (right matrix)
          for i from 0
          until (eql foo matrix)
          do (setf (aref col-array i) foo
                   foo (right foo)))
    (rcv-loop
     (let* ((triplet (list r c v))
            (r-c (make-instance 'column-element
                                :column-header (aref col-array (+ (* 9 r)
                                                                  c))
                                :triplet triplet))
            (r-v (make-instance 'column-element
                                :column-header (aref col-array (+ (* 9 9)
                                                                  (* 9 r)
                                                                  v))
                                :triplet triplet))
            (c-v (make-instance 'column-element
                                :column-header (aref col-array (+ (* 2 9 9)
                                                                  (* 9 c)
                                                                  v))
                                :triplet triplet))
            (b-v (make-instance 'column-element
                                :column-header (aref col-array (+ (* 3 9 9)
                                                                  (* 9 (+ (* 3 (floor r 3))
                                                                          (floor c 3)))
                                                                  v))
                                :triplet triplet)))
       (insert-rightof r-c r-v)
       (insert-rightof r-v c-v)
       (insert-rightof c-v b-v)
       (setf (aref row-array (+ (* 9 (+ (* 9 r) c)) v)) r-c)))
    (values matrix row-array)))

Ok, this is the big one. I first create the row of header elements, adding a short descriptive name in each element. I collect them in a array as well, so it is easy to reference the headers when adding the elements later.

Each row in the matrix has 4 elements. One for the row/col constraint, one for row/value, one for col/value and one for block/value. Remember that the make-instance method for column-elements automatically inserts the element in the right column if I pass it a :column-header argument.

Finally I arrange the elements in a row and put a pointer to the row in my row-array.


(defun cover-col (col-header)
  (do ((elem (down col-header) (down elem)))
      ((eql elem col-header) col-header)
    (do ((row-elem (right elem) (right row-elem)))
        ((eql row-elem elem) row-elem)
      (remove-vertical row-elem)))
  (remove-horizontal col-header))

Knuth calls the hiding of columns/rows for covering, so I'm going to use the same name. This function takes a column header, and for each row that thas an element in this column, it will hide the row vertically so that this row can not be chosen later. So if col-header is the header of the first column, this will remove all rows with a element in the first column.

(defun uncover-col (col-header)
  (do ((elem (up col-header) (up elem)))
      ((eql elem col-header) col-header)
    (do ((row-elem (right elem) (right row-elem)))
        ((eql row-elem elem) row-elem)
      (replace-vertical row-elem)))
  (replace-horizontal col-header))

Uncovering the column again. Both this and the previous function test for the equality of the current and initial element to detect that the circular list has been traversed.


(defun cover-row (elem)
  (do ((column-iter (right elem) (right column-iter)))
      ((eql elem column-iter))
    (cover-col (column-header column-iter))))

(defun uncover-row (elem)
  (do ((column-iter (left elem) (left column-iter)))
      ((eql elem column-iter))
    (uncover-col (column-header column-iter))))

When we choose a move we need to remove the colliding moves for all the columns where the row have elements. This functions takes care of this by iterating over the row, and calling cover-col for each of the columns.


(defun try-elem (elem array solution)
  (push elem solution)
  (cover-row elem)
  (let ((foo (find-solution array solution)))
    (uncover-row elem)
    (pop solution)
    foo))

Why the let? I want to return the value returned by find-solution to the caller because the return value signals whether a solution has been found.


(defun find-smallest-col (array)
  (do ((head-iter (right array) (right head-iter))
       (cur-min (right array) (if (< (size head-iter) (size cur-min))
                                head-iter
                                cur-min)))
      ((or (= 1 (size cur-min))
           (eql head-iter array)) cur-min)))

Loop through the headers and find the column with the smallest number of elements. If we find a column with only one element, we can stop searching.


(defun find-solution (array solution)
  (if (eql (left array) array)
    (progn
      (print-solution solution)
      t)
    (let ((col-header (find-smallest-col array)))
      (prog2
          (cover-col col-header)
          (do ((row-elem (down col-header) (down row-elem)))      
              ((or (eql row-elem col-header)
                   (try-elem row-elem array solution))
               (not (eql row-elem col-header))))
        (uncover-col col-header)))))

Try to find a solution. If all the headers are hidden, we have found a solution, and can print it out. Return t to signal that the solution is found. If a solution is not found yet, find the column with the smallest number of elements and try out each element to see if we find a solution. If a solution is found try-elem will return true, and we can stop searching. Make sure to pass on the true value to the calling function.


(let ((solc 0))
  (defun print-solution (solution)
    (format t "  ~a  " (incf solc))
    (loop with arr = (make-array '(9 9))
          for s in solution
          do (destructuring-bind (i j d) (triplet s)
               (setf (aref arr i j) (+ 1 d)))
          finally (loop for i from 0 to 8
                        do (format t "~{~A~}"
                                   (loop for j from 0 to 8
                                         collect (aref arr i j))))))
  
  (defun reset-count ()
    (setf solc 0)))

Very simple way of printing out solutions. I just build up the resulting matrix from the triplets in the solution list and then print out the resulting matrix on a single line. I also maintain a counter that I can reset counting how many puzzles has been solved.


(defun read-starting-position (pos-string)
  (with-input-from-string (s pos-string)
    (remove-if 'null
               (apply 'append
                      (loop for r from 0 to 8
                            collect (loop for c from 0 to 8
                                          collect (let ((v (- (char-code (read-char s))
                                                              (char-code #\0))))
                                                    (when (> v 0)
                                                      (list r c (1- v))))))))))

Our starting positions is just a string of numbers where 0 means an open square. This function collects a list of triplets (r c v) for all non-null positions. Note the complete lack of any validation of input. I have to subtract 1 from the value because I use the value as a index into the matrix and indexes start at 0.


(defun solve-file (filename)
  (reset-count)
  (multiple-value-bind (matrix rowarr) (sudoinit)
    (with-open-file (s filename)
      (do ((string (read-line s nil 'eof) (read-line s nil 'eof)))
          ((eql string 'eof) nil)
 (format t "~%~a" string)
        (solve string matrix rowarr)))))

(defun solve (pos-string &optional matrix rowarr)
  (unless (and matrix rowarr)
    (multiple-value-bind (foo bar) (sudoinit)
      (setf matrix foo)
      (setf rowarr bar)))
  (let ((solution ()))
    (dolist (triplet (read-starting-position pos-string))
      (destructuring-bind (i j d) triplet
        (push (aref rowarr (+ (* 9 (+ (* 9 i) j)) d)) solution)
        (cover-col (column-header (aref rowarr (+ (* 9 (+ (* 9 i) j)) d))))
        (cover-row (aref rowarr (+ (* 9 (+ (* 9 i) j)) d)))))
    (find-solution matrix solution)
    (do ((elem (pop solution) (pop solution)))
        ((null elem) matrix)
      (uncover-row elem)
      (uncover-col (column-header elem)))))

Finally, getting to run the whole thing. solve-file just reads in a file a line at a time and passes it on to the solve funtion for processing. I create the matrix and the row-array in solve-file as well, otherwise I would have to recreate the matrix for every puzzle.

The solve function is quite simple, read in the starting position. apply the moves in the starting position as if it is moves we have searched so far. Pass everything on to find-solution, and clean up the matrix for reuse after find-solution returns.

So how good is it? Remember a long time ago, I said something about the number of possibilities being a very large number? On my machine, which by no means is a top of the line machine, this program manages to solve about 200 puzzles/second. And that is without doing much in the way of optimizing.

2 comments:

David said...

Just to give you a point of comparison, I implemented a brute force method (using recursion with backtracking) and my C program running on an Intel dual-core 2.2GHz (E4500) computer solves a 9x9 Sudoku puzzle in about 30ms. This was far easier than researching algorithms to do the job. I like your investigation of the Knuth dancing links algorithm.

Using my solver, I generated Sudoku puzzles to solve which I make available at my site http://this1that1whatever.com/miscellany/sudoku/.

Unknown said...

The Most Successful Sites for Crypto, Casino & Poker - Goyang
Goyang Casino & Poker filmfileeurope.com is one of the most famous and well known crypto bsjeon.net gambling sites, founded in 2012. They are popular 출장샵 because goyangfc of their 1xbet app great