いろんな方法で数独解いてみる(SAT, CSP編) #Python

数独を色々な論理プログラミングで解いてみて、その一長一短なところをみてみましょう。
紹介するのは次の5つです

SAT - CNF(連言標準形)の論理式をとく
CSP - 変数の定義域と制約から解を求める
ASP - 一階述語論理式をとく(次回)
Prolog - ASPと似てるやつ(次回)
自前プログラム - 比較用、rust(次回)

数独のルールの定義

ルールを定義して、プログラムが定義に乗っ取っているかを確認します

一つのマスに数字は一つ
一つの行に同じ数字は一つ
一つの列に同じ数字は一つ
一つのブロックに同じ数字は一つ
入力を満たす

入力は、「セルの行x列yに値iが入る」という条件をタプル(x,y,i)で表し、入力全体はそのタプルのリストで表されるとします。

SAT

SATソルバーはCNF(連言標準形)と呼ばれる論理式を充足する真偽値が存在するかを求めるアルゴリズムです。CNFは節と呼ばれるリテラルの論理和($x_1 \vee \neg x_2\dots \vee x_n$など)の集合の論理積で表されるものです。

それでは、ルールをSATに変換していきましょう。

定義

「セルの行x列yに値iが入る」とき$g_{x,y,i}$が真であるとします。逆に「セルの行x列yに値iが入らない(他の値が入る)」というとき$g_{x,y,i}$は偽となります。

ルール1(一つのマスに数字は一つ)

安直に考えるとこんな感じになります

(g_{x,y,1} \wedge \neg g_{x,y,2} \wedge \dots \wedge\neg g_{x,y,9}) \vee(\neg g_{x,y,1} \wedge g_{x,y,2} \wedge \dots \wedge\neg g_{x,y,9}) \vee \dots \vee\\
(\neg g_{x,y,1} \wedge \neg g_{x,y,2} \wedge \dots\wedge g_{x,y,9})

でもこれはCNFの形ではないので、分配法則(?)のように展開して簡略化すると次のようになります

(g_{x,y,1} \vee g_{x,y,2} \vee \dots g_{x,y,9}) \wedge (\neg g_{x,y,1} \vee \neg g_{x,y,2})\wedge (\neg g_{x,y,1} \vee \neg g_{x,y,3}) \wedge \dots \wedge\\
(\neg g_{x,y,8} \vee \neg g_{x,y,9})

ルール2,3,4(一つの行(列、ブロック)に同じ数字は一つ)

「行xに数字iが一つ以上存在する」という条件は次のように書けます

g_{x,1,i} \vee g_{x,2,i} \vee \dots g_{x,9,i}

これに「行xに数字iが2つ以上存在しない」という条件を加えればルール2になるわけですが、ルール1と組み合わせるとその条件は必要無くなります。なぜなら、数字1~9が9個のセルにそれぞれ1つ以上あるわけですから、鳩の巣原理的に、それぞれの数字は1つしか存在しえないわけです。

ちなみに、「行xに数字iが2つ以上存在しない」という条件を表すとこうなります。

(\neg g_{x,1,i} \vee \neg g_{x,2,i})\wedge (\neg g_{x,1,i} \vee \neg g_{x,3,i}) \wedge \dots \wedge(\neg g_{x,8,i} \vee \neg g_{x,9,i})

列やブロックについては同様であるので、省略

ルール5(入力を満たす)

「セルの行x列yに値iが入る」という条件は$g_{x,y,i}$に真を割り当てるということになります。

実装

ということでpysatのminisatを使って解きます。

from pysat.solvers import Minisat22
from itertools import product, combinations


def grid(i, j, k):
  return i * 81 + j * 9 + k + 1

def sudoku_sat(arr):
  m = Minisat22()

  # ルール１
  for i, j in product(range(9), range(9)):
    m.add_clause([grid(i, j, k) for k in range(9)])
    for k1, k2 in combinations(range(9), 2):
      m.add_clause([-grid(i, j, k1), -grid(i, j, k2)])

  # ルール2,3
  for i in range(9):
    for k in range(9):
      m.add_clause([grid(i, j, k) for j in range(9)])

  for j in range(9):
    for k in range(9):
      m.add_clause([grid(i, j, k) for i in range(9)])

  # ルール4
  for p, q in product(range(3), range(3)):
    for k in range(9):
      m.add_clause([grid(i, j, k) for i, j in product(range(p*3, p*3+3), range(q*3, q*3+3))])

  # ルール5
  for a in arr:
    m.add_clause([grid(a[0], a[1], a[2] - 1)])
  if not m.solve():
    return None

  model = m.get_model()
  return [
    [
      [k + 1 for k in range(9) if model[grid(i, j, k) - 1] > 0][0]
      for j in range(9)
    ]
    for i in range(9)
  ]

CSP

~~大学の課題のコピペです~~

CSP(制約プログラミング)は次の$V, D, C$からなる問題を解くアルゴリズムです。

$V= \{v_0, v_1,\dots,v_n\}$ - 変数の集合
$D$ - それぞれの変数の定義域。$D_i$は変数$v_i$の定義域を示す
$C$ - 変数同士のの制約(例: $v_1 \neq v_2$)の集合

前提

変数$v_{x,y}$が行x列yのセルの値であるとします。つまり:

$V = \{v_{0,0}, v_{0,1},\dots,v_{8,8}\}$
$D_{x,y} = \{1,2,\dots,9\}$

この前提のみでルール１は満たすことができます。

ルール2,3,4(一つの行(列、ブロック)に同じ数字は一つ)

ルールを書き下すと次のようになります

(v_{x,0}, v_{x,1})\in\{(1,2), (1,3), \dots (9,8)\} \\
(v_{x,0}, v_{x,2})\in\{(1,2), (1,3), \dots (9,8)\} \\
\vdots\\
(v_{x,7}, v_{x,8})\in\{(1,2), (1,3), \dots (9,8)\}

ですが、通常CSPにはAllDifferentという便利機能がありますので、それを利用すると簡単に記述できます

\text{AllDifferent}(v_{x,0}, v_{x,1},\dots,v_{x,8})

列、ブロックも同様です。

ルール5(入力を満たす)

ルール5は定義域から組み込むこともできますが、数独の根本のルールの部分と入力を分離したいという目的のため、制約で表すとします。「セルの行x列yに値iが入る」という制約は次のようになります

v_{x,y} \in \{i\}

実装

python-constraintを使って解きます。

from constraint import Problem, AllDifferentConstraint, InSetConstraint

def sudoku_csp(arr):

  def grid(i, j):
    return '{}_{}'.format(i, j)

  p = Problem()

  for i, j in product(range(9), range(9)):
    p.addVariable(grid(i, j), range(1, 10))

  for i in range(9):
    p.addConstraint(AllDifferentConstraint(), [grid(i, j) for j in range(9)])
    p.addConstraint(AllDifferentConstraint(), [grid(j, i) for j in range(9)])

  for m, n in product(range(3), range(3)):
    p.addConstraint(
      AllDifferentConstraint(),
      [grid(i, j) for i, j in product(range(m*3, m*3+3), range(n*3, n*3+3))]
    )

  for a in arr:
    p.addConstraint(InSetConstraint([a[2]]), [grid(a[0], a[1])])

  solution = p.getSolution()
  return [
    [
      solution[grid(i, j)]
      for j in range(9)
    ]
    for i in range(9)
  ]

最後に

この調子でアドベントカレンダーやったら死にそうなので、次から手を抜きます。