Algorithm to the exact cover problem: Dancing Links

2021/01/25

Introduction:

舞蹈链(Dancing Links), 也叫 DLX, 是由 Donald Knuth 提出的数据结构。Knuth 归功于 Hiroshi Hitotsumatsu 与 Kōhei Noshita 在 1979 的研究1, 但是 Knuth 的论文2让舞蹈链流行。

 

首先在学术界已经证明“覆盖问题”是 NPC 问题,没有多项式复杂度的算法可以解决,比如最经典的问题当属所谓的“N皇后问题”以及“数独”了,清晰记得大概半年前还不会 DLX,数独问题只能对着《算法竞赛进阶指南》中晦涩又复杂的搜索和剪枝技巧望而却步。而且用朴素的搜索来写意味着需要自己构造一个缜密的搜索思路和搜索框架,以及最优的剪枝技巧,这是很困难的。如果用 DLX 来解决数独问题,也许只需要把 DLX 的代码模板交上去就能 AC 了。

 

覆盖问题分为“精确覆盖问题”和“重复覆盖问题”,避免篇幅过长,这篇博客只介绍精确覆盖问题以及代码模板的应用。之后的博客中会补充重复覆盖问题。

精确覆盖问题:

给定一个 \(N×M\) 的数字矩阵 \(A\),矩阵中的元素 \(A_{i,j} \in \lbrace 0,1 \rbrace\)。

请问,你能否在矩阵中找到一个行的集合,使得这些行中,每一列都有且仅有一个数字 \(1\)。

输入格式:

第一行包含两个整数 \(N\) 和 \(M\)。

接下来 \(N\) 行,每行包含 \(M\) 个整数(\(0\) 或 \(1\)),表示完整的数字矩阵。

输出格式:

如果能找到满足条件的行的集合,则在一行中依次输出这些行的编号(行编号 \(1 \sim N\))。

如果方案不唯一,则以任意顺序输出任意方案即可。否则,输出 No Solution!

数据范围:

\(1 \le N,M \le 500\) \( \ \) 时/空限制: \(1s / 64MB\)

数据保证矩阵中 \(1\) 的数量不超过 \(5000\)。

输入样例1:

3 3
0 1 0
0 0 1
1 0 0

输出样例1:

1 2 3

输入样例2:

4 4
0 0 0 1
1 0 0 0
1 1 0 1
0 1 0 0

输出样例2:

No Solution!

 

接下来介绍使用 \(DLX\) 搜索的算法流程:

首先是建图

DLX 建图之后的状态(图源OI-WIKI)

其次是搜索框架

补充说明:

DLX模板:

#include <bits/stdc++.h>
using namespace std;

const int N = 5510; //所有1的数量加上第0行的head表头 注意这个模板没有列的表头

int n, m;
//l每个点的左边的节点 r右边的节点 d下边的 u上边的 s该列有多少个1 row点的行坐标 col列坐标
int l[N], r[N], u[N], d[N], s[N], row[N], col[N], idx;
int ans[N], top;

void init() { //初始化表头 表头一共以m + 1个点,最左边一个是哨兵 比列数多一列
    for (int i = 0; i <= m; i++) {
        l[i] = i - 1, r[i] = i + 1;
        u[i] = d[i] = i;
    }
    l[0] = m, r[m] = 0; //循环链表
    idx = m + 1; //点的内存分配指针
}

void add(int &hh, int &tt, int x, int y) { //在坐标(x,y)插入一个点
    row[idx] = x, col[idx] = y, s[y]++;
    //每次都是插到第0行表头的下边一行 是头插
    //hh tt每一行最左边2个指针,插入到这2个指针中间
    u[idx] = y, d[idx] =  d[y], u[d[y]] = idx, d[y] = idx;
    r[hh] = l[tt] = idx, r[idx] = tt, l[idx] = hh;
    tt = idx ++; //tt向左移动到新加入的点的位置
}

void remove(int p) { //删除第p列 以及该列有1的所有的行
    r[l[p]] = r[p], l[r[p]] = l[p];  //删除列只需要删除表头即可
    //下边的这个二重循环,并没有删除第p列的元素,只删除了表头
    for (int i = d[p]; i != p; i = d[i])
        for (int j = r[i]; j != i; j = r[j]) {
            //删除一行 修改该行的每个元素的上下指针即可
            s[col[j]] --;
            u[d[j]] = u[j], d[u[j]] = d[j];
        }
}

