精确覆盖问题和 DLX 算法

前言
¶

很早就想要补一下舞蹈链和精确覆盖算法，却一直各种拖延，趁着最近有空，又重新翻开了刘汝佳的书。大学的时候看过几次，甚至照着书里的思路手敲了一遍并通过了例题，但对于算法的原理一直有些不求甚解。很早以前，一位同学告诉我说“现在看不懂的东西不用勉强，以后慢慢就会懂了”，后来也真的在不断印证这句话；但我很担心随着年纪的增长，记忆力和学习能力不断退化之后，恐怕这个 flag 会逐渐倒下。所以趁着眼下尚能理解进去，尽量用自己的语言做一下记录。

精确覆盖问题
¶

有一些由整数 $1 \sim n$ 中的数字组成的集合 $S_{1}, S_{2}, \dots, S_{m}$ ，要求选择若干个集合 $S_{i}$ ，使得 $1 \sim n$ 中每个整数都在选出的集合中的某个出现且恰好仅出现一次。举个栗子：

不妨假设 $n = 7, m = 6$ ，集合为：
$\begin{aligned} S_{1} & = {1, 4, 7} \\ S_{2} & = {1, 4} \\ S_{3} & = {4, 5, 7} \\ S_{4} & = {3, 5, 6} \\ S_{5} & = {2, 3, 6, 7} \\ S_{6} & = {2, 7} \end{aligned}$
则一个精确覆盖为 ${S_{2}, S_{4}, S_{6}}$ ，因为 ${1, 4}$ , ${3, 5, 6}$ , ${2, 7}$ 无重复、无遗漏地包含了 $1 \sim 7$ 中的所有整数。

我们可以用一个 $m \times n$ 的 $01$ 矩阵来表示集合，其中， $0$ 表示不包含， $1$ 表示包含。比如第 $(i, j)$ 个位置若为 $0$ ，则说明 $S_{i}$ 中不包含 $j$ 。上文中的栗子用矩阵表示如下所示：

$$	$1$	$2$	$3$	$4$	$5$	$6$	$7$
$S_{1}$	$1$	$0$	$0$	$1$	$0$	$0$	$1$
$S_{2}$	$1$	$0$	$0$	$1$	$0$	$0$	$0$
$S_{3}$	$0$	$0$	$0$	$1$	$1$	$0$	$1$
$S_{4}$	$0$	$0$	$1$	$0$	$1$	$1$	$0$
$S_{5}$	$0$	$1$	$1$	$0$	$0$	$1$	$1$
$S_{6}$	$0$	$1$	$0$	$0$	$0$	$0$	$1$

则精确覆盖问题可重新表述为：在一个 $m \times n$ 的 $01$ 矩阵中，选择若干行，对这些行做向量加法，得到的结果为 $(1, 1, \dots, 1)$ 。即：

选出的行里，不存在某列同时在两行中值均为 $1$ ；
所有选出的行叠加在一起的结果覆盖所有列（每一列的值都不为 $0$ ）。

算法 X
¶

算法 X（Algorithm X），其实就是回溯，可能是专门针对覆盖问题提出的算法吧。算法描述如下：

每次选择一个没有被删除列，然后枚举该列为 $1$ 的所有行，尝试删除这些行，递归搜索后再恢复这些行。尝试删除行时，还要将该行中所有值为 $1$ 的列也一并删除^[1]，恢复时也同样。
若没有可以选择的列了（即所有列都被删除了），说明已经找到一个精确覆盖的解了^[2]；若还有列但是没有行了^[3]说明原问题无解。

舞蹈链
¶

舞蹈链（Dancing Link），一种支持快速删除、恢复列和行的数据结构。

舞蹈链是一个十字型双向链表结构，链表中的每个节点对应上述 $01$ 矩阵中的一个 $1$ 。另外还有 $n + 1$ 个虚拟节点，其中每列最上方一个虚拟节点作为该列链表的头指针，而所有虚拟节点的最前方有一个虚拟节点，作为虚拟节点的头指针，它也是舞蹈链的头指针。如下图所示（图片来源于网络）：

