I have finished my cleanup of the dancing links based solver for Sudoku. I don't have time to compare performance with the other programs that have been posted recently, or to do more profiling of my code. For those who will say "It is ugly, imperative, and ugly!" please remember this is a conversion of Knuth's c-code, which depended on five non-trivial goto jumps because he did not have tail recursion. And the whole point of the algorithm are the imperative unlink & relink routines acting on a sparse binary matrix. This does not use any clever logic, but it does pick the tightest constraint at every step. This means if there is only one obvious possibility for a position/row/column/block then it will immediately act on it. The literate source file is attached. [ My clever logic solver may eventually be cleaned up as well. ] -- Chris A Sukodku solver by Chris Kuklewicz (haskell (at) list (dot) mightyreason (dot) com) I compile on a powerbook G4 (Mac OS X, ghc 6.4.2) using ghc -optc-O3 -funbox-strict-fields -O2 --make -fglasgow-exts This is a translation of Knuth's GDANCE from dance.w / dance.c http://www-cs-faculty.stanford.edu/~uno/preprints.html http://www-cs-faculty.stanford.edu/~uno/programs.html http://en.wikipedia.org/wiki/Dancing_Links I have an older verison that uses lazy ST to return the solutions on demand, which was more useful when trying to generate new puzzles to solve.
module Main where
import Prelude hiding (read) import Control.Monad import Control.Monad.Fix import Data.Array.IArray import Control.Monad.ST.Strict import Data.STRef.Strict import Data.Char(intToDigit,digitToInt) import Data.List(unfoldr,intersperse)
new = newSTRef {-# INLINE new #-} read = readSTRef {-# INLINE read #-} write = writeSTRef {-# INLINE write #-} modify = modifySTRef {-# INLINE modify #-}
Data types to prevent mixing different index and value types
type A = Int newtype R = R A deriving (Show,Read,Eq,Ord,Ix,Enum) newtype C = C A deriving (Show,Read,Eq,Ord,Ix,Enum) newtype V = V A deriving (Show,Read,Eq,Ord,Ix,Enum) newtype B = B A deriving (Show,Read,Eq,Ord,Ix,Enum)
Sudoku also has block constraints, so we want to look up a block index in an array:
lookupBlock :: Array (R,C) B lookupBlock = listArray bb [ toBlock ij | ij <- range bb ] where ra :: Array Int B ra = listArray (0,pred (rangeSize b)) [B (fst b) .. B (snd b)] toBlock (R i,C j) = ra ! ( (div (index b j) 3)+3*(div (index b i) 3) )
The values for an unknown location is 'u'. The bound and range are given by b and rng. And bb is a 2D bound.
u = V 0 -- unknown value b :: (Int,Int) b = (1,9) -- min and max bounds rng = enumFromTo (fst b) (snd b) -- list from '1' to '9' bb = ((R (fst b),C (fst b)),(R (snd b),C (snd b)))
A Spec can be turned into a parsed array with ease:
type Hint = ((R,C),V) newtype Spec = Spec [Hint] deriving (Eq,Show)
type PA = Array (R,C) V
parse :: Spec -> PA parse (Spec parsed) = let acc old new = new in accumArray acc u bb parsed
The dancing links algorithm depends on a sparse 2D node structure. Each column represents a constraint. Each row represents a Hint. The number of possible hints is 9x9x9 = 271
type (MutInt st) = (STRef st) Int
The pointer types:
type (NodePtr st) = (STRef st) (Node st) type (HeadPtr st) = (STRef st) (Head st)
The structures is a 2D grid of nodes, with Col's on the top of columns and a sparse collection of nodes. Note that topNode of Head is not a strict field. This is because the topNode needs to refer to the Head, and they are both created monadically.
type HeadName = (Int,Int,Int) -- see below for meaning
data Head st = Head {headName:: !HeadName ,topNode:: (Node st) -- header node for this column ,len:: !(MutInt st) -- number of nodes below this head ,next,prev:: !(HeadPtr st) -- doubly-linked list }
data Node st = Node {getHint:: !Hint ,getHead:: !(Head st) -- head for the column this node is in ,up,down,left,right :: !(NodePtr st) -- two doubly-linked lists }
instance Eq (Head st) where a == b = headName a == headName b
instance Eq (Node st) where a == b = up a == up b
To initialize the structures is a bit tedious. Knuth's code reads in the problem description from a data file and builds the structure based on that. Rather than short strings, I will use HeadName as the identifier. The columns are (0,4,5) for nodes that put some value in Row 4 Col 5 (1,2,3) for nodes that put Val 3 in Row 2 and some column (2,7,4) for nodes that put Val 4 in Col 7 and some row (3,1,8) for nodes that put Val 8 in some (row,column) in Block 1 The first head is (0,0,0) which is the root. The non-root head data will be put in an array with the HeadName as an index.
headNames :: [HeadName] headNames = let names = [0,1,2,3] in (0,0,0):[ (l,i,j) | l<-names,i<-rng,j<-rng]
A "row" of left-right linked nodes is a move. It is defined by a list of head names.
type Move = [(Hint,HeadName)]
Initial hints are enforced by making them the only legal move for that location. Blank entries with value 'u = V 0' have a move for all possible values [V 1..V 9].
parseSpec :: Spec -> [Move] parseSpec spec = let rowsFrom :: Hint -> [Move] rowsFrom (rc@(R r,C c),mv@(V v')) = if mv == u then [ rsyms v | v <- rng ] else [ rsyms v' ] where (B b) = lookupBlock ! rc rsyms :: A -> Move rsyms v = map ( (,) (rc,V v) ) [(0,r,c),(1,r,v),(2,c,v),(3,b,v)] in concatMap rowsFrom (assocs (parse spec))
mkDList creates doubly linked lists using a monadic smart constructor and the recursive "mdo" notation as documented at http://www.haskell.org/ghc/docs/latest/html/users_guide/syntax-extns.html#md... http://www.cse.ogi.edu/PacSoft/projects/rmb/ For more fun with this, see the wiki page at http://haskell.org/hawiki/TyingTheKnot
mkDList :: (MonadFix m) => (b -> a -> b -> m b) -> [a] -> m b mkDList _ [] = error "must have at least one element" mkDList mkNode xs = mdo (first,last) <- go last xs first return first where go prev [] next = return (next,prev) go prev (x:xs) next = mdo this <- mkNode prev x rest (rest,last) <- go this xs next return (this,last)
toSimple takes a function and a header node and iterates (read . function) until the header is reached again, but does not return the header itself.
toSingle step header = loop =<< (read . step) header where loop y = if header/=y then liftM (y:) (read (step y) >>= loop) else return []
forEach is an optimization of (toSimple step header >>= mapM_ act)
forEach step header act = loop =<< (read . step) header where loop y = if header/=y then (act y >> (read (step y)) >>= loop) else return ()
Now make the root node and all the head nodes. This also exploits mdo:
makeHeads :: [HeadName] -> (ST st) (Head st) makeHeads names = mkDList makeHead names where makeHead before name after = mdo ~newTopNode <- liftM4 (Node ((R 0,C 0),V 0) newHead) (new newTopNode) (new newTopNode) (new newTopNode) (new newTopNode) newHead <- liftM3 (Head name newTopNode) (new 0) (new after) (new before) return newHead
The Head nodes will be places in an array for easy lookup while building moves:
type HArray st = Array HeadName (Head st) hBounds = ((0,1,1),(3,9,9)) type Root st = (Head st,HArray st)
The addMove function creates the (four) nodes that represent a move and adds them to the data structure. The HArray in Root makes for a fast lookup of the Head data.
addMove :: forall st. (Root st) -> Move -> (ST st) (Node st) addMove (_,ha) move = mkDList addNode move where addNode :: (Node st) -> (Hint,HeadName) -> (Node st) -> (ST st) (Node st) addNode before (hint,name) after = do let head = ha ! name let below = topNode head above <- read (up below) newNode <- liftM4 (Node hint head) (new above) (new below) (new before) (new after) write (down above) newNode write (up below) newNode modify (len head) succ l <- read (len head) seq l (return newNode)
Create the column headers, including the fast lookup array. These will be resused between puzzles.
initHA :: (ST st) (Root st) initHA = do root <- makeHeads headNames heads <- toSingle next root let ha = array hBounds (zip (map headName heads) heads) return (root,ha)
Take the Root from initHA and a puzzle Spec and fill in all the Nodes.
initRoot :: (Root st) -> Spec -> (ST st) () initRoot root spec = do let moves = parseSpec spec mapM_ (addMove root) moves
Return the column headers to their condition after initHA
resetRoot :: (Root st) -> (ST st) () resetRoot (root,ha) = do let heads@(first:_) = elems ha let resetHead head = do write (len head) 0 let node = topNode head write (down node) node write (up node) node reset (last:[]) = do write (prev root) last write (next root) first reset (before:xs@(head:[])) = do resetHead head write (prev head) before write (next head) root reset xs reset (before:xs@(head:after:_)) = do resetHead head write (prev head) before write (next head) after reset xs reset (root:heads)
getBest iterates over the unmet constraints (i.e. the Head that are reachable from root). It locates the one with the lowest number of possible moves that will solve it, aborting early if it finds 0 or 1 moves.
getBest :: (Head st) -> (ST st) (Maybe (Head st)) getBest root = do first <- read (next root) if first == root then return Nothing else do let findMin m best head | head == root = return (Just best) | otherwise = do l <- read (len head) if l <= 1 then return (Just head) else if l < m then findMin l head =<< read (next head) else findMin l best =<< read (next head) findMin 10 first first
The unlink and relink operations are from where Knuth got the name "dancing links". So long as "a" does not change in between, the relink call will undo the unlink call. Similarly, the unconver will undo the changes of cover and unconverOthers will undo coverOthers.
unlink :: (a->STRef st a) -> (a->STRef st a) -> a -> (ST st) () unlink prev next a = do before <- read (prev a) after <- read (next a) write (next before) after write (prev after) before
relink :: (a->STRef st a) -> (a->STRef st a) -> a -> (ST st) () relink prev next a = do before <- read (prev a) after <- read (next a) write (next before) a write (prev after) a
cover :: (Head st) -> (ST st) () cover head = do unlink prev next head let eachDown rr = forEach right rr eachRight eachRight nn = do unlink up down nn modify (len $ getHead nn) pred forEach down (topNode head) eachDown
uncover :: (Head st) -> (ST st) () uncover head = do let eachUp rr = forEach left rr eachLeft eachLeft nn = do modify (len $ getHead nn) succ relink up down nn forEach up (topNode head) eachUp relink prev next head
coverOthers :: (Node st) -> (ST st) () coverOthers node = forEach right node (cover . getHead)
uncoverOthers :: (Node st) -> (ST st) () uncoverOthers node = forEach left node (uncover . getHead)
A helper function for gdance:
choicesToSpec :: [(Node st)] -> Spec choicesToSpec = Spec . (map getHint)
This is the heart of the algorithm. I have altered it to return only the first solution, or produce an error if none is found. Knuth used several goto links to do what is done below with tail recursion.
gdance :: (Head st) -> (ST st) Spec -- [Spec] gdance root = let forward choices = do maybeHead <- getBest root case maybeHead of Nothing -> if null choices then error "No choices in forward" -- return [] -- for [Spec] else do -- nextSols <- recover choices -- for [Spec] return $ (choicesToSpec choices) -- :nextSols -- for [Spec] Just head -> do cover head startRow <- readSTRef (down (topNode head)) advance (startRow:choices)
advance choices@(newRow:oldChoices) = do let endOfRows = topNode (getHead newRow) if (newRow == endOfRows) then do uncover (getHead newRow) if (null oldChoices) then error "No choices in advace" -- return [] -- for [Spec] else recover oldChoices else do coverOthers newRow forward choices
recover (oldRow:oldChoices) = do uncoverOthers oldRow newRow <- readSTRef (down oldRow) advance (newRow:oldChoices)
in forward []
Convert a text board into a Spec
parseBoard :: String -> Spec parseBoard s = Spec (zip rcs vs'check) where rcs :: [(R,C)] rcs = [ (R r,C c) | r <- rng, c <- rng ] isUnset c = (c=='.') || (c==' ') || (c=='0') isHint c = ('1'<=c) && (c<='9') cs = take 81 $ filter (\c -> isUnset c || isHint c) s vs :: [V] vs = map (\c -> if isUnset c then u else (V $ digitToInt c)) cs vs'check = if 81==length vs then vs else error ("parse of board failed\n"++s)
This is quite useful as a utility function which partitions the list into groups of n elements. Used by showSpec.
groupTake :: Int->[a]->[[a]] groupTake n b = unfoldr foo b where foo [] = Nothing foo b = Just (splitAt n b)
Make a nice 2D ascii board from the Spec (not used at the moment)
showSpec :: Spec -> String showSpec spec = let pa = parse spec g = groupTake 9 (map (\(V v) -> if v == 0 then '.' else intToDigit v) $ elems pa) addV line = concat $ intersperse "|" (groupTake 3 line) addH list = concat $ intersperse ["---+---+---"] (groupTake 3 list) in unlines $ addH (map addV g)
One line display
showCompact spec = map (\(V v) -> intToDigit v) (elems (parse spec))
The main routine is designed to handle the input from http://www.csse.uwa.edu.au/~gordon/sudoku17
main = do all <- getContents let puzzles = zip [1..] (map parseBoard (lines all)) root <- stToIO initHA let act :: (Int,Spec) -> IO () act (i,spec) = do answer <- stToIO (do initRoot root spec answer <- gdance (fst root) resetRoot root return answer) print (i,showCompact answer) mapM_ act puzzles