用 SAT 解数独：把恰好一个写成 CNF#
本章的内容#
本章用 SAT 求解器解数独，从一个小推理谜题开始，逐步过渡到数独，看问题如何被翻译成布尔式。
SAT 与 CNF 的含义，后者即合取范式
把析取范式变成 CNF 的两种做法，用 SymPy 与手工推导
一组之中只有一个为真这条约束的写法
用
pysat解数独，用assumptions给出初始盘面
本章的重点是最后一件事：把“恰好一个”写成 CNF，其余部分都是它的重复应用。
可满足性问题#
可满足性问题，即 SATisfiability problem，简称 SAT，问的是能否给出一组变量的值，使给定的布尔式为真。
SAT 是最早被证明为 NP 完全的问题，与计算科学的理论基础、算法、人工智能、硬件设计等许多领域都相关，而实用的 SAT 求解器速度极快，因此把待解问题转换成 SAT 这一手法被广泛使用。NP 完全意味着没有已知的多项式时间算法，实用求解器靠冲突学习与启发式把常见情形做得很快。
本章用 pysat 解数独，先用一个小的推理谜题，确认布尔式与 CNF 的写法，从这里可以看出 SAT 表述问题的方式。推理谜题只有四个变量，数独有 729 个，两者写法完全一样，只是规模不同。
推理谜题
四名嫌疑人作了如下的供述。
甲：我不是犯人
乙：丁是犯人
丙：乙是犯人
丁：乙想陷害我
四人之中只有一个人说了真话，犯人是哪一位？
布尔变量 A、B、C、D 分别表示甲、乙、丙、丁是犯人，而 ~A 表示甲不是犯人，只要把各人的供述分别写成布尔式即可。
嫌疑人 |
供述 |
布尔式 |
|---|---|---|
甲 |
我不是犯人 |
S1 = ~A |
乙 |
丁是犯人 |
S2 = D |
丙 |
乙是犯人 |
S3 = B |
丁 |
乙想陷害我 |
S4 = ~D |
说真话的只有一个人，因此 S1〜S4 之中恰好一个为真。这一步把口语里的“只有一个人说真话”翻译成了计数条件，后面的变形都围绕它展开。
恰好一个为真，先直白地写出来，其中 & 是合取，| 是析取。
S1 & ~S2 & ~S3 & ~S4 |
~S1 & S2 & ~S3 & ~S4 |
~S1 & ~S2 & S3 & ~S4 |
~S1 & ~S2 & ~S3 & S4
这是用 OR 连起来的若干 AND，称为析取范式，即 DNF。SAT 求解器只接受合取范式，即 CNF，也就是把 OR 连成的子句再用 AND 并起来，两者的差别只在括号的层次，因此要把恰好一个为真换一种写法。
不能有两个以上为真 →
~(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 这样恒为真的子句也不会被消去。数独这类较大的问题，这种膨胀就直接体现为求解时间。因此实务上先让 SymPy 化简，再交给求解器。
两种 CNF 应当得到相同的答案。
from sympy.logic.inference import satisfiable
print(satisfiable(to_cnf(dnf, simplify=True)))
print(satisfiable(cnf_raw))
{A: True, B: False}
{A: True, B: False, D: False}
satisfiable() 返回使式子为真的变量赋值。这里 A 为真，因此甲是犯人，而 C 没有出现在式子中，取哪个值都可以，即没有约束。
用 pysat 求解#
CNF 可以用嵌套的列表表示，一个整数对应一个布尔变量，正数表示该变量为真，负数表示取反，因此 1 就是 A，-1 就是 ~A，读起来很直接。
把上面的 CNF 原样写成数值列表即可，此时 A=1, B=2, C=3, D=4。
用
append_formula()把 CNF 加入求解器，可以反复调用用
solve()求解用
get_model()取出解，正数表示真
编号与布尔变量的对应关系要自己记住，SAT 求解器只认整数。
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〜9 各出现一次
每个 3×3 区块中 1〜9 各出现一次
黑色的数字是题面，灰色的数字是答案。
SAT 的变量只有 True 与 False 两个值，因此把格 (r, c) 中填入数字 n 做成一个布尔变量，用三维数组 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)
数独的约束可以整理成四条。
# |
约束 |
分组的做法 |
|---|---|---|
1 |
每格只有一个数字 |
把数字所在的轴，即第 2 轴作为一组 |
2 |
每行中同一数字只有一个 |
按行与数字分组 |
3 |
每列中同一数字只有一个 |
按列与数字分组 |
4 |
每个区块中同一数字只有一个 |
按区块与数字分组 |
这四条都是同一个形状，即某一组之中恰好一个为真。把四种分组化成同一种形状之后，写约束的代码只需要一个函数。
用 CNF 表达恰好一个#
变量 A、B、C 中只有一个为真，这条约束可以分成两个条件。
至少有一个为真：\(A | B | C\)
不能有两个以上为真：\(\neg(A \& B) \land \neg(A \& C) \land \neg(B \& C)\)
用德摩根律，第二个条件变成 \( (\neg A | \neg B) \land (\neg A | \neg C) \land (\neg B | \neg C) \)。
若元素有 N 个，则有一个全用 OR 连起来的子句，与两两组合对应的 \(N(N-1)/2\) 个子句，即 ~X | ~Y 这样的子句，两两组合用 itertools.combinations 生成。N 一大，子句数目按平方增长，数独里子句最多的部分正在这里。
get_conditions(bools) 以 bools 的最后一轴为组，返回这条约束，做法分三步。
对最后一轴的元素数
n，造出全部两两组合把每一组用 OR 连成一个子句，即至少一个为真
对每个组合造出
~X | ~Y的子句，即两个不能同时为真
from itertools import combinations
def get_conditions(bools):
"""把最后一轴作为一组，生成“恰好一个为 True”的 CNF"""
n = bools.shape[-1]
pairs = np.array(list(combinations(range(n), 2)))
conditions = []
# 1. 至少一个为 True
conditions.extend(bools.reshape(-1, n).tolist())
# 2. 两个不能同时为 True
conditions.extend((-bools[..., pairs].reshape(-1, 2)).tolist())
return conditions
下面用小例子确认，[[1, 2, 3], [4, 5, 6]] 表示有两个组，每组各含三个变量。两组的子句互不相交，可以分别检查。
[1, 2, 3]、[4, 5, 6]：每组之中至少有一个为真[-1, -2]、[-1, -3]、[-2, -3]：1〜3 中不能有两个同时为真[-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，每格一个数字：直接传入，最后一轴就是数字
约束 2，每行一个数字：把列与数字对调
约束 3，每列一个数字：把行与数字对调
c1 = get_conditions(bools) # 每个格子只有一个数字
c2 = get_conditions(np.swapaxes(bools, 1, 2)) # 每行中同一数字只有一个
c3 = get_conditions(np.swapaxes(bools, 0, 2)) # 每列中同一数字只有一个
约束 4，即区块的约束，略为复杂一些。bools 的轴是 (行, 列, 数字)，因此要把行与列分解成区块编号与区块内的位置。
reshape(3, 3, 3, 3, 9)，即(区块行, 行内, 区块列, 列内, 数字)swapaxes(1, 2)，即(区块行, 区块列, 行内, 列内, 数字)reshape(9, 9, 9)，即(区块编号, 区块内的 9 格, 数字)swapaxes(1, 2)，即(区块编号, 数字, 区块内的 9 格)，组落在最后一轴
这样就得到每个区块、每个数字中，9 格之内只有一个为真的约束。三次 reshape 与两次 swapaxes 之后，组的轴落在最后，与行、列两种情形对齐。
# 按 (区块, 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)) # 每个区块中同一数字只有一个
conditions = c1 + c2 + c3 + c4
print(len(conditions), "个子句")
11988 个子句
解一个空盘面#
先在没有提示的状态下求解，满足数独约束的答案有无数个，求解器只给出其中的一个。空盘面能解，说明约束没有写错，这一步相当于先验证模型再填题目。
模型是变量编号的列表，因此要还原成 (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 的每个子句都只含一个布尔变量，即 [x]，表示 x 必须为真，写进求解器就等于把这一格固定住。这样的子句只有一个文字，求解器会在预处理阶段直接定下这个变量。
验算解#
确认得到的盘面是否真的满足数独的条件。
def check_sudoku(grid):
"""每行、每列、每个区块是否都是 1〜9 各一次"""
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 这样只含一个变量的子句，也可以把希望为真的变量列表交给 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，即合取范式，析取范式不能直接传入
恰好一个为真由至少一个为真与所有两两组合的
~X | ~Y合起来构成 CNF数独把格与数字的组合做成布尔变量，再排上格、行、列、区块四类约束
用
swapaxes()把要作为组的轴移到最后，同一个函数就能反复使用初始盘面可以用
hints给出，也可以用assumptions给出，反复求解时用assumptions更方便