使用四个数组 $L$ , $R$ , $U$ , $D$ 分别表示舞蹈链中节点的左、右、上、下四个方向的指针，下标为节点的编号；同时列编号作为该列的虚拟节点， $0$ 表示舞蹈链的头指针对应的虚拟节点。

删除行：只需要修改将该行中所有节点的上下行的下、上指针互指就好了。

dancing-link.01.ts

1
2
3
4
5
// 设要删除的某行中某个节点编号为 R[i]
for (let j = R[i]; j !== i; j = R[j]) {
  U[D[j]] = U[j]
  D[U[j]] = D[j]
}

恢复行：只需要将该行中所有节点的上下行的下、上指针分别指向该该节点就好了。

dancing-link.01.ts

1
2
3
4
5
// 设要删除的某行中某个节点编号为 R[i]
for (let j = R[i]; j !== i; j = R[j]) {
  U[D[j]] = j
  D[U[j]] = j
}

删除列和恢复列也类似，此处略去。

DLX 算法
¶

使用了舞蹈链的算法 X 通常被称为 DLX 算法。此处给出 DLX 算法的 Typescript 实现^[4]。

dlx.ts | 240 lines.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
/**
 * The algorithm X that applied the dancing-link, it is also called as "DLX".
 * It is used to solve the exact-cover problem.
 *
 * Dancing-link: A cross doubly linked list, each column has a virtual node as
 * the head pointer, and at the top of all virtual nodes there is an additional
 * virtual node as the head pointer of the virtual node, which is also the head
 * pointer of the entire dancing-link. In the implementation of using an array
 * to simulate a linked list, the virtual node is represented by a column
 * number, and the head pointer of the dancing-link can be represented by 0.
 *
 * @see https://me.guanghechen.com/post/algorithm/dlx/
 */
export interface DLX {
  /**
   * Initialize the dancing-link.
   * @param totalColumns   number of columns
   */
  init(totalColumns: number): void

  /**
   * Release memory variables.
   */
  destroy(): void

  /**
   * Add a row to the dancing-link.
   *
   * It should be noted that after solving the exact-cover problem, the
   * result is a list of selected row numbers, so the row number should be
   * specified as a value that can carry information.
   *
   * @param r         the row number
   * @param columns   columns on the row
   */
  addRow(rowNo: number, columns: ReadonlyArray<number>): void

  /**
   * Try to find a precise coverage.
   *
   * When a solution is found, return the row numbers of all selected rows,
   * otherwise return null.
   */
  solve(): number[] | null
}

/**
 * Generate an object that encapsulates the DLX algorithm.
 *
 * @param MAX_N   maximum number of nodes in the dancing-link
 * @returns
 */
