概要
SATソルバの実装minisatをpython経由で呼び出して数独を解きました。
数独の81個のマスそれぞれに(1～9という数に対応する)9個のフラグをつけた729bitを考えます。数独の解は729bitのうち81箇所が1で残りが0です。数独の各種ルールも各bit間に対する制約として表現できるのでSATソルバで解けるというわけです。
1.実装
1-1. マスとユニット
81個のマスに0～81の番号を振ります。左上が0番目、1段目は0～8を割り振って右下が80番目になります。ユニットとは1～9の数がぴったり入るマスの集合として定義します。たて・よこ・矩形の27個のユニットがあります。例えば2段目のよこのユニットは [9,10,11,12,13,14,15,16,17] という配列で表します。
1-2. 命題変数
各マスに対して、「nが入るか (n=1,2,...9)」を意味する命題変数を用意します。つまりマスは81個あるので81×9である729個の命題変数を用意します。$2^{729}$通りの真理値割り当ての方法から、解となる割り当てを見つけるという具合に数独を充足可能性問題に落とし込むというわけです。
i番目のマスに数nが入るか？を表す命題変数を i*9+(n-1)+1 番目と符号化すれば、命題変数は 1～729 という識別番号になります。後述するように制約の記述には非ゼロの自然数とする必要があるのでこの方法が自然かと思います。
1-3. 制約A
解であるための条件の一つは「全てのマスの値が1～9のどれかただひとつに定まっている」です。「ただひとつ」とは、「少なくとも一つ」かつ「任意の二つの数m,nに対応する命題変数が同時にTrueにならない」と分解できます。前者は9つの命題変数の選言がTrueになると言い換えられます。後者はドモルガンの法則を使えば「任意の二つの数m,nに対応する命題変数の否定の選言がFalseになる」と節で表現できます（コードコメントの方がわかりやすいかも）。
1-4. 制約B
解であるためのもう一つの条件は「任意のユニット uについて、任意の数 mについて、m が u の中にただ一つ配置されている」と表現できます。制約Aと同様に「ただ一つ」というのを二つの制約に分解してあげればよいです。
1-5. 制約C
制約A,Bは任意の数独の問題に適用できる制約です。特定の問題を解く際は、これらに加え初期盤面を制約の形にします。これは明快です。
1-6. MiniSat
Ubuntuでは sudo apt install minisat で /usr/bin/minisat というバイナリがインストールされます。SATソルバです。入出力ファイルを明示して呼び出すことで、充足解を得られます:
/usr/bin/minisat input.txt output.txt
入出力フォーマットは、ひと昔前に記述したので興味ある人はどうぞ。コードではpythonの os.system() を使うことでプログラムの中から呼び出しをしています。
1-7. 解く問題
世界一難しい数独 を解きます。
この問題では、節の数は11988個でした。
2. コード
python3で実行します。カレントディレクトリに input.txt と output.txt が生成されます。コンソールに seems good とでれば動作しているはずです。/usr/bin/minisat がインストールされていることを前提にしています:
# coding:utf8
import os
# ユニット（計27)
UNITS = []
for i in range(9):
UNITS.append([i+9*j for j in range(9)]) # タテ
UNITS.append([9*i+j for j in range(9)]) # ヨコ
for i in [0,3,6,27,30,33,54,57,60]:
UNITS.append([i+j for j in [0,1,2,9,10,11,18,19,20]]) # 四角
# 命題変数 (729=81*9個)
def varOf(i,m):
"""マスiに数mが入るかを格納する命題変数番号(1,2,...)を返す"""
return i*9+m+1
def invVar(n):
return ((n-1)//9, (n-1)%9)
CONST = []
# 制約A: 各セルに注目
for i in range(81):
# どれか一つの値が入っているべき
CONST.append([varOf(i,m) for m in range(9)])
# 任意の二つの数m,nに対して同時にフラグが立たない
# ¬(var(i,m)⋀var(i,n)) 変形して ¬var(i,m)⋁¬var(i,n)
for m in range(9):
for n in range(m+1,9):
CONST.append([-varOf(i,m), -varOf(i,n)])
# 制約B: 各ユニットに注目
for u in UNITS:
for m in range(9):
# 数mがユニットに含まれているべき
CONST.append([varOf(i,m) for i in u])
# 任意の二つのマスi,jが同じ値をもたない すなわち
# ¬(var(m,i)⋀var(m,j)) 変形して ¬var(m,i)⋁¬var(m,j)
for i in range(9):
for j in (range(i+1,9)):
CONST.append([-varOf(i,m), -varOf(j,m)])
# 制約C: 初期盤面
Q = '85...24..72......9..4.........1.7..23.5...9...4...........8..7..17..........36.4.'
for i, m in enumerate(Q):
if m not in "123456789":
continue
CONST.append([varOf(i, int(m)-1)]) # -1 to 0-indexed
# ソルバ呼び出し
with open("input.txt", "w", encoding="utf8") as ofh:
ofh.write("p cnf {} {}\n".format(729, len(CONST)))
for cst in CONST:
ofh.write(" ".join([str(c) for c in cst] + ["0"]))
ofh.write("\n")
os.system("/usr/bin/minisat input.txt output.txt")
with open("output.txt") as ifh:
if ifh.readline() != "SAT\n":
raise Exception("ouch")
ans = ifh.readline()
data = [0]*81
for x in ans.split(" ")[:-1]: # 最後の1つは区切り文字
n = int(x)
if n < 0:
continue
i, m = invVar(n)
data[i] = m+1
ANS="859612437723854169164379528986147352375268914241593786432981675617425893598736241"
if "".join(map(str, data)) == ANS:
print("seems good.")