096. 数独(Su Doku)

数独一种非常流行的解迷游戏，它的起源并不是清楚，但是发明它的功劳部分归功于欧拉，他发明了一种类似的并且更加困难的解迷游戏，叫做拉丁方块。数独游戏的目标是把一个九乘九的网格中空白或者是零的位置替换成其它数字，使得网络中的每行、每列以及三乘三的小方块中只包含一至九九个数字。下面是一个数独的例子，左边是迷题，右边是答案：

一个良好构建的数独应该只有一个唯一解，能够用逻辑分析来解决，尽管必须用一些猜测与检验的方法来排队一些选项(对此有很多有争议的意见)。搜索的复杂性决定了问题的难度， 上面的这个示例很简单，因此可以通过直接的逻辑推理来解决。

文本文件sudoku.txt 包含50个难度各不相同的数独问题，但每一个都有唯一解(文件中的第一个问题就是上面的示例 )。请你通过解决所有五十个数独问题，找到每个解答中网格左上角的三位数字之和(例如，上面示例网格答案的左上角找到的三位数字是483)。

分析：数独是一个非常流行的解迷游戏，通常需要我们使用大量的逻辑推理来得到答案，可以对人的思维能力有很大的提升。同时解数独也是一个非常经典的编程问题，可以用多种经典的算法来解决。

一、解决数独问题的三种思路

我们最先想到的可能是把各种之前归纳总结出来的解题技巧用编程语言表达出来，再用它们来求解数独。这种方法的问题在于相应的解题规则多达数十条，全都用代码写出来会非常繁琐臃肿，并且这些方法通常都有一定的局限性，不能确保求解出所有的数独问题。

第二种可以想到的方法就是回溯法，它本质是一种深度优先搜索，是一种聪明的暴力解法。我们通过空白网络中不断尝试不同的数字，如果满足要求则继续向前尝试，如果不满足则退回到原来网络重新尝试，直到我们尝试出一个可行的答案，把数独问题解决。这种方法的优点是可以确保解答出所有数独问题，缺陷是效率较低，对于部分数独问题，可能耗时很久。因此，我们可以把第一种和第二种方法结合起来，结合逻辑推理与回溯法的优点，即确保算法能够解决所有数独问题，也使得算法有较高的效率。这种路径的代表大家可能参见这个网站，我觉得他基本上做到了这种方法的极致。

我这里重点讲解了第三种方法，也是我使用的方法，也就是把数独问题转化为一个精确覆盖问题(exact cover problem)，再通过高效的dancing links算法加以解决，最终得到原数独问题的解。相比而言，这种方法的通用性更强，除开解决数独问题以外，还可以解决其它的约束满足问题，比如著名的八皇后问题。此外，这种算法的性能也更好，相比于大部分逻辑推理加回溯法的方法，能够有五十倍以上的性能提升。下面我具体讲解以下这种算法的思路：

二、精确覆盖问题

精确覆盖问题是指假设一个全集\(X\)中若干子集的集合为\(S\)，如果\(S\)的子集\(S^*\)满足\(X\)中的每一个元素在\(S^*\)中恰好出现一次，那么集合\(S^*\)就是\(X\)的精确覆盖。简单来说，就是要找集合\(X\)的一系列子集，使得这些子集两两之间不相交，并且这些子集的并集恰好等于集合\(X\)。上面的定义可能有些难以理解，可以看一个具体的例子：假如有一个集合\(X=\{1, 2, 3, 4, 5, 6, 7\}\)，它有一系列的子集如下：

• \(A=\{1,4,7\}\)
• \(B=\{1,4\}\)
• \(C=\{4,5,7\}\)
• \(D=\{3,5,6\}\)
• \(E=\{2,3,6,7\}\)
• \(F=\{2,7\}\)

求解精确覆盖问题就是要在这些子集中找到一系列集合，使得它们之间两两不相交同时它们的并集要等于原集合\(X\)。经过简单的尝试，我们可以发现集合\(S^*=\{B,D,F\}\)满足条件，所以集合\(S^*\)就是集合\(X\)一个精确覆盖。这里的例子因为很简单，所以我们只需要简单尝试一下就能够找到答案，但是如果问题更为复杂一些，元素和子集的数是更多，我们就很难简单尝试出答案，需要一种系统的解决方法。

对于精确覆盖问题，除开上面集合形式的表示方法以外，我们还可以用矩阵的形式来表示，比如上面的例子我们可以表示成如下矩阵的形式：

这里矩阵中的行标题为原集合\(X\)中包含的元素，也可以理解为可供我们选择的元素或条件的范围，列标题为所有待选的子集。以行为单位看，可以看到每个子集包含的元素，也就是单元格为一的元素，比如集合\(A\)就包含元素1,4,7；以列为单位看，可以看到每个元素包含在那些集合中，比如数字一就包含在集合\(A\)与集合\(B\)中。通过这种转化，我们求解精确覆盖问题就变成在这个矩阵中找到一系列的行，使得这些行对应的每一列都只有一个数字一，也就是表格中标红的\(B,D,F\)三行。

把需要求解的精确覆盖问题用矩阵形式表示出来后，我们还需要一个高效的算法对这个矩阵进行处理，只保留那些满足条件的行，最终得到一个较小的矩阵。拜高德纳大神所赐，我们有了这样一个算法，高大神把它称为Algorithm X，有时候也叫做dancing links algorithm，这主要这个算法用到一个专门的数据结构dancing links，中文有时称做舞蹈链或者跳舞链。从上面举的那个简单例子我们可以看出，求解精确覆盖问题实际上就是对这个矩阵的行与列进行处理，Algorithm X对这种思路进行了细化，明确了处理的步骤，提高的了算法的效率。算法的具体执行过程我就不详细介绍了，大家可以参考这个维基百科，里面有非常详细的解释；也可以参考一个这篇文章，里面也舞蹈链算法做了一个非常直观的解释。