export function createDLX(MAX_N: number): DLX {
  // The number of nodes in the dancing-link (including the virtual nodes on
  // the column).
  let sz: number

  // the number of columns in the dancing-link.
  let totalColumns: number

  const selectedRowNos: number[] = new Array(MAX_N) // list of row numbers of selected rows
  let countOfSelectedRows: number // the number of selected rows

  const count: number[] = new Array(MAX_N) // lhe number of nodes of a column in the dancing-link
  const row: number[] = new Array(MAX_N) // the row number of a node in the dancing-link
  const col: number[] = new Array(MAX_N) // the column number of a node in the dancing-link
  const L: number[] = new Array(MAX_N) // left pointer of cross-link list
  const R: number[] = new Array(MAX_N) // right pointer of cross-link list
  const U: number[] = new Array(MAX_N) // up pointer of cross-link list
  const D: number[] = new Array(MAX_N) // down pointer of cross-link list
  return { init, destroy, addRow, solve }

  /**
   * @see DLX#init
   * @public
   */
  function init(_totalColumns: number): void {
    totalColumns = _totalColumns
    sz = _totalColumns + 1

    // Resize arrays.
    if (selectedRowNos.length < sz) {
      selectedRowNos.length = sz
      count.length = sz
      row.length = sz
      col.length = sz
      L.length = sz
      R.length = sz
      U.length = sz
      D.length = sz
    }

    for (let i = 0; i < sz; ++i) {
      L[i] = i - 1
      R[i] = i + 1
      U[i] = i
      D[i] = i
    }
    R[_totalColumns] = 0
    L[0] = _totalColumns

    count.fill(0, 0, sz)
  }

  /**
   * @see DLX#destroy
   * @public
   */
  function destroy(): void {
    selectedRowNos.length = 0
    count.length = 0
    row.length = 0
    col.length = 0
    L.length = 0
    R.length = 0
    U.length = 0
    D.length = 0
  }

  /**
   * @see DLX#addRow
   * @public
   */
  function addRow(r: number, columns: ReadonlyArray<number>): void {
    const first = sz
    for (let i = 0; i < columns.length; ++i, ++sz) {
      const c = columns[i]
      row[sz] = r
      col[sz] = c
      count[c] += 1

      // Connect left and right nodes
      L[sz] = sz - 1
      R[sz] = sz + 1

      // Connect top and bottom nodes,
      // c is the virtual node on the c-th column, and is also the head pointer
      // of the linked list of the column, so at this time U[c] is the last
      // element of the column
      D[sz] = c
      D[U[c]] = sz
      U[sz] = U[c]
      U[c] = sz
    }

    // Since this is a circular linked list, the first and last columns of the
    // current row are connected to each other.
    R[sz - 1] = first
    L[first] = sz - 1
  }

  /**
   * @see DLX#solve
   * @public
   */
  function solve(): number[] | null {
    if (!algorithmX(0)) return null
    return selectedRowNos.slice(0, countOfSelectedRows)
  }

  /**
   * Remove a column from the dancing-link.
   * @param c   column number
   * @private
   */
  function removeColumn(c: number): void {
    L[R[c]] = L[c]
    R[L[c]] = R[c]
    for (let i = D[c]; i !== c; i = D[i]) {
      for (let j = R[i]; j !== i; j = R[j]) {
        U[D[j]] = U[j]
        D[U[j]] = D[j]
        count[col[j]] -= 1
      }
    }
  }

  /**
   * Restore a previously deleted column
   * @param c   column number
   * @private
   */
  function restoreColumn(c: number): void {
    for (let i = U[c]; i !== c; i = U[i]) {
      for (let j = L[i]; j !== i; j = L[j]) {
        count[col[j]] += 1
        U[D[j]] = j
        D[U[j]] = j
      }
    }
    L[R[c]] = c
    R[L[c]] = c
  }

  /**
   * Algorithm X.
   *
   * Recursively solve the problem of precise coverage, enumerate which rows are
   * selected in the recursive process, remove the selected rows and all the
   * columns on the rows, and restore these rows and columns during the
   * backtrack.
   *
   * @param dep   recursion depth
   * @private
   */
  function algorithmX(dep: number): boolean {
    // Find a solution when the dancing-link is empty.
    if (R[0] === 0) {
      // Record the length of the solution.
      countOfSelectedRows = dep
      return true
    }

    /**
     * Optimization: Find the column with the least number of nodes, and try to
     * cover from this column.
     */
    let c = R[0]
    for (let i = R[0]; i !== 0; i = R[i]) {
      if (count[i] < count[c]) c = i
    }

    // Remove this column.
    removeColumn(c)
    for (let i = D[c]; i !== c; i = D[i]) {
      selectedRowNos[dep] = row[i]
      for (let j = R[i]; j !== i; j = R[j]) removeColumn(col[j])

      // Recursively processing.
      if (algorithmX(dep + 1)) return true

      // Backtrack.
      for (let j = L[i]; j !== i; j = L[j]) restoreColumn(col[j])
    }
    // Backtrack.
    restoreColumn(c)

    return false
  }
}

我把它封装在了 @algorithm.ts/dlx 中，你可以通过 npm 包导入它：

1
2
3
4
5
6
7
8
9
10
11
12
13
import { createDLX } from 'algorithm.ts/dlx'

// 创建一个基于至多有 1000 个节点的舞蹈链的 dlx 算法
const dlx = createDLX(1000)

