Project Euler の96問目が数独を解く問題で、手筋を調べて解くプログラムを書こうかと思ったが、SATソルバを使うとゲームのルールごとに解法を考えてプログラムを作らなくても、制約を与えれば解を求めることができる、ということを知ったので使って見ることにした。

参考にしたのはHaskellとSATソルバーで数独を解く - メタサブカル病のページ。このページのプログラムを実行すると、数独を解くための制約を標準出力に出力する。それをMiniSatというオープンソースのSATソルバに与えて、その出力をパースして結果をわかりやすく9x9に表示している。Haskellだけで閉じてるわけじゃなくて、SATソルバは外部プログラムになっている。

MiniSatに与える制約は、乗法標準形(Conjunctive Normal Form、orでつないだ命題を、andでつないだもの)という形式で与える必要がある。論理変数を正の数で表して、そのまま正の数を使えばその論理変数が真の時、符号反転して負の数を使えば偽を表す。プログラムが出力する制約は次のようなファイルになる:

p cnf 729 17607
5 0
...
9 8 7 6 5 4 3 2 1 0
90 89 88 87 86 85 84 83 82 0
171 170 169 168 167 166 165 164 163 0
...

最初の行で変数の数と節の数を明示していて、続く各行がorでつないだ命題で、全体をandでつないだものが制約となる(各行最後の0は行の終わりを示しているだけで、意味はない)。1つの論理変数は真か偽かしか持てないので、数独では各マスに対して変数を9個用意して、変数1が真なら1が入る、変数2が真なら2が入る、…というように用いている。なので9x9マスx9変数 = 729個変数がある。節は17577+問題に埋め込まれている数字の個数となる。

で、外部プログラムを使わずにHaskellでSATを解くライブラリがないか調べたところいろいろあるらしく、どれがいいのかよくわからなかった。gatlin / sat.hs - gistはファイル1つだけでソースも短いので、これを試してみたところ、速いかどうかは他と比べてないのでわからないが、ちゃんと動いて問題が解けた。うーむ便利。

関数solveに与える引数はMiniSatとほぼ同じで、整数の二重リストで、外側のリストが節全体、内部の各リストが節、数値が論理変数の真(正)または偽(負)となっている。戻り値はMaybeで、Nothingなら解なし、あれば論理変数の値のリストで、ソートはされてない。

このSATソルバを使って、Haskellだけで解けるコード:

-- This code requires SAT solver module,
-- e.g. https://gist.github.com/gatlin/1755736

import Data.List (sort)
import Sat (solve)

main = do
  case solveSudoku problem of
    Nothing -> putStrLn "Impossible"
    Just answer -> putStr $ unlines $ map (concatMap show) $ pack 9 answer

pack n [] = []
pack n xs = (take n xs) : pack n (drop n xs)

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]

solveSudoku :: [Integer] -> Maybe [Integer]
solveSudoku field = solve ts >>= return . toAnswer
  where ts = map (map cellToVar) $ makeProblem field ++ sudokuConstraints
        toAnswer = map varToNum . sort . filter (> 0)
        varToNum v = (v - 1) `mod` 9 + 1
        cellToVar (x,y,n) | n > 0 = var x y n
                          | n < 0 = -(var x y (-n))
        var x y n = (y-1) * 9 * 9 + (x-1) * 9 + (n-1) + 1

makeProblem :: [Integer] -> [[(Integer, Integer, Integer)]]
makeProblem l = map mp $ filter ((/= 0) . snd) $ zip allGrid l
  where mp ((x,y), n) = [(x, y, n)]

sudokuConstraints :: [[(Integer, Integer, Integer)]]
sudokuConstraints = (cellConstraints ++ horzConstraint ++ vertConstraint ++
                     boxConstraint)
  where horzConstraint = concatMap hrzn allGrid
        vertConstraint = concatMap vtcl allGrid
        boxConstraint = concatMap box allGrid

cellConstraints :: [[(Integer, Integer, Integer)]]
cellConstraints = map f allGrid
  where f (x,y) = map (\n -> (x,y,n)) [1..9]

hrzn (x, y) = [[(x, y, -n), (x, y', -n)] | n <- [1..9], y' <- [1..9], y' /= y]
vtcl (x, y) = [[(x, y, -n), (x', y, -n)] | n <- [1..9], x' <- [1..9], x' /= x]
box (x, y) = [[(x, y, -n), (x', y', -n)] | 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)

allGrid = [(x,y) | y <- [1..9], x <- [1..9]]

各列や行に同じ数字が入らないという制約は、論理式「P→Q」が「¬P ∨ Q」と表せることを使って、例えば論理式Aを「マスaに1が入る」、論理式Bを「マスbに1が入る」とすると、「A→¬B(AならばBじゃない)」を「¬A ∨ ¬B」という形で構築している。

ソルバの結果から数独の解答を取り出すのは簡単で、論理変数が各マスにつき9個x横方向に9マスx縦方向に9マス、1から729まで順番に割り振られているので、真となっている変数(正の値)だけ取り出してソートして、9で割ったあまりが左上のマスから順に入る答えの数になる(1オリジンなところだけ注意)。

5秒くらいかかることはあるにせよ、変数729個、制約17,600節を、現実的な時間で解けるってすごいな。パズルのルールから手筋を考えずに、一般的に解けるって便利すぎる。