HaskellとSATソルバーで数独を解く
「システム検証基礎演習」の課題で、SATソルバーを使って数独の解答を出力せよというものがあったのでそのレポートを載せてみます。
どう考えても以下のページの解説の方がわかりやすいし実装もスマートだと思うけれど。
http://d.hatena.ne.jp/ku-ma-me/20080108/p1
- -
SATソルバーを用いて、数独を解いた。SATソルバーに与える連言標準形はHaskellプログラムで作成し、またSATソルバーの出力した一つの解を別のHaskellプログラムで目に見える形とした。
■原理
命題変数の値はtrueかfalseだけである。従って、1から9までの値を持つ数独の一つのマスを表すには、9つの命題変数を使うのがスマートであろう。命題変数p1からp9までがある一つのマスの数字を表し、もしこのマスに3が含まれるのであれば、p3がtrueで他はfalseになるとする。
数独のルールから、全てのマスには何らかの数字を埋めなければならないから、
p1 ∨ p2 ∨ p3 ∨ p4 ∨ p5 ∨ p6 ∨ p7 ∨ p8 ∨ p9
が成り立つ。一方別のマスがq1からq9までの命題変数で表され、やはりこのマスにも何らかの数字を埋めなければならないから、
q1 ∨ q2 ∨ q3 ∨ q4 ∨ q5 ∨ q6 ∨ q7 ∨ q8 ∨ q9
も成り立つ。全てのマスに対してこのような論理式を作る。
次に、同じ行、同じ列、同じブロックには、全て違った数字が選ばれるという制限を論理式で書く。なお、この制限を加えることにより、一つのマスを表す9つの命題変数の中の複数が同時にtrueとなることは自動的になくなる。一つのマスが状態の揺らぎを持ってしまうと、その列（行、ブロック）の他のマスの数字をそれぞれ重ならないように選択することができなくなってしまうからである。これは鳩ノ巣原理で証明できる。
さて、x列y行のマスを、p(x,y,1)からp(x,y,9)の命題変数で表すとする。x列y行のマスが1であった場合、その行の他のマスが1ではないのであれば、次の節たちが全て満たされなければならない。
p(x,y,1) → ¬(1,y,1) p(x,y,1) → ¬(2,y,1) ： p(x,y,1) → ¬(x-1,y,1) p(x,y,1) → ¬(x+1,y,1) ： p(x,y,1) → ¬(8,y,1) p(x,y,1) → ¬(9,y,1)
or節をandで繋げる連言標準形に直すと、
（¬p(x,y,1) ∨ ¬p(1,y,1) ） ∧ （¬p(x,y,1) ∨ ¬p(2,y,1) ） ： ∧ （¬p(x,y,1) ∨ ¬p(x-1,y,1) ） ∧ （¬p(x,y,1) ∨ ¬p(x+1,y,1) ） ： ∧ （¬p(x,y,1) ∨ ¬p(8,y,1) ） ∧ （¬p(x,y,1) ∨ ¬p(9,y,1) ）
となる。このような節を、1だけではなく2から9までの数字に対しても作る。また、その列の他のマスと数字が重ならないことは次のように書ける。
（¬p(x,y,1) ∨ ¬p(x,1,1) ） ∧ （¬p(x,y,1) ∨ ¬p(x,2,1) ） ： ∧ （¬p(x,y,1) ∨ ¬p(x,y-1,1) ） ∧ （¬p(x,y,1) ∨ ¬p(x,y+1,1) ） ： ∧ （¬p(x,y,1) ∨ ¬p(x,8,1) ） ∧ （¬p(x,y,1) ∨ ¬p(x,9,1) ）
これもやはり1だけではなく、2から9までの数字に対しても同様に作る。更にブロック内の他のマスと数字が重ならないことを次のように書く。尚これは2列3行のマスの例である。
（¬p(2,3,1) ∨ ¬p(1,1,1) ） ∧ （¬p(2,3,1) ∨ ¬p(2,1,1) ） ∧ （¬p(2,3,1) ∨ ¬p(3,1,1) ） ∧ （¬p(2,3,1) ∨ ¬p(1,2,1) ） ∧ （¬p(2,3,1) ∨ ¬p(2,2,1) ） ∧ （¬p(2,3,1) ∨ ¬p(3,2,1) ） ∧ （¬p(2,3,1) ∨ ¬p(1,3,1) ） ∧ （¬p(2,3,1) ∨ ¬p(3,3,1) ）
当然ながら、これも2から9までの数字に対しても同様に作る。
ここまでの、同じ行、列、ブロック内に違った数字が表れないという制限は、一つのマスに対してのものだから、更にこれらの節を全てのマスに対しても作らなければならない。
これで数独のルールは記述できたと言える。残るは固有の「問題」に拠る制限である。これは非常に単純で、単に命題変数だけの節を作ればいい。例えば2列3行のマスに既に6が埋まっていたのであれば、
p(2,3,6)
これだけである。埋まっているマスのぶんだけ、このような節を作る。
上記より必要な節の数は、
9*9 + 8*9*3*9*9 + 既に埋まっているマスの数 = 17577 + α
となる。
■実装
SATソルバーに与える連言標準形を作るHaskellプログラムは以下である。
-- solving sudoku for SAT solver import List problem = [5,3,0,0,7,0,0,0,0, 6,0,0,1,9,5,0,0,0, 0,9,8,0,0,0,0,6,0, 8,0,0,0,6,0,0,0,3, 4,0,0,8,0,3,0,0,1, 7,0,0,0,2,0,0,0,6, 0,6,0,0,0,0,2,8,0, 0,0,0,4,1,9,0,0,5, 0,0,0,0,8,0,0,7,9] index (_, _, 0) = "0\n" index (x, y, n) | n > 0 = (show $ (y-1) * 9 * 9 + (x-1) * 9 + (n-1) + 1) ++ " " | n < 0 = "-" ++ (show $ (y-1) * 9 * 9 + (x-1) * 9 + (-n-1) + 1) ++ " " allGrid = [(x,y) | x <- [1..9], y <- [1..9] ] main = do let a = makeProblem problem ++ g ++ concatMap (\(x,y) -> concat $ hrzn(x,y) ++ vtcl(x,y) ++ box(x,y)) allGrid putStrLn $ "p cnf " ++ (show $ 9*9*9) ++ " " ++ (show $ cc a) putStr $ concatMap index $ a -- make problem restrict makeProblem l = mp 1 1 l where next x y l | x == 9 = mp 1 (y+1) l | otherwise = mp (x+1) y l mp x y (n:l) | n>0 = (x, y, n):((0,0,0):(next x y l)) | otherwise = next x y l mp _ _ _ = [] -- input number to all grids g = concatMap g' allGrid where g' (x, y) = g'' 9 where g'' 0 = [(x, y, 0)] g'' n = (x, y, n):(g'' $ n-1) -- restrict number for sudoku rule hrzn (x, y) = [[(x, y, -n), (x, y', -n), (0,0,0)] | n <- [1..9], y' <- [1..9]\\[y]] vtcl (x, y) = [[(x, y, -n), (x', y, -n), (0,0,0)] | n <- [1..9], x' <- [1..9]\\[x]] box (x, y) = [[(x, y, -n), (x', y', -n), (0,0,0)] | n <- [1..9], x' <- [(boxS x)..(boxS x + 2)], y' <- [(boxS y)..(boxS y + 2)], (x/=x' || y/=y')] where boxS n = n - ((n-1) `mod` 3) -- count clause cc ((_, _, 0):l) = 1 + (cc l) cc (_:l) = cc l cc _ = 0
このプログラムでやっていることは、ほぼ原理で説明したことそのままである。内部的に命題変数(列,行,数字)を表す(Int, Int, Int)の型を用いているが、これはややトリッキーな使い方をされており、最後の数字が負になっていれば、否定であることを表す。さらに0であれば、節の終端であることを示す。このようなデータを、関数indexにおいて、SATソルバーが扱う自然数のみで識別される命題変数へと変換している。
さて、このプログラムの出力（ここれはファイル名はsudokuとする）をMiniSatに与えてみた。
> ./MiniSat sudoku sudoku.out （略） SATISFIABLE
充足されたことがわかる。さて解はsudoku.outを見ればいいわけだが、命題変数が全て自然数になっており非常にわかりにくい。そこで、次のようなHaskellプログラムを用意した。
import Text.ParserCombinators.Parsec main = do src <- getContents putStr $ concatMap indexr $ readResult src readResult src = case parse (do string "SAT" many space num <- many readVar return num) [] src of Right madeSrc -> madeSrc Left err -> error $ show err readVar :: Parser Int readVar = do string "-" n <- many1 digit many1 space return $ negate $ read n <|> do n <- many1 digit many1 space return $ read n indexr m | m > 0 = let y = ((m-1) `div` (9 * 9)) + 1 m' = ((m-1) `mod` (9 * 9)) x = (m' `div` 9) + 1 n = (m' `mod` 9) + 1 in show n ++ if x==9 then "\n" else [] indexr _ = []
これは単純に、関数indexの逆関数を適用するプログラムである。このプログラムにsudoku.outを与えてみると、次のようになる。
>main < sudoku.out 534678912 672195348 198342567 859761423 426853791 713924856 961537284 287419635 345286179
これは我々一般人の理解できる数独の解答である。コンピュータが本問題に勝利したことが明らかとなった。
- -
Haskellは良い言語だけれど、やっぱりIOが面倒ですね。SATソルバーの出力を読み込むのに、わざわざParsecを使う方がラクだ、と思ってしまうほどに。
そういえば8クイーンの解答を出力させるプログラムも書いたので、それも載せてみます。数独のプログラムも微妙だけれど、こちらはヨリひどいプログラムです。
data XY = Horizon | Vertical deriving Eq main = do putStr "p cnf 64 1472\n" putStr $ (p Horizon 8) ++ (p Vertical 8) ++ (concat.concat $ map restrict (allPoint 8 8)) index x y = show $ (x-1) * 8 + y p _ 0 = "" p xy n = (p' 8) ++ (p xy $ n - 1) where p' 0 = "0\n" p' m = (if xy==Horizon then index n m else index m n) ++ " " ++ (p' $ m - 1) allPoint x y = [(x', y') | x' <- [1..x], y' <- [1..y]] restrict (x, y) = restrictX(x, y) ++ restrictY(x, y) ++ restrictN(x, y) where restrictX (x, y) = ["-" ++ (index x y) ++ " -" ++ (index x' y) ++ " 0\n"| x' <- [1..8], x'/=x ] restrictY (x, y) = ["-" ++ (index x y) ++ " -" ++ (index x y') ++ " 0\n"| y' <- [1..8], y'/=y ] restrictN (x, y) = ["-" ++ (index x y) ++ " -" ++ (index x' y') ++ " 0\n"| x' <- [1..8], x'/=x, y' <- [1..8], y'/=y, ((x'-x) == (y'-y) || (x'-x) == (y-y') )]
