A while back there was a long thread about Sudoku solvers (some of which ended
up on http://haskell.org/haskellwiki/Sudoku ). I contributed my brute-force
dancing links solver at the time, and mentioned that I had a by-logic solver
that, while a bit slow, was as good as most of those being discussed.
At the time the code for my solver was too ugly to post. Attached is a cleaned
up version.
I have gone back and rewritten it, and come to the conclusion: There are only
two deduction algorithms: "subsets" and "blocks". These subsume the other types
of propagation and deduction. So I made this version as a "minimalist" example
instead of going for performance.
The "subsets" algorithm can be applied to all 6 permutations of row column and
value, as well as 1 special case of value and block indices.
The "blocks" algorithm can be applied 4 ways (in two flavors and to either
permutation of row/column or column/row).
There are newtypes for row, column, value, block index, and sub-block index.
The state is held in an array of type DiffArray (R,C,V) Bool
The actual computation is a series of concat/map/filter/group/sort
operations on the assocs's of the array.
The choice of which permutation is handled by leaning on the type system to
reify the type into appropriate view,shuffle, and unshuffle functions.
It should solve exactly the same number of puzzles as my older version, where "I
sent the 36628 line sudoku17 puzzle through it and it could solve 31322 of the
puzzles, leaving 5306 resistant."
--
Chris Kuklewicz
module Main (main) where
import Data.Ix(inRange,range)
import Data.Char(intToDigit,digitToInt)
import Deduce(deduce,lo,hi)
loC = intToDigit lo
hiC = intToDigit hi
unsetC = pred loC
main = do
all <- getContents
let puzzles = zip [1..] (map parseBoard (lines all))
act (i,p) = do p' <- deduce p
return (i,length p,length p')
mapM_ (\ip -> act ip >>= print) puzzles
parseBoard :: String -> [(Int,Int,Int)]
parseBoard s = map toHint justSet
where rcs = [ (r,c) | r <- range (lo,hi), c <- range (lo,hi) ]
isHint vC = inRange (loC,hiC) vC
justSet = filter (isHint . snd) (zip rcs s)
toHint ((r,c),vC) = (r,c,digitToInt vC)
{- By Chris Kuklewicz -}
module Deduce (deduce,lo,hi) where
{- The exposed function deduce takes a list of (row,column,value)
tuples that are the known parts of the solutions and returns a
(hopefully longer) list in the same format. The indices can be any
enumerated type in the range (lo,hi).
-}
import Data.Array.Diff (assocs,(!),(//),ixmap,range,inRange,accumArray,DiffArray,Ix)
import Data.List(sortBy,groupBy,transpose,(\\))
import Control.Monad(liftM,guard)
default ()
-- Typesafe values for indices
-- This machinery allows for more type safety than if R,C,V,B,D were all Int or Char
type E = Int
newtype R = R E deriving (Eq,Ord,Ix,Enum,Show) -- Row index
newtype C = C E deriving (Eq,Ord,Ix,Enum,Show) -- Column index
newtype V = V E deriving (Eq,Ord,Ix,Enum,Show) -- Value index
newtype B = B E deriving (Eq,Ord,Ix,Enum,Show) -- 3x3 Block index
newtype D = D E deriving (Eq,Ord,Ix,Enum,Show) -- Inside 3x3 Block index
lo,hi :: (Enum a) => a
lo = toEnum 1
hi = toEnum 9
fullRange :: (Enum a,Ix a) => [a]
fullRange = range (lo,hi)
rcToBD (R r) (C c) = let (rq,rm) = quotRem (r-lo) 3
(cq,cm) = quotRem (c-lo) 3
b = lo + ( 3*rq + cq )
d = lo + ( 3*rm + cm )
in (B b, D d)
bdToRC (B b) (D d) = let (bq,bm) = quotRem (b-lo) 3
(dq,dm) = quotRem (d-lo) 3
r = lo + ( 3*bq + dq )
c = lo + ( 3*bm + dm )
in (R r, C c)
-- Typeclasses and Data for "shuffle" and "unshuffle"
class (Show x, Ix x, Enum x, Ord x) => IE x
instance IE R; instance IE C; instance IE V; instance IE B; instance IE D
data Perms a b c = Perms { shuffle' :: (R,C,V) -> (a,b,c)
, unshuffle' :: (a,b,c) -> (R,C,V) }
-- Reify the types "a b c" to a value of type Perms
class (IE a, IE b, IE c) => Perm a b c where perm :: Perms a b c
instance Perm R C V where perm = Perms id id
instance Perm R V C where perm = Perms (\ (r,c,v) -> (r,v,c)) (\ (r,v,c) -> (r,c,v))
instance Perm C V R where perm = Perms (\ (r,c,v) -> (c,v,r)) (\ (c,v,r) -> (r,c,v))
instance Perm C R V where perm = Perms (\ (r,c,v) -> (c,r,v)) (\ (c,r,v) -> (r,c,v))
instance Perm V R C where perm = Perms (\ (r,c,v) -> (v,r,c)) (\ (v,r,c) -> (r,c,v))
instance Perm V C R where perm = Perms (\ (r,c,v) -> (v,c,r)) (\ (v,c,r) -> (r,c,v))
-- Special cases
instance Perm B D V where perm = Perms (\ (r,c,v) -> let (b,d) = rcToBD r c in (b,d,v))
(\ (b,d,v) -> let (r,c) = bdToRC b d in (r,c,v))
instance Perm V B D where perm = Perms (\ (r,c,v) -> let (b,d) = rcToBD r c in (v,b,d))
(\ (v,b,d) -> let (r,c) = bdToRC b d in (r,c,v))
shuffle :: (Perm a b c) => (R,C,V) -> (a,b,c)
shuffle = shuffle' perm
unshuffle :: (Perm a b c) => (a,b,c) -> (R,C,V)
unshuffle = unshuffle' perm
-- Array types, values and functions
type Index = (R,C,V)
type Cell = Bool
on,off :: Cell
on = True -- Means this might be part of the solution
off = False -- Means this cannot be part of the solution
boundsCube :: (Perm a b c) => ((a,b,c),(a,b,c))
boundsCube = ((lo,lo,lo),(hi,hi,hi))
type Cube = DiffArray Index Cell
emptyCube :: Cube
emptyCube = accumArray const on boundsCube []
type View a b c = DiffArray (a,b,c) Cell
{-# INLINE view #-}
view :: (Perm a b c) => Cube -> View a b c
view cube = ixmap boundsCube unshuffle cube
type Hints = [(Index,Cell)]
isOn :: (Perm a b c) => View a b c -> [(a,b,c)]
isOn = map fst . filter snd . assocs
-- The goal is to create functions that turn the current Cube into a
-- list of Hints. These Hints will be purely subtractive: they all
-- turn a Cell from 'on' to 'off'.
type Rule = Cube -> Hints
-- Small utility functions
fst3 (x,_,_) = x; snd3 (_,x,_) = x; thd3 (_,_,x) = x
fst4 (x,_,_,_) = x; snd4 (_,x,_,_) = x; thd4 (_,_,x,_) = x; fth4 (_,_,_,x) = x
by un bi = (\ left right -> (un left) `bi` (un right))
sortWith un = sortBy (by un compare)
groupWith un = groupBy (by un (==))
groupSort un = groupWith un . sortWith un
atLeastOne = not . null
atLeastTwo (_:_:_) = True; atLeastTwo _ = False
exactlyOne [_] = True; exactlyOne _ = False
oneOrTwo [_] = True; oneOrTwo [_,_] = True; oneOrTwo _ = False
{- ruleBlock1 : When operating on Perm V R C:
Given a value V, look along each row and see which blocks that value
may occupy. Find a row R1 for which the value is allowed in exactly
one block B1 (an no other blocks). This occupies [(B1,d)] in
B1. Eliminate V from the other locations in B1.
Given a value V, look along each row and see which blocks that value
may occupy. Find two rows [R1,R2] for which the value only is allowed
in exactly the same two blocks [B1,B2] (and no others). These occupy
[(B1,d11s)] in R1, [(B2,d12s)] in R1, [(B1,d21s)] in R2, and
[(B2,d22s)] in R2. Eliminate V from the other locations in B1 and the
other locations in B2.
Works for R and C reversed, of course.
-}
{- ruleBlock2 : When operating on Perm V R C:
Given a value V, look inside each block and see which rows that value
may occupy. Find a block B1 for which only one row R1 is occupied
(and no other rows). This occupies [(R1,c11s)] in B1. Eliminate V
from all the other c's in row R1.
Given a value V, look inside each block and see which rows that value
may occupy. Find two blocks [B1,B2] for which the value only is allowed
in exactly the same two rows [R1,R2] (and no others). These occupy
[(R1,c11s)] in B1, [(R1,c12s)] in B2, [(R2,c21s)] in B1, and [(R2,c22s)] in
B2 for some C's. Eliminate V from the other rows in R1 and the
other columns in R2.
Works for R and C reversed, of course.
-}
{- ruleBlockP
There is enough similarity between ruleBlock1 and ruleBlock2 to parameterize
over expand and contract.
The V index is special because it is "orthogonal" to the 3x3 blocks.
The line
sVsA2 = map (filter atLeastTwo) $ sVsA -- drop unique ones
is used to prevent propagating solved constraints, as rule4P should
already do this in its k==1 case.
(assocs view) becomes useful hints through a chain of map, sort, group,
filter, and concat operations. The assemble function also takes
care to remove redundant hints by consulting the view.
-}
{-# INLINE ruleBlockP #-}
ruleBlockP :: forall a b c x y . (IE a, IE b, IE c, Perm V x y)
=> ( (V,x,y) -> (V,a,b,c) ) -- "expand"
-> (V -> b -> [c] -> [Index]) -- "contract"
-> View V x y -> Hints
ruleBlockP expand contract view =
let allOn :: [(V,a,b,c)]
allOn = map expand . isOn $ view
sV :: [ [(V,a,b,c)] ]
sV = groupWith fst4 $ allOn -- group by V
sVsA :: [[ [(V,a,b,c)] ]]
sVsA = map (groupSort snd4) $ sV -- group by a
sVsAsB,sVsAsB2 :: [[[ [(V,a,b,c)] ]]]
sVsAsB = map (map (groupSort thd4)) $ sVsA -- group by b
-- These filters are to remove empty and redundantly full possibilities
sVsAsB2 = filter atLeastOne . map (filter oneOrTwo) $ sVsAsB
sVsAsBgB :: [[ ([b],[[(V,a,b,c)]]) ]]
-- The filter is to remove solved parts of the puzzle (punt to ruleSubsetP)
sVsAsBgB = map (map getAllB . filter (atLeastTwo . concat)) $ sVsAsB2
where getAllB :: [[(V,a,b,c)]] -> ([b],[[(V,a,b,c)]])
getAllB vabcss = (map (thd4 . head) vabcss, vabcss) -- length (concat vabcss) >= 2
useful :: [ [([b],[[(V,a,b,c)]])] ]
useful = concatMap (filter exactlySame . groupSort fst) $ sVsAsBgB
where exactlySame :: [([b],[[(V,a,b,c)]])] -> Bool
exactlySame sas@((sbs,_):_) = length sas == length sbs
assemble :: [ ([b],[[(V,a,b,c)]]) ] -> Hints
assemble stuff = [ (rcv,off) | rcv <- ixs, (view ! shuffle rcv) /= off ]
where byVB :: [[(V,a,b,c)]]
byVB = map concat . transpose . map snd $ stuff -- Regroup by identical 'b'
-- byVB = groupSort thd4 . concat . concat . map snd $ stuff -- equivalent
act :: [(V,a,b,c)] -> [Index]
act allVB@((v,_,b,_):_) = contract v b (fullRange \\ map fth4 allVB)
ixs :: [Index]
ixs = concatMap act byVB
in concatMap assemble useful
{-# INLINE ruleBlock1 #-}
ruleBlock1 :: (Perm V x y) => View V x y -> Hints
ruleBlock1 = ruleBlockP expand contract
where expand all@(v,x,y) = let (r,c,_) = unshuffle all
(b,d) = rcToBD r c
in (v,x,b,d)
contract v b ds = map undo ds
where undo d = let (r,c) = bdToRC b d in (r,c,v)
{-# INLINE ruleBlock2 #-}
ruleBlock2 :: (Perm V x y) => View V x y -> Hints
ruleBlock2 = ruleBlockP expand contract
where expand all@(v,x,y) = let (r,c,_) = unshuffle all
(b,_) = rcToBD r c
in (v,b,x,y)
contract v x ys = map undo ys
where undo y = unshuffle (v,x,y)
ruleBlocks :: [Rule]
ruleBlocks = [ (\ known -> ruleBlock1 (view known :: View V R C) )
, (\ known -> ruleBlock1 (view known :: View V C R) )
, (\ known -> ruleBlock2 (view known :: View V R C) )
, (\ known -> ruleBlock2 (view known :: View V C R) )
]
{- Given a list of locations, such as for the 9 columns of a row,
look at the allowed values at each location. Find a subset of k
columns for which the union of their allowed values [V..] has
length k. Then eliminate [V..] from the (9-k) other columns.
This clearly finds a list of N columns each with the same N values
if such a thing exists, so it subsumes rule1P.
This is fully symmetric in R C and V and depends on the constaints in
R and V but not C. So the (View B D V) case also works.
A useful property of this rule is that once there is only one way
to place a value in a row or column or block then it will propagate
that solution to the related contraints. This is the case when
minK is 1.
-}
{-# INLINE ruleSubsetP #-}
ruleSubsetP :: forall a b c.(Perm a b c) => View a b c -> Hints
ruleSubsetP view =
let allOn :: [(a,b,c)]
allOn = isOn $ view
sAsB :: [[ [(a,b,c)] ]]
sAsB = map (groupWith snd3) . groupWith fst3 $ allOn
sAsBgC :: [ [(Int,[c],[(a,b,c)])] ]
sAsBgC = map (sortWith fst3 . map (\ abcs -> (length abcs
,map thd3 abcs
,abcs) ) ) $ sAsB
makeChains :: [(Int,[c],[(a,b,c)])] -> Int -> [(a,[b],[c])]
makeChains input k = filter check . map toChain . subsets k . upToK $ input
where upToK :: [(Int,[c],[(a,b,c)])] -> [(Int,[c],[(a,b,c)])]
upToK = takeWhile ((k>=).fst3)
toChain :: [(Int,[c],[(a,b,c)])] -> (a,[b],[c])
toChain vals = ( fst3 . head . thd3 . head $ vals -- record "a" for easy retrieval later
, map (snd3 . head . thd3) $ vals -- the bs, length [b] == length vals == k by property of subsets k
, combine . map snd3 $ vals) -- union of the cs at each b in bs
check :: (a,[b],[c]) -> Bool
check (_,_,cs) = (k == length cs) -- check that length [c] == k as well
getUseful :: [(Int,[c],[(a,b,c)])] -> [(a,[b],[c])]
getUseful [] = []
getUseful input = concatMap (makeChains input) [minK .. maxK]
-- assertion: (length input) == (length . combine . map snd3 $ input)
where minK,maxK :: Int
minK = fst3 . head $ input
maxK = pred . length $ input
useful :: [(a,[b],[c])]
useful = concatMap getUseful sAsBgC
{- The chains are (a,bs,cs) such that
cs == nub . sort $ [ z | (x,y,z) <- allOn, x==a, y `elem` bs]
and
length bs == length cs.
Thus all the (a,b in bs,c) in the final puzzle have distict c in cs.
For location (a,b' not in bs,c') cannot have c' in cs,
thus if c' is in cs then (a,b',c') should be turned off.
-}
assemble :: (a,[b],[c]) -> Hints
assemble (a,inBs,inCs) = do -- List Monad
b <- fullRange \\ inBs
c <- inCs
let abc :: (a,b,c)
abc = (a,b,c)
guard (view ! abc)
return (unshuffle abc,off)
in concatMap assemble useful
-- All subsets of length 'k', order is stable
subsets :: Int -> [a] -> [[a]]
subsets 0 _ = [[]]
subsets _ [] = []
subsets k (x:xs) = (fmap (x:) (subsets (pred k) xs)) ++ subsets k (xs)
-- Hopefully efficient merge of (list of (sorted lists)), unique values only
combine :: (Ord a) => [[a]] -> [a]
combine [] = []
combine [x] = x
combine xs = let (a,b) = split xs
in merge (combine a) (combine b)
where split [] = ([],[])
split [a] = ([a],[])
split (x:y:cs) = let (a,b) = split cs in (x:a,y:b)
merge a [] = a
merge [] b = b
merge a@(x:a') b@(y:b') = case compare x y of
EQ -> x : merge a' b'
LT -> x : merge a' b
GT -> y : merge a b'
{-# INLINE eachPerm #-}
eachPerm :: (forall x y z. (Perm x y z) => View x y z -> Hints) -> [ Rule ]
eachPerm rule =
[ (\ known -> rule ( view known :: View R C V ) )
, (\ known -> rule ( view known :: View R V C ) )
, (\ known -> rule ( view known :: View C V R ) )
, (\ known -> rule ( view known :: View C R V ) )
, (\ known -> rule ( view known :: View V R C ) )
, (\ known -> rule ( view known :: View V C R ) )
]
ruleSubsets :: [ Rule ]
ruleSubsets = eachPerm ruleSubsetP ++ [ (\ known -> ruleSubsetP (view known :: View B D V) ) ]
allRules :: [ Rule ]
allRules = ruleSubsets ++ ruleBlocks
-- Applying the rules
{- The evolution strategy is simple:
Apply each rule in turn, keeping track whether or not there were any changes.
If all rules cause no change then it is done evolving.
-}
evolve :: Cube -> Cube
evolve cube =
let (cube',changed) = foldl step (cube,False) allRules
in if changed then evolve cube' else cube'
where step orig@(known,_) rule =
case rule known of
[] -> orig
hints -> (known // hints,True)
toCube :: (Monad m,Enum e) => [(e,e,e)] -> m Cube
toCube locs = do hints <- liftM concat $ mapM setLoc locs
return (emptyCube // hints)
where setLoc:: (Enum e,Monad m) => (e,e,e) -> m Hints
setLoc i@(re,ce,ve) = mapM checkM (rs ++ cs ++ vs)
where r = toEnum $ fromEnum re
c = toEnum $ fromEnum ce
v = toEnum $ fromEnum ve
rs = [((r',c,v),off) | r' <- fullRange, r' /= r]
cs = [((r,c',v),off) | c' <- fullRange, c' /= c]
vs = [((r,c,v'),off) | v' <- fullRange, v' /= v]
checkM hint = if check hint then return hint
else fail "Input location is out of range"
check ((r,c,v),_) = and [ inRange (lo,hi) r
, inRange (lo,hi) c
, inRange (lo,hi) v ]
fromCube :: (Enum e) => Cube -> [(e,e,e)]
fromCube cube = map head . filter exactlyOne $
[ [ eee | v <- fullRange, cube ! (r,c,v)
, let eee = (toEnum $ fromEnum r
,toEnum $ fromEnum c
,toEnum $ fromEnum v) ]
| r <- fullRange, c <- fullRange ]
consistent :: Cube -> Bool
consistent known = and [ consistentView ( view known :: View R C V )
, consistentView ( view known :: View C V R )
, consistentView ( view known :: View R V C )
, consistentView ( view known :: View V B D )
]
where
consistentView :: (Perm a b c) => View a b c -> Bool
consistentView view = and [ atLeastOne [ () | c <- fullRange, view ! (a,b,c) ]
| a <- fullRange, b <- fullRange ]
checkCube :: (Monad m) => String -> Cube -> m Cube
checkCube msg cube = if consistent cube then return cube else fail msg
deduce :: (Monad m,Enum e) => [(e,e,e)] -> m [(e,e,e)]
deduce locs = toCube locs >>=
checkCube "Inconsistent locations passed in" >>=
return . evolve >>=
checkCube "Inconsistent cube deduced from input" >>=
return . fromCube
test :: [(E,E,E)]
test = [(1,8,1),(1,9,2),(2,5,3),(2,6,5),(3,4,6),(3,8,7),(4,1,7),(4,7,3),(5,4,4),(5,7,8),(6,1,1),(7,4,1),(7,5,2),(8,2,8),(8,8,4),(9,2,5),(9,7,6)]
testC :: IO Cube
testC = toCube test
check :: IO ()
check = do c <- testC
print (consistent (evolve c))