有两点补充说明一下：一是Algorithm X本质上仍是一种使用了回溯的深度优先搜索算法；二是dancing links这种数据结构是Algorithm X之所以较为高效的关键，因为精确覆盖问题对应的表示矩阵一般是一个规模较大的稀疏矩阵，值为零的元素远远多于值为一的元素，所以直接对矩阵的行与列进行删除处理非常的低效，如果使用舞蹈链这种数据结构就可以避免这个问题。

三、数独问题求解

之前已经提到了可以把数独问题转化为精确覆盖问题求解，这种思路能够成立需要有两个条件：一是精确覆盖问题本身能够高效的求解，否则我们转化就没有意义了，这个问题上面已经解决了，精确覆盖问题确实存在一个高效的求解算法；二是数独问题能够高效地转化为精确覆盖问题，如果这种转换不够高效，整个算法的效率也会受到拖累。在此之前，我们需要看一下如何把数独问题转化为精确覆盖问题。

我们知道数独问题是一类典型的约束满足问题，即求满足特定约束条件的解，比如对于数独来说，它的约束条件主要有四个：

• 每一行必须包含且仅能包含数字一至九；
• 每一列必须包含且仅能包含数字一至九；
• 每一个三乘三的小网格必须包含且仅能包含数字一至九；
• 每个单元格只能填写一个数字。

这四个条件中前三个是显性的，最后一个是隐含的，但也是非常关键的条件。

精确覆盖问题实际上也可以看作是一种约束满足问题，在之前的例子里，列标题可以看作是需要满足的限制或约束条件，而行标题可以看作是待选的问题的解，我们的目标就是保留一些特定的行，使他们两两之间互不排斥并且合起来可以满足全部的约束条件。因此对于数独问题，它的列标题应该就是上面的四个约束条件，我们一个一个来看，首先每一行必须包含且仅能包含数字一至九，因为总共有九行且每行有九个选择(数字的顺序不重要，所以是组合而不是排列，故为\(C^1_9\))，所以总共有八十一个选择。同理，剩下三个约束条件每个也有八十个选择，所以列标题的选择项共有\(81\times4=324\)个选择。对于列标题，它表示每个单元格可供选择的数字，因为总共有八十一个单元格，每个单元格有九个选择，所以总共有\(9^3=729\)个选择。因此九乘九数独转化为精确覆盖问题的结果是一个\(724\)行和\(324\)列的矩阵。

关于具体如何把数独转化为精确覆盖问题的矩阵，大家可以参见这个知乎回答。关于一个九乘九数独的完整转化矩阵，大家可以参见这个网址。具体到Algorithm X的实现，我直接借用了这个网址中一个非常高效的实现，它使用python内置的字典实现了类似舞蹈链的功能。对代码略做修改，我们就可以这个算法直接求解题目中给出的五十个数独问题，求解五十个数独只需要花大约两百毫秒。以下为代码：

# time cost = 221 ms ± 2.83 ms

from itertools import product

def exact_cover(X, Y):
    X = {j: set() for j in X}
    for i, row in Y.items():
        for j in row:
            X[j].add(i)
    return X, Y

def select(X, Y, r):
    cols = []
    for j in Y[r]:
        for i in X[j]:
            for k in Y[i]:
                if k != j:
                    X[k].remove(i)
        cols.append(X.pop(j))
    return cols

def deselect(X, Y, r, cols):
    for j in reversed(Y[r]):
        X[j] = cols.pop()
        for i in X[j]:
            for k in Y[i]:
                if k != j:
                    X[k].add(i)

def solve(X, Y, solution):
    if not X:
        yield list(solution)
    else:
        c = min(X, key=lambda c: len(X[c]))
        for r in list(X[c]):
            solution.append(r)
            cols = select(X, Y, r)
            for s in solve(X, Y, solution):
                yield s
            deselect(X, Y, r, cols)
            solution.pop()

def solve_sudoku(size, grid):
    R, C = size
    N = R * C
    X = ([("rc", rc) for rc in product(range(N), range(N))] +
         [("rn", rn) for rn in product(range(N), range(1, N + 1))] +
         [("cn", cn) for cn in product(range(N), range(1, N + 1))] +
         [("bn", bn) for bn in product(range(N), range(1, N + 1))])
    Y = dict()
    for r, c, n in product(range(N), range(N), range(1, N + 1)):
        b = (r // R) * R + (c // C) # Box number
        Y[(r, c, n)] = [
            ("rc", (r, c)),
            ("rn", (r, n)),
            ("cn", (c, n)),
            ("bn", (b, n))]
    X, Y = exact_cover(X, Y)
    for i, row in enumerate(grid):
        for j, n in enumerate(row):
            if n:
                select(X, Y, (i, j, n))
    for solution in solve(X, Y, []):
        for (r, c, n) in solution:
            grid[r][c] = n
        yield grid

def import_grids():
    grids = []
    with open('data/pe096.txt') as f:
        data = f.readlines()
    for i in range(1,492,10):
        grid = [[int(x) for x in row.strip()] for row in data[i:i+9]]
        grids.append(grid)
    return grids

def main():
    ans = 0
    grids = import_grids()
    for grid in grids:
        g = list(solve_sudoku((3, 3), grid))[0]
        ans += (g[0][0]*100+g[0][1]*10+g[0][2])
    return ans