void resume(int p) { //恢复一次remove的所有删除
    for (int i = u[p]; i != p; i = u[i])
        for (int j = l[i]; j != i; j = l[j]) {
            u[d[j]] = j, d[u[j]] = j;
            s[col[j]] ++;
        }
    r[l[p]] = p, l[r[p]] = p;
}

bool dfs() {
    if (!r[0]) return true; //第0行链表头右边已经被删完了 说明所有的列都已经被覆盖
    int p = r[0];
    for (int i = r[0]; i; i = r[i]) //找出1最少的一列 注意不要混淆r[]和row[]
        if (s[i] < s[p])
            p = i;
    remove(p); //先删去这个表头 表示这一列已经被覆盖 因为搜的时候只会找表头

    for (int i = d[p]; i != p; i = d[i]) { //遍历当前1最少的列的每一个1所在的行,每一行都横向遍历删除列
        ans[++top] = row[i]; //记录答案
        //删除该行除了当前被覆盖的列以外,其余所有1所在的列以及该列有1的所有行
        for (int j = r[i]; j != i; j = r[j]) remove(col[j]);
        if (dfs()) return true;
        for (int j = l[i]; j != i; j = l[j]) resume(col[j]);
        top --;
    }

    resume(p); //恢复现场 恢复表头
    return false;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    cin >> n >> m;
    init();
    for (int i = 1; i <= n;  i++) { //读取所有点
        int hh = idx, tt = idx; //头插指针
        for (int j = 1; j <= m; j++) {
            int x;
            cin >> x;
            if (x) add(hh, tt, i, j); //是1就插入
        }
    }

    if (dfs()) {
        for (int i = 1; i <= top; i++) cout << ans[i] << ' ';
        cout << endl;
    }
    else cout << "No Solution!" << endl;

    return 0;
}

 

DLX解决数独问题:

请你将一个 \(16 \times 16\) 的数独填写完整,使得每行、每列、每个 \(4 \times 4\)十六宫格内字母 \(A \sim P\) 均恰好出现一次。

保证每个输入只有唯一解决方案。

输入格式:

输入包含多组测试用例。

每组测试用例包括 \(16\) 行,每行一组字符串,共 \(16\) 个字符串。

第i个字符串表示数独的第 \(i\) 行。

字符串包含字符可能为字母\(A \sim P\) 或 \(“-”\)(表示等待填充)。

测试用例之间用单个空行分隔,输入至文件结尾处终止。

输出格式:

对于每个测试用例,均要求保持与输入相同的格式,将填充完成后的数独输出。

每个测试用例输出结束后,输出一个空行。

输入样例:

--A----C-----O-I
-J--A-B-P-CGF-H-
--D--F-I-E----P-
-G-EL-H----M-J--
----E----C--G---
-I--K-GA-B---E-J
D-GP--J-F----A--
-E---C-B--DP--O-
E--F-M--D--L-K-A
-C--------O-I-L-
H-P-C--F-A--B---
---G-OD---J----H
K---J----H-A-P-L
--B--P--E--K--A-
-H--B--K--FI-C--
--F---C--D--H-N-

输出样例:

FPAHMJECNLBDKOGI
OJMIANBDPKCGFLHE
LNDKGFOIJEAHMBPC
BGCELKHPOFIMAJDN
MFHBELPOACKJGNID
CILNKDGAHBMOPEFJ
DOGPIHJMFNLECAKB
JEKAFCNBGIDPLHOM
EBOFPMIJDGHLNKCA
NCJDHBAEKMOFIGLP
HMPLCGKFIAENBDJO
AKIGNODLBPJCEFMH
KDEMJIFNCHGAOPBL
GLBCDPMHEONKJIAF
PHNOBALKMJFIDCEG
IAFJOECGLDPBHMNK

思路:

考虑如何将数独转换为精确覆盖问题,仔细思考一下覆盖问题的本质,就是选择一些当前已经有的条件,然后满足一个条件。 对于数独而言,就是“选择一些格子的填法,来满足数独的约束条件”,链表的行是已经有的条件,链表的列是待满足的条件。

我们用链表的所有行代表所有格子可能填的情况,一共有 \(16 \times 16\) 个格子,每个格子都能填 \(16\) 种情况,也就有 \(16 \times 16 \times 16\) 行。

列就是所有的约束条件,一共有 \(4\) 个约束条件:

  1. 每个格子只能填一个数,一共 \(16 \times 16\) 列
  2. 每行都是 \(0 \sim 15\) 只填一个,同样 \(16 \times 16\) 列
  3. 列和条件 \(2\) 相同, 也是 \(16 \times 16\) 列
  4. 每个宫也是只能 \(0 \sim\ 15\) 各填一个,同样 \(16 \times 16\) 列