// 初始化 dlx 算法，总共有 1000 列
dlx.init(1000)

// 添加行，此处为伪代码
dlx.addRow(...)

// 尝试找到一个精确覆盖
dlx.solve()

求解数独问题
¶

数独是精确覆盖问题的一个特例。为了套用 DLX 算法框架，首先需要弄清楚如何构建 $01$ 矩阵。一般来说，可以将列对应成约束，而将行对应到策略，即选择某个策略时能够满足哪些约束。比如考虑经典 $x^{2} \times x^{2}$ 数独游戏，其有如下类型的约束：

$S l o t (a, b)$ : 第 $a$ 行第 $b$ 列格子要有数字；
$R o w (a, b)$ : 第 $a$ 行要有数字 $b$ ；
$C o l (a, b)$ : 第 $a$ 列要有数字 $b$ ；
$S u b (a, b)$ : 第 $a$ 个子方阵中要有方阵 $c$ ；

一共有 $x^{2} \times x^{2} \times 4$ 中约束。再考虑可选择的策略，即在第 $r$ 行第 $c$ 列填入数字 $v$ ，且可以满足约束 $S l o t (r, c)$ 、 $R o w (r, v)$ 、 $C o l (c, v)$ 、 $S u b (⌊ \frac{r}{x} ⌋ \times x + ⌊ \frac{c}{x} ⌋, v)$ 。

代码如下：

sudoku.ts | 83 lines.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
import { createDLX } from './dlx'

/**
 * Sudoku constraints.
 */
export enum SudokuConstraint {
  SLOT = 0, // Slot(a, b) 表示第 a 行 b 列个格子上要有数字
  ROW = 1, // Row(a, b) 表示第 a 行要有数字 b
  COL = 2, // Col(a, b) 表示第 a 列要有数字 b
  SUB = 3, // Sub(a, b) 表示第 a 个子方阵要有数字 b
}

export interface SudokuSolver {
  /**
   * 数独的谜题格子，从 0 开始填，若某个格子未被填，则将其置为 -1
   * @param puzzle
   */
  solve(puzzle: number[][]): boolean
}

/**
 * @param SUDOKU_SIZE_SQRT 数独子方阵大小（即数独行数的平方根)
 */
export function createSudokuSolver(SUDOKU_SIZE_SQRT: number): SudokuSolver {
  const ebs = 1e-6
  const SUDOKU_SIZE = SUDOKU_SIZE_SQRT * SUDOKU_SIZE_SQRT
  const SUDOKU_SIZE_SQUARE = SUDOKU_SIZE * SUDOKU_SIZE
  const MAX_NODES = SUDOKU_SIZE_SQUARE * 4

  let codeA: number, codeB: number, codeC: number
  const columns: number[] = new Array<number>(4)
  const solver = createDLX(MAX_NODES)
  return { solve }

  function solve(puzzle: number[][]): boolean {
    solver.init(MAX_NODES)
    for (let r = 0; r < SUDOKU_SIZE; ++r) {
      for (let c = 0; c < SUDOKU_SIZE; ++c) {
        const w = puzzle[r][c]

        // (r,c) 所属的子方阵编号
        const s =
          Math.floor(r / SUDOKU_SIZE_SQRT + ebs) * SUDOKU_SIZE_SQRT +
          Math.floor(c / SUDOKU_SIZE_SQRT + ebs)
        for (let v = 0; v < SUDOKU_SIZE; ++v) {
          if (w === -1 || w === v) {
            columns[0] = encode(SudokuConstraint.SLOT, r, c)
            columns[1] = encode(SudokuConstraint.ROW, r, v)
            columns[2] = encode(SudokuConstraint.COL, c, v)
            columns[3] = encode(SudokuConstraint.SUB, s, v)
            solver.addRow(encode(r, c, v), columns)
          }
        }
      }
    }

    const answer: number[] | null = solver.solve()
    if (answer === null) return false

    for (const code of answer) {
      decode(code)
      // eslint-disable-next-line no-param-reassign
      puzzle[codeA][codeB] = codeC
    }

    return true
  }

  function encode(a: number, b: number, c: number): number {
    return a * SUDOKU_SIZE_SQUARE + b * SUDOKU_SIZE + c + 1
  }

  function decode(code: number): void {
    let c = code - 1
    codeC = c % SUDOKU_SIZE

    c = Math.floor(c / SUDOKU_SIZE + ebs)
    codeB = c % SUDOKU_SIZE

    c = Math.floor(c / SUDOKU_SIZE + ebs)
    codeA = c
  }
}

