数独をSATで解く ― 「ちょうど1つ」をCNFで表す#
この章で学ぶこと#
数独を、SAT（充足可能性問題）ソルバーで解きます。問題をブール式に翻訳する流れを、小さな推理パズルから始めて、数独まで段階的に進めます。
SAT と CNF（連言標準形）とは何か
選言標準形（DNF）を CNF に変換する方法（SymPy と手作り）
「グループの中で1つだけが True」という制約の書き方
pysatで数独を解き、assumptionsで初期盤面を与える
充足可能性問題（SAT）#
充足可能性問題（SATisfiability problem、SAT）は、「与えられたブール式を真にする変数の値の組み合わせはあるか？」という問題です。
SAT は最初にNP完全であることが証明された問題で、計算機科学の基礎理論、アルゴリズム、人工知能、ハードウェア設計など、多くの分野に関係しています。そして、実用的な SAT ソルバーは非常に高速なので、「解きたい問題を SAT に変換する」という手法が広く使われています。
この章では、pysat で数独を解きます。まず小さな推理パズルで、ブール式と CNF の書き方を確認しましょう。
推理パズル
4人の容疑者が、次のように供述しました。
甲：「私が犯人ではない」
乙：「丁が犯人だ」
丙：「乙が犯人だ」
丁：「乙は私を陥れようとしている」
4人のうち、真実を話しているのは1人だけです。犯人は誰でしょう。
甲・乙・丙・丁について、「その人が犯人である」を表すブール変数を A, B, C, D とします。~A は「甲は犯人ではない」です。
各供述は、次のようなブール式になります。
容疑者 |
供述 |
ブール式 |
|---|---|---|
甲 |
私が犯人ではない |
S1 = ~A |
乙 |
丁が犯人だ |
S2 = D |
丙 |
乙が犯人だ |
S3 = B |
丁 |
乙は私を陥れようとしている |
S4 = ~D |
真実を話している人が1人だけなので、S1〜S4 のうちちょうど1つが真になります。
「ちょうど1つが真」を、まず素直に書いてみます（& は論理積、| は論理和）。
S1 & ~S2 & ~S3 & ~S4 |
~S1 & S2 & ~S3 & ~S4 |
~S1 & ~S2 & S3 & ~S4 |
~S1 & ~S2 & ~S3 & S4
これは「OR でつないだ AND の集まり」で、選言標準形（DNF, Disjunctive Normal Form）と呼びます。
SAT ソルバーが受け取れるのは連言標準形（CNF, Conjunctive Normal Form）だけです。CNF は「OR でつないだ節を、AND で並べたもの」です。そこで、「ちょうど1つが真」を別の形で書き直します。
2つ以上が真になってはいけない →
~(Si & Sj)全部が偽でもいけない →
~( ~S1 & ~S2 & ~S3 & ~S4 )
~(S1 & S2) & ~(S1 & S3) & ~(S1 & S4) &
~(S2 & S3) & ~(S2 & S4) & ~(S3 & S4) &
~(~S1 & ~S2 & ~S3 & ~S4)
ド・モルガンの法則（~(X & Y) = ~X | ~Y）で変形すると、CNF になります。
変換結果（S1 = ~A, S2 = D, S3 = B, S4 = ~D を代入）
~S1 | ~S2 → A | ~D
~S1 | ~S3 → A | ~B
~S1 | ~S4 → A | D
~S2 | ~S3 → ~D | ~B
~S2 | ~S4 → ~D | D ← 常に真（何も制約しない）
~S3 | ~S4 → ~B | D
S1 | S2 | S3 | S4 → ~A | D | B | ~D
~D | D は常に真なので、この節は何も制約しません。それでも害はないので、そのまま入れています。
SymPy で DNF から CNF へ#
手で変換するのは間違えやすいので、SymPy に任せます。
from sympy import symbols
from sympy.logic.boolalg import to_cnf
A, B, C, D = symbols("A:D")
# 各人の供述
S1 = ~A
S2 = D
S3 = B
S4 = ~D
# ちょうど1人だけが真実を話す（DNF）
dnf = ((S1 & ~S2 & ~S3 & ~S4) |
(~S1 & S2 & ~S3 & ~S4) |
(~S1 & ~S2 & S3 & ~S4) |
(~S1 & ~S2 & ~S3 & S4))
# simplify=True：式を簡単にしてから返す
print("CNF:", to_cnf(dnf, simplify=True))
# simplify=False：分配法則でそのまま展開する
cnf_raw = to_cnf(dnf, simplify=False)
print("simplify=False だと", len(cnf_raw.args), "個の節になる")
CNF: A & ~B
simplify=False だと 41 個の節になる
simplify=True にすると、式を簡単にしてから返してくれます。結果は A & ~B、つまり「甲が犯人で、乙は犯人ではない」です。
simplify=False にすると、式を分配法則でそのまま展開するので、節の数が何十個にも膨らみます。A | ~A のような常に真の節も消されません。数独のような大きな問題では、この膨らみがそのまま解く時間になります。
どちらの CNF でも、答えは同じになるはずです。
from sympy.logic.inference import satisfiable
print(satisfiable(to_cnf(dnf, simplify=True)))
print(satisfiable(cnf_raw))
{A: True, B: False}
{A: True, D: False, B: False}
satisfiable() は、式を真にする変数の割り当てを返します。A が真なので、甲が犯人です。C は式に出てこないので、どちらでもよい（制約なし）になっています。
pysat で解く#
CNF は、入れ子のリストで表せます。整数1つがブール変数1つに対応し、正の数はその変数が真、負の数は偽（否定）を意味します。たとえば 1 は A、-1 は ~A です。
上の CNF を、そのまま数値のリストにします（A=1, B=2, C=3, D=4）。
append_formula()で CNF をソルバーに追加する（何度呼んでもよい）solve()で解を求めるget_model()で解を取り出す（正の数が真）
from pysat.solvers import Solver
sat = Solver()
problem = [[1, -4], [1, -2], [1, 4], # ~S1 | ~S2 など
[-4, -2], [-4, 4], [-2, 4], # ~D | D は常に真
[-1, 4, 2, -4]] # S1 | S2 | S3 | S4
sat.append_formula(problem)
print(sat.solve())
print(sat.get_model())
True
[1, -2, -3, -4]
[1, -2, -3, -4] は「A が真、B・C・D は偽」です。SymPy と同じ結論になりました。
数独#
数独は、9×9 のマスに1〜9の数字を入れ、次の条件を満たすパズルです。
各行に1〜9が1回ずつ
各列に1〜9が1回ずつ
各3×3ブロックに1〜9が1回ずつ
黒い数字が問題、灰色の数字が答えです。
SAT の変数は True / False の2値なので、「マス (r, c) に数字 n が入る」を1つのブール変数にします。3次元配列 bools を使います。
第0軸 … 行
第1軸 … 列
第2軸 … 数字（1〜9）
たとえば bools[4, 1, 3] が真なら、「5行目・2列目に数字4が入る」という意味です。数独は1始まりですが、配列の添字は0始まりであることに注意してください。
bools には、SAT の変数番号を入れます（変数番号は1始まり）。
import numpy as np
# 9×9×9 = 729 個の変数（1〜729）
bools = np.arange(1, 9 * 9 * 9 + 1).reshape(9, 9, 9)
数独の制約は、4つに整理できます。
# |
制約 |
グループの作り方 |
|---|---|---|
1 |
各マスには1つの数字だけ |
数字の軸（第2軸）をグループにする |
2 |
各行に同じ数字は1つだけ |
行と数字でグループにする |
3 |
各列に同じ数字は1つだけ |
列と数字でグループにする |
4 |
各ブロックに同じ数字は1つだけ |
ブロックと数字でグループにする |
どれも「あるグループの中で、ちょうど1つだけが True」という形です。
「ちょうど1つ」をCNFで表す#
変数 A, B, C のうち1つだけが真、という制約は、2つの条件に分けられます。
少なくとも1つは真: \(A | B | C\)
2つ以上が真にならない: \(\neg(A \& B) \land \neg(A \& C) \land \neg(B \& C)\)
ド・モルガンの法則で、2番目は \( (\neg A | \neg B) \land (\neg A | \neg C) \land (\neg B | \neg C) \) になります。
要素が N 個なら、「すべてを OR でつないだ節」1つと、「2つの組すべてについて ~X | ~Y」が \(N(N-1)/2\) 個です。2つの組は itertools.combinations で作れます。
get_conditions(bools) は、bools の最終軸をグループとして、この制約を返します。
最終軸の要素数
nに対して、2つの組すべてを作る各グループを OR でつないだ節を作る（少なくとも1つ真）
各組について
~X | ~Yの節を作る（2つ同時に真にならない）
from itertools import combinations
def get_conditions(bools):
"""最終軸をグループとして、「ちょうど1つだけ True」の CNF を作る"""
n = bools.shape[-1]
pairs = np.array(list(combinations(range(n), 2)))
conditions = []
# 1. 少なくとも1つは True
conditions.extend(bools.reshape(-1, n).tolist())
# 2. 2つ同時に True にならない
conditions.extend((-bools[..., pairs].reshape(-1, 2)).tolist())
return conditions
小さな例で動きを確かめます。[[1, 2, 3], [4, 5, 6]] は「2つのグループがあり、それぞれ3つの変数を持つ」という意味です。
[1, 2, 3],[4, 5, 6]… それぞれのグループで、少なくとも1つは真[-1, -2],[-1, -3],[-2, -3]… 1〜3 のうち2つが同時に真にならない[-4, -5],[-4, -6],[-5, -6]… 4〜6 についても同じ
get_conditions(np.array([[1, 2, 3], [4, 5, 6]]))
[[1, 2, 3],
[4, 5, 6],
[-1, -2],
[-1, -3],
[-2, -3],
[-4, -5],
[-4, -6],
[-5, -6]]
数独の制約を作る#
get_conditions() は最終軸をグループにするので、行・列・ブロックの制約は、swapaxes() でグループにしたい軸を最後に移動してから渡します。
制約1（各マスに1つ）… そのまま（最終軸が数字）
制約2（各行に1つ）… 列と数字を入れ替える
制約3（各列に1つ）… 行と数字を入れ替える
c1 = get_conditions(bools) # 各マスには1つの数字だけ
c2 = get_conditions(np.swapaxes(bools, 1, 2)) # 各行で同じ数字は1つだけ
c3 = get_conditions(np.swapaxes(bools, 0, 2)) # 各列で同じ数字は1つだけ
制約4（ブロック）は少し複雑です。bools は (行, 列, 数字) なので、行と列を「ブロック番号」と「ブロック内の位置」に分解します。
reshape(3, 3, 3, 3, 9)…(ブロック行, 行内, ブロック列, 列内, 数字)swapaxes(1, 2)…(ブロック行, ブロック列, 行内, 列内, 数字)reshape(9, 9, 9)…(ブロック番号, ブロック内の9マス, 数字)swapaxes(1, 2)…(ブロック番号, 数字, ブロック内の9マス)← グループは最後の軸
これで「ブロックごと・数字ごとに、9マスのうち1つだけが真」という制約になります。
# ブロックごとに (ブロック, 9マス, 数字) の形にする
blocks = np.swapaxes(bools.reshape(3, 3, 3, 3, 9), 1, 2).reshape(9, 9, 9)
c4 = get_conditions(np.swapaxes(blocks, 1, 2)) # 各ブロックで同じ数字は1つだけ
conditions = c1 + c2 + c3 + c4
print(len(conditions), "個の節ができました")
11988 個の節ができました
空の盤面を解く#
まず、ヒントが何もない状態で解いてみます。数独の制約だけを満たす答えが1つ出てきます（ヒントがないので、答えは無数にあります）。
モデルは変数番号のリストなので、(9, 9, 9) に戻して、真になっている変数の数字（最終軸の添字 + 1）を取り出します。
def format_solution(model, shape=(9, 9, 9)):
"""モデル（変数番号のリスト）を 9×9 の数字の盤面にする"""
grid = np.array(model).reshape(shape)
return np.where(grid > 0)[2].reshape(9, 9) + 1
sat = Solver()
sat.append_formula(conditions)
print(sat.solve())
format_solution(sat.get_model())
True
array([[2, 4, 6, 1, 3, 5, 8, 7, 9],
[1, 3, 5, 7, 8, 9, 2, 4, 6],
[8, 7, 9, 2, 4, 6, 1, 3, 5],
[5, 6, 3, 4, 1, 2, 7, 9, 8],
[4, 1, 2, 9, 7, 8, 5, 6, 3],
[7, 9, 8, 5, 6, 3, 4, 1, 2],
[6, 5, 4, 3, 2, 1, 9, 8, 7],
[3, 2, 1, 8, 9, 7, 6, 5, 4],
[9, 8, 7, 6, 5, 4, 3, 2, 1]])
実際の問題を解く#
問題は文字列で与えます。0 が空きマスです。空きマスではないところは「その数字の変数が真」という制約を追加すれば、その数字に固定できます。
sudoku_str = """
000000185
007030000
000021400
800000020
003905600
050000004
004860000
000040300
931000000"""
puzzle = np.array([[int(x) for x in line]
for line in sudoku_str.strip().split()])
# ヒントのマスについて「その数字である」を制約にする
r, c = np.where(puzzle != 0)
v = puzzle[r, c] - 1
hints = [[int(x)] for x in bools[r, c, v]]
sat = Solver()
sat.append_formula(conditions + hints)
sat.solve()
solution = format_solution(sat.get_model())
solution
array([[3, 6, 2, 7, 9, 4, 1, 8, 5],
[4, 1, 7, 5, 3, 8, 2, 6, 9],
[5, 9, 8, 6, 2, 1, 4, 3, 7],
[8, 7, 9, 4, 1, 6, 5, 2, 3],
[2, 4, 3, 9, 7, 5, 6, 1, 8],
[1, 5, 6, 3, 8, 2, 7, 9, 4],
[7, 2, 4, 8, 6, 3, 9, 5, 1],
[6, 8, 5, 1, 4, 9, 3, 7, 2],
[9, 3, 1, 2, 5, 7, 8, 4, 6]])
hints の各節は、ブール変数1つだけの節（[x]）です。これは「x は必ず真」という意味になります。
解を検算する#
出てきた盤面が、本当に数独の条件を満たしているか確かめます。
def check_sudoku(grid):
"""各行・各列・各ブロックが 1〜9 の1回ずつになっているか"""
target = set(range(1, 10))
rows = [set(grid[i, :]) for i in range(9)]
cols = [set(grid[:, j]) for j in range(9)]
blocks = [set(grid[i:i+3, j:j+3].ravel())
for i in range(0, 9, 3) for j in range(0, 9, 3)]
return all(g == target for g in rows + cols + blocks)
print(check_sudoku(solution))
# ヒントと矛盾していないかも確認する
print(bool((puzzle[puzzle > 0] == solution[puzzle > 0]).all()))
True
True
assumptions でヒントを与える#
hints のように「変数1つだけの節」を追加する代わりに、solve() の assumptions に真にしたい変数のリストを渡すこともできます。
assumptions は解を求めたあとにクリアされるので、同じソルバーを使い回して、別の問題を続けて解けます。何問も解くときはこちらが便利です。
sat = Solver()
sat.append_formula(conditions) # 数独のルールだけ入れておく
sat.solve(assumptions=bools[r, c, v].tolist()) # ヒントは仮定として渡す
solution2 = format_solution(sat.get_model())
print(check_sudoku(solution2))
solution2
True
array([[3, 6, 2, 7, 9, 4, 1, 8, 5],
[4, 1, 7, 5, 3, 8, 2, 6, 9],
[5, 9, 8, 6, 2, 1, 4, 3, 7],
[8, 7, 9, 4, 1, 6, 5, 2, 3],
[2, 4, 3, 9, 7, 5, 6, 1, 8],
[1, 5, 6, 3, 8, 2, 7, 9, 4],
[7, 2, 4, 8, 6, 3, 9, 5, 1],
[6, 8, 5, 1, 4, 9, 3, 7, 2],
[9, 3, 1, 2, 5, 7, 8, 4, 6]])
まとめ#
SAT は「ブール式を真にする割り当てがあるか」を解く問題。NP完全だが、実用的なソルバーは速い
ソルバーが受け取れるのは CNF（連言標準形）。DNF はそのままでは渡せない
「ちょうど1つが真」は「少なくとも1つ」＋「2つの組すべてで
~X | ~Y」で CNF になる数独は「マス × 数字」をブール変数にして、マス・行・列・ブロックの4種類の制約を並べる
swapaxes()でグループにしたい軸を最後に移すと、同じ関数を使い回せる初期盤面は
hints（変数1つの節）でもassumptionsでも与えられる。使い回すならassumptions