所以有 \(4096\) 行,\(1024\) 列,但是 \(1\) 的个数是稀疏的,建图时我们要计算最多的可能的 \(1\) 的个数。 只考虑每一行中 \(1\) 可能出现的最多次数即可,每一行代表一个格子选择的字母,因为只是“一个格子,一个字母”,所以在 \(4\) 种列的约束条件中,每一种都最多只会出现一次。举例子就是一个格子,不可能在“每个格子恰好填一个数”这一类中出现 \(2\)次,也不可能在“每行 \(0 \sim 15\) 恰好出现一次”中出现 \(2\) 次,其余的约束条件同理。 所以一行中“最多”只会出现 \(4\) 次 \(1\),那么最多的 \(1\)的点个数就是 \(4096 \times 4\),加上表头的 \(1024 + 1 \)个点,大约是 \(17000\) 多个点,开 \(18000\) 即可。

代码:

#include <bits/stdc++.h>
using namespace std;

const int N = 18000; //每行最多4个点,4096+1024+1个1最多

int m = 16 * 16 * 4; //列数
int u[N], d[N], l[N], r[N], s[N], col[N], row[N], idx;
int ans[N], top;

struct Op {
    int x, y; //坐标
    char z; //该坐标填的字母
} op[N];

char g[20][20];

void init() {
    for (int i = 0; i <= m; i++) {
        l[i] = i - 1, r[i] = i + 1;
        s[i] = 0;
        d[i] = u[i] = i;
    }
    l[0] = m, r[m] = 0;
    idx = m + 1;
}

void add(int &hh, int &tt, int x, int y) {
    row[idx] = x, col[idx] = y, s[y] ++;
    u[idx] = y, d[idx] = d[y], u[d[y]] = idx, d[y] = idx;
    r[hh] = l[tt] = idx, r[idx] = tt, l[idx] = hh;
    tt = idx ++;
}

void remove(int p) {
    r[l[p]] = r[p], l[r[p]] = l[p];
    for (int i = d[p]; i != p; i = d[i])
        for (int j = r[i]; j != i; j = r[j]) {
            s[col[j]] --;
            u[d[j]] = u[j], d[u[j]] = d[j];
        }
}

void resume(int p) {
    for (int i = u[p]; i != p; i = u[i])
        for (int j = l[i]; j != i; j = l[j]) {
            u[d[j]] = j, d[u[j]] = j;
            s[col[j]] ++;
        }
    r[l[p]] = p, l[r[p]] = p;
}

bool dfs() {
    if (!r[0]) return true;
    int p = r[0];
    for (int i = r[0]; i; i = r[i])
        if (s[i] < s[p])
            p = i;
    remove(p);
    for (int i = d[p]; i != p; i = d[i])
    {
        ans[ ++ top] = row[i];
        for (int j = r[i]; j != i; j = r[j]) remove(col[j]);
        if (dfs()) return true;
        for (int j = l[i]; j != i; j = l[j]) resume(col[j]);
        top -- ;
    }
    resume(p);
    return false;
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(0);

    while (cin >> g[0]) {
        for (int i = 1; i < 16; i++) cin >> g[i];
        init();

        for (int i = 0, n = 1; i < 16; i++) //n是行号
            for (int j = 0; j < 16; j++) {
                int a = 0, b = 15;
                if (g[i][j] != '-') a = b = g[i][j] - 'A';//一个位置已经填好了那么就让a和b都是已经填的字母
                for (int k = a; k <= b; k++, n ++) { //枚举该位置填什么
                    int hh = idx, tt = idx;
                    op[n] = {i, j, k + 'A'};
                    add(hh, tt, n, i * 16 + j + 1); //加1是因为行列下标都从1开始
                    add(hh, tt, n, 256 + i * 16 + k + 1);
                    add(hh, tt, n, 256 * 2 + j * 16 + k + 1);
                    add(hh, tt, n, 256 * 3 + (i / 4 * 4 + j / 4) * 16 + k + 1); //每个宫,宫从0开始到15
                }
            }
        dfs();
        for (int i = 1; i <= top; i++) {
            auto t = op[ans[i]];
            g[t.x][t.y] = t.z;
        }
        for (int i = 0; i < 16; i++) puts(g[i]);
        puts("");
    }
    return 0;
}

 

扩展题目:

参考文献: