数独解法小探

Posted 2020-07-29 kaige

数独的游戏要求在一个9X9的格子内填入1~9的数字，使得每一行，每一列，以及九个3X3的子区域内都没有重复的数字。如何用程序的方法来解这个问题呢？

 

稍作思索，我写出了第一种解法。从事后查询维基百科¹来看，这种方法可以称之为回溯法。思路很简单，依次扫描每一个待填数字的空格：

1. 在第一个空格里面填上“1”，检查这个数字是否合法（其所在的行、列，以及3X3的子区域里不存在重复的数字）。如果合法，则前进到第二个格子。否则，在这个格子里继续试2，3，… ，直到合法为止。

2. 在第二个格子里面继续填数字，从“1”开始试起，直到找到一个合法的数字，继续进到下一格。

3. 如果在某个格子里，从“1”到“9”都不合法，这说明前面某个格子填错了。这时就回退到上一格，从还没试的数字里面继续试。例如，如果上一格已填的数字是3，就继续试4，5，6，… 是否合法。如果找到一个合法的数字，则又前进到下一格。如果找不到，说明前面还有格子也填错了，则继续回退到更前面一格，… 如此反复。

4. 如果这个数独是有解的，我们总会遇到“每个格子都碰巧填对了”的情况。如果这个数独是没有解的，那么我们会遇到“第一个格子试了1~9所有的数字都不行”的情况。

用C++实现了一下²。求解上面的例子很快，基本上是秒杀。

 

按理到这儿就结束了。不过，看上去这个解法只能给出数独的一个解，有没有什么方法能够求出某个数独的所有解呢？

数独的游戏就是往9X9个格子里填入1~9的数字，那么一共有9X9X9 … X9 = 9⁸¹ 种排列组合的情况。如果把81个格子摊平，写成一个81位长的数字，那么从0 0 0 0 0 … 0 到 9 9 9 9 9 … 9 就穷尽了所有的情况。但是显然，我们不可能去验证这所有9⁸¹种情况。目前Intel i7处理器每秒约可执行2.38 X 10¹¹条指令³，即使我们的程序只有9⁸¹条指令，它也无法在我们的有生之年算出结果。

有没有办法减少待验证的解的数目呢？第一个想法是我们只需要排列组合空格子里的数字，那么对于本文的例子而言，这将是9⁵¹种情况——虽然减少了很多，但对我们的CPU而言这仍是一个不切实际的数字。第二个想法是，对于这51个空格，每个格子也没必要穷尽1~9里面所有的数字。考虑每个格子在同行、同列、同3X3子区域里非空格子里的数字，它可能的取值个数实际上是小于9的，也许有的格子的值就是唯一确定的!于是很快的实现了这个想法。但是，很遗憾的发现，对于本文的例子，我们只是把待验证解的数量降到了约3X10²³。虽然相比于9⁵¹这又是一次飞跃，不过这仍然无法导致一个可用的程序。

难道这种排列组合的方法就没有前途吗？

经过思考，我找到了答案：可以先对每一个3X3的子区域找到所有可能的解。整张数独表的解则由9个3X3子区域的解排列组合而成。求单个3X3子区域的解应该是很快的，特别是在我们已经根据全局数独表限制了每个格子可能取值的范围的情况下。求出9个子区域的所有可能解，再对这9个子区域的解作排列组合也应该是很快的。带着这样的想法，我实现了第二版：所有待验证解的数目降到了2X10⁸，程序运行后，花了10秒钟左右就吐出了答案。

还能更快点吗？当然！方法是中间再加一层：9个3X3的子区域排成了三行，可以先对每一行找到所有合法的解，然后再对三行的解进行排列组合得到整张表的解。很快实现了这一想法：对于本文的例子，和回溯法一样，也是秒杀。最后总的待验证解的数目降到了20000左右。

这种解法可称作排列组合法，思路总结如下：

1. 对于每一空格子，考虑其所在行、列以及子区域，找出所有可能取值的列表S1

2. 对每一个3X3的子区域，对其包含的所有空格子的取值S1进行排列组合，找出该子区域的所有解的列表S2.

3. 9个3X3的子区域排成了三行。对每一行，对其包含的三个子区域的解S2排列组合，找出这一行的解的列表S3.

4. 对三行的解S3排列组合，找出整张数独表的解。

至此，排列组合法的研究就告一段落了。虽然在实现上还可以有一些优化，例如加上并行的处理，不过在方法的思路上也就大致如此了。

 

还有没有其他的解数独的思路呢？搜索了一下维基百科，发现回溯法已有提及，排列组合法倒似未见描述。除此之外，还有三种将数独问题抽象成不同数学问题的方法值得一提。

第一种是把数独的问题抽象为一个精确覆盖问题，再用解精确覆盖问题的算法如舞蹈链算法去解它⁴。

具体而言，我们是将解数独的问题转化为求一个“Exact Hitting Set”的问题。什么叫“Exact Hitting Set”呢？例如，我们有如下集合：

A = {1, 4, 7};

B = {1, 4};

C = {4, 5, 7};

D = {3, 5, 6};

E = {2, 3, 6, 7};

F = {2, 7};

其中每个集合都是集合X = {1, 2, 3, 4, 5, 6, 7}的一个子集。于是，集合X^*={1, 2, 5}被称为是{A, B, C, D, E, F}的一个Exact Hitting Set。这是因为A, B, C, D, E, F都恰好只包含X^*中的一个元素。

怎么把数独的问题抽象为求一个Exact Hitting Set的问题呢？首先，数独的问题可以表述为在9X9格子内填入1~9的数字，那么所有的格子加起来就有9X9X9=729种取值。不妨记为:

X= {R1C1#1, R1C1#2,…, R9C9#9}

数独的游戏规则可以描述为如下四条规则：

1. 每个格子只能填1个1~9之间的数。例如对于格子R1C1其所有取值的可能性的集合为:

R1C1 = {R1C1#1, R1C1#2, R1C1#3, R1C1#4, R1C1#5, R1C1#6, R1C1#7, R1C1#8, R1C1#9}

2. 每一行中每个数字必须出现且仅出现一次。例如对于第一行，数字1出现的位置有这些可能: R1#1 = {R1C1#1, R1C2#1, R1C3#1, R1C4#1, R1C5#1, R1C6#1, R1C7#1, R1C8#1, R1C9#1}

3. 每一列中每个数字必须出现且仅出现一次。同上，对于第一列，数字1可能出现的位置为：C1#1 = {R1C1#1, R2C1#1, R3C1#1, R4C1#1, R5C1#1, R6C1#1, R7C1#1, R8C1#1, R9C1#1}

4. 每一个3X3子区域里，每个数字必须出现且仅出现一次。同上对于第一个3X3子区域，数字1可能出现的位置为：

B1#1 = {R1C1#1, R1C2#1, R1C3#1, R2C1#1, R2C2#1, R2C3#1, R3C1#1, R3C2#1, R3C3#1}

每条规则中每个数字可能出现的位置都是一个集合。我们有9行，9列，81个格子，9个待填的数字，那么我们有81+81+81+81=324个集合。这些集合都是X的子集。则解数独相当于找到X的一个子集X^*，使X^*中的每一个元素在这324个子集中出现且仅出现一次——这正是一个Exact Hitting Set!

 

第二种是将数独的问题表达为一个优化的问题，再用求解优化问题的算法去解。问题的核心是定义一个评价函数：将数独表里待填的数字当作自变量，将当前整个表格与有效解局面的区别程度当作函数值，数独问题即转化为一个优化问题：当数独表里填哪些数字时，评价函数的值最小？下面以经典的优化算法模拟退火法为例，简单介绍一下这类方法的工作原理⁷。

首先我们给数独表里的待填数字，也就是优化问题中的自变量，定义一个初始值。方法是在每个3X3的子区域的空格内，随机填入1~9中的数字，使得1~9在这个子区域内没有重复。其次我们给数独表定义一个操作，让它可以从自变量的当前值“搜索”得到一个新的值。我们将操作定义为：在9X9数独表内随机选取一个填好数字的空格，将它的值与其所在的3X3子区域内另一随机选取的空格的值进行交换。可见，在这样定义的初始值和操作方法下，我们得到的数独表永远满足“在3X3子区域内没有重复”这一条件。于是，数独表的评价函数就可以简化定义为：当前所有行、列在1~9中缺少的数字的个数之和。例如，在下左图中，第一列只缺少数字9，所以缺少的数字个数是1。第二列缺少数字6，8，所以缺少的数字个数是2，…，整张表所有行、列加起来缺少的数字个数是34。而在下右的有效解图中，所有行列缺少的数字个数之和是0。

 定义好了初始值、操作方法以及评价函数，我们就可以用标准的模拟退火算法来解这个问题了。将数独表的初始值记为X₀，操作方法记为Op(X)，评价函数记为F(X)，则算法流程如下⁶：

1. 选定一个初始的“温度”值t，及其最小可取值t_min

2. 运用定义的操作方法Op，得到一个新的值X₁。X₁ = Op(X₀)

3. 计算评价函数F(X₀)和F(X₁)，计算差值dF = F(X₁) - F(X₀)

4. 如果dF < 0，这意味着新的值更接近一个有效解。则将数独表的状态更新到X₁，也即令X₀=X₁。

5. 如果dF > 0，我们仍然以一定的概率将数独表的状态更新到X₁，这是为了使我们的搜索能跳出评价函数的局部最小值。这个概率条件通常定义为：e^(-dF/t) > random(0, 1)。即当自然指数e^(-dF/t) 大于0和1之间的一个随机数值时，X₀=X₁。因为-dF/t总是小于0的，所以e^(-dF/t)总是在0和1之间。这意味着温度值t越大时，将状态更新到X₁的概率也越大。随着温度值减小，这种概率也减小。

6. 以某种策略降低温度值，如t = R * t，R是0到1之间的一个常数。如果 t < t_min，退出，否则重回步骤2。

这种随机搜索的方法将带着我们从一个随机给定的初始值，最终搜索到数独表的一个有效解！感觉有点神奇吧？有论文显示，这种算法对9X9的数独有100%的成功率，对16X16，25X25的数独也有相当高的成功率⁷。

 

第三种方法把数独看成是一个约束求解问题，然后用约束编程的方法去解。这其实是一个很自然的理解方式：把数独表里的每一个空格看成一个变量，这些变量的可取值范围是1~9之间的整数。“每行、每列及每个3X3子区域内的变量均各不相等”就是这些变量要满足的约束。所以解数独的问题也就是一个约束求解问题：当这些变量取什么值时，所有的约束都能同时被满足？

利用一些约束编程的语言如Prolog，我们很容易写出数独的求解程序。和普通的编程语言不同，利用约束编程语言，我们只要将这个约束问题描述出来，就可以自动得到问题的解，而不需指定求解的具体步骤。下面给出一段利用Prolog语言求解数独的程序。为简化起见，这个数独是4X4的⁸。

待求解的数独见下左图。首先将数独表的每一个空格都标记为一个变量，如下右图所示。

则求解的程序如下：

sudoku(Puzzle, Solution) :-
  Solution = Puzzle,

  Puzzle = [A, B, C, D,
            E, F, G, H,
            I, J, K, L,
            M, N, O, P],
  
  fd_domain(Solution, 1, 4),
  
  Row1 = [A, B, C, D],
  Row2 = [E, F, G, H],
  Row3 = [I, J, K, L],
  Row4 = [M, N, O, P],      
  
  Col1 = [A, E, I, M],
  Col2 = [B, F, J, N],
  Col3 = [C, G, K, O],
  Col4 = [D, H, L, P],      
  
  Square1 = [A, B, E, F],
  Square2 = [C, D, G, H],
  Square3 = [I, J, M, N],
  Square4 = [K, L, O, P],      
  
  valid([Row1, Row2, Row3, Row4,
         Col1, Col2, Col3, Col4,
         Square1, Square2, Square3, Square4]).
 
valid([]).
valid([Head | Tail]) :- fd_all_different(Head), valid(Tail).

这段程序首先将4X4数独表定义成16个变量组成的一个集合。然后定义每一行、每一列、每一2X2子区域内变量组成的子集。最后定义各个子集要满足的约束：其包含的变量之间互不相等(all_different)。

运行时，只要按格式输入待解的数独表：

| ?- sudoku([_, _, 4, _,
             _, 2, _, _,
             _, _, 1, _,
             _, _, 3, _],
             Solution).

就可以得到结果:

S = [3,1,4,2,4,2,3,1,2,4,1,3,1,3,2,4]

看上去似乎很简单。不过，约束编程到底是如何求数独的解的呢?

一种方法是将约束系统表达为图的形式⁹。

将变量当成图的顶点，将可取值的范围写在顶点旁边。

 

将变量之间的约束当成图的边。在这里我们只有一种约束：两个变量不相等(all_different)。如果某两个变量之间有“不相等”的约束，我们就在代表它们的顶点之间加一条边。例如，加上“每行的变量互不相等”的约束，我们的图就变成这样：

 

加上“每列的变量互不相等”，我们的图会变成这样：

 

将每行、每列以及每个2X2小区域要满足的约束叠加起来，最终的图是这样：

 

下一步就是更新各个变量的可取值范围。当某个变量的取值范围有变时，我们就在这个变量的顶点旁边标一个星号（*）。首先更新代表非空格子的变量的取值范围。我们有四个非空格子C，F，K，N，所以把它们的取值范围设为格子里的数字，然后给顶点标上星号。

 

接下来更新标有星号顶点的相邻顶点的变量取值范围。以顶点C为例，因为它只能取值4了，所以它所有的相邻顶点就不能再取值4了。于是将A，B，D，O，H，G的取值范围变成 {1, 2, 3} ，因为它们的变量取值范围有变，所以把它们的顶点标上星号，…，以此类推。

 

就这样，继续从标有星号的顶点出发，更新它们相邻顶点的取值范围。如果任何顶点的取值范围变了，就又把它们标为星号。注意到每一次更新顶点的取值范围时，我们都只考虑当前的一对顶点，而不需考虑其它的顶点。所以这其实是一个简单的重复过程。

如果当前的数独表只有唯一解，那么我们的更新过程会收敛到如下结果，也就是数独的解！

 

至此，利用约束编程解数独的原理就介绍完了。当然这只是一个很粗浅的介绍，实际应用中还有很多地方要做特殊处理和优化，在这里就不赘述了。

我们这篇介绍数独解法的小文章到这儿也该结束了。我们一共介绍了五种方法来求解数独：回溯法，排列组合法，精确覆盖问题法，模拟退火法以及约束编程法。小小一个数独，竟可以从这么多的角度来看待和分析，不由得让人感叹思维之奇，数字之妙啊！

参考资料

1. https://en.wikipedia.org/wiki/Sudoku_solving_algorithms

2. https://github.com/kaige/Sudoku/

3. https://en.wikipedia.org/wiki/Instructions_per_second#Millions_of_instructions_per_second

4. https://en.wikipedia.org/wiki/Exact_cover#Sudoku

5. Rhyd Lewis. Metaheuristics can Solve Sudoku Puzzles.

6. http://www.cnblogs.com/heaad/archive/2010/12/20/1911614.html

7. Perez, Meir and Marwala, Tshilidzi. Stochastic Optimization Approaches for Solving Sudoku

8. http://www.ybrikman.com/writing/2012/02/16/seven-languages-in-seven-weeks-prolog_16/

9. https://www.cl.cam.ac.uk/teaching/0809/Prolog/Prolog08ML6R2.pdf