我把求解数独的算法封装在了 @algorithm.ts/sudoku 中，下面是求解 $3^{2} \times 3^{2}$ 数独谜题的示例^[5]：

sudoku live demo

// import { SudokuSolver, createSudokuBoardData } from '@algorithm.ts/sudoku'
const solver = new SudokuSolver({ childMatrixWidth: 3 })

render(<SudokuLiveSolver />)

function SudokuLiveSolver() {
  const [puzzle, setPuzzle] = React.useState(() => ([
     5,  0,  6,  7,  8, -1, -1,  3,  2,
     1,  4,  2,  3, -1, -1,  8,  7,  6,
     7,  3, -1, -1, -1, -1,  4,  0, -1,
     8,  7,  3,  0,  6, -1, -1,  4,  1,
     2, -1,  5,  1,  4,  3,  0,  8,  7,
    -1,  1, -1,  8, -1, -1, -1, -1, -1,
    -1,  8,  0, -1, -1, -1,  5,  1,  4,
    -1,  5,  1,  4,  0,  6, -1,  2,  8,
    -1, -1,  7,  5,  1,  8, -1,  6, -1
  ]))

// Resolve sudoku puzzle.
  const solution = React.useMemo(() => {
    // Solve a sudoku puzzle and write the result into the `solution` array.
    const solution = createSudokuBoardData(solver.size)
    solver.solve(puzzle, solution)
    return solution
  }, [puzzle])

return (
    <SudokuSolverDemo
      puzzle={puzzle}
      solution={solution}
      onPuzzleChange={setPuzzle}
    />
  )
}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
// import { SudokuSolver, createSudokuBoardData } from '@algorithm.ts/sudoku'
const solver = new SudokuSolver({ childMatrixWidth: 3 })

render(<SudokuLiveSolver />)

function SudokuLiveSolver() {
  const [puzzle, setPuzzle] = React.useState(() => ([
     5,  0,  6,  7,  8, -1, -1,  3,  2,
     1,  4,  2,  3, -1, -1,  8,  7,  6,
     7,  3, -1, -1, -1, -1,  4,  0, -1,
     8,  7,  3,  0,  6, -1, -1,  4,  1,
     2, -1,  5,  1,  4,  3,  0,  8,  7,
    -1,  1, -1,  8, -1, -1, -1, -1, -1,
    -1,  8,  0, -1, -1, -1,  5,  1,  4,
    -1,  5,  1,  4,  0,  6, -1,  2,  8,
    -1, -1,  7,  5,  1,  8, -1,  6, -1
  ]))

  // Resolve sudoku puzzle.
  const solution = React.useMemo(() => {
    // Solve a sudoku puzzle and write the result into the `solution` array.
    const solution = createSudokuBoardData(solver.size)
    solver.solve(puzzle, solution)
    return solution
  }, [puzzle])

  return (
    <SudokuSolverDemo
      puzzle={puzzle}
      solution={solution}
      onPuzzleChange={setPuzzle}
    />
  )
}

 5  0  6  7  8  4  1  3  2
 1  4  2  3  5  0  8  7  6
 7  3  8  6  2  1  4  0  5
 8  7  3  0  6  5  2  4  1
 2  6  5  1  4  3  0  8  7
 0  1  4  8  7  2  6  5  3
 6  8  0  2  3  7  5  1  4
 3  5  1  4  0  6  7  2  8
 4  2  7  5  1  8  3  6  0

刘汝佳《算法竞赛入门经典──训练指南》 P406 6.3.3 精确覆盖问题和 DLX 算法

前言¶