在 Hackergame 赛后自己去补了点基础，现在写一点 Dancing Links X 算法的学习笔记

算法原理

算法原理参见 Dancing Links - OI Wiki 。

这个算法的原理已经在各种文章中被重述了 114514 次，我并没有信心讲得比他们好，不过还是简要描述一下。

代码实现

代码内包含大量注释（代码基本就是 OI Wiki 的）

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
struct DLX {
    static const int MAXSIZE = 1e5 + 10;
    int n, m, tot, first[MAXSIZE + 10], siz[MAXSIZE + 10];
    int L[MAXSIZE + 10], R[MAXSIZE + 10], U[MAXSIZE + 10], D[MAXSIZE + 10];
    int col[MAXSIZE + 10], row[MAXSIZE + 10];
    int crd, stk[N];

    /*
    建立一个以 column 为元素的循环链表，同时列方向上建立自环。
    编号 0 节点是，列节点的链表头，编号 [1, c] 都是表示列的节点，是一个列方向上的链表头
    行指针 first 指向行方向上的第一个元素
    */
    void build(const int &r, const int &c) {
        n = r, m = c;
        for (int i = 0; i <= c; ++i) {
            L[i] = i - 1, R[i] = i + 1;
            U[i] = D[i] = i;
        } 
        L[0] = c, R[c] = 0, tot = c; 
        memset(first, 0, sizeof(first));
        memset(siz, 0, sizeof(siz));
    }

    /*
    在 r 行 c 列插入元素
    row 和 col 数组储存每一个矩阵元素的行列坐标，siz 储存第 c 列实际元素个数
    首先将新节点插入在列节点 c 的正下方
    随后将新节点插入到行指针 first指向的位置上，注意行方向上也是循环链表结构
    */
    void insert(const int &r, const int &c) {
        col[++tot] = c, row[tot] = r, ++siz[c]; 
        D[tot] = D[c], U[D[c]] = tot, U[tot] = c, D[c] = tot;
        if (!first[r])
            first[r] = L[tot] = R[tot] = tot;
        else {
            R[tot] = R[first[r]], L[R[first[r]]] = tot;
            L[tot] = first[r], R[first[r]] = tot;
        }
    }

    /*
    移除列 c 以及与列 c 连接的所有行上的元素
    remove 操作并不在实际上移除任何元素，而是令两侧的指针跨越被移除元素，通过更新列大小 siz
    被移除的行和列依然能够被对应的行指针和列节点
    */
    void remove(const int &c) {
        int i, j;
        L[R[c]] = L[c], R[L[c]] = R[c];
        for (i = D[c]; i != c; i = D[i])
            for (j = R[i]; j != i; j = R[j])
                U[D[j]] = U[j], D[U[j]] = D[j], --siz[col[j]];
    }

    /*
    恢复列 c 以及与列 c 连接的所有行上的元素
    remove 操作的逆操作
    */
    void recover(const int &c) {
        int i, j;
        for (i = U[c]; i != c; i = U[i])
            for (j = L[i]; j != i; j = L[j])
                U[D[j]] = D[U[j]] = j, ++siz[col[j]];
        L[R[c]] = R[L[c]] = c;
    }

    bool dance(int dep = 1, bool find_one_solution = true) { // 搜索第一个解
        if (!R[0]) {
            crd = dep - 1;
            return find_one_solution;
        }
        int i, j, c = R[0];
        for (i = R[0]; i != 0; i = R[i])
            if (siz[i] < siz[c])
                c = i;
        remove(c);
        for (i = D[c]; i != c; i = D[i]) {
            stk[dep] = row[i];
            for (j = R[i]; j != i; j = R[j]) remove(col[j]);
            if (dance(dep + 1)) return 1;
            for (j = L[i]; j != i; j = L[j]) recover(col[j]);
        }
        recover(c);
        return 0;
    }
} solver;

才发现 Luogu 有一个模板题目 P4929

建模应用

虽然对 DLX 算法的时间复杂度分析是不容易的，但是我们可以注意到问题的复杂度很大程度上由矩阵中 1 的数量来决定。矩阵中的 1 实际上代表了两个不同决策（行）之间通过某种特征（列）产生的排斥关系。因此我们在建模的时候应该尽可能降低 1 节点的数量，让矩阵尽可能小而稀疏。

Luogu P1784 数独

题目地址 P1784

喜闻乐见的模板题。我们首先决定状态也就是决策，这里的决策应该是 在第 r 行第 c 列插入数字 n ，因此总共有 729 种状态

在精确覆盖问题的建模中，需要关注两个决策之间的排斥关系，在本题中有四种不同的排斥关系：

  1. 相同列相同数字存在唯一
  2. 相同行相同数字存在唯一
  3. 相同九宫格相同数字存在唯一
  4. 相同行相同列存在唯一

其中前三种是九宫格的规则比较容易想到，最后一点是每个格子只能填一个数字的天然要求。综合考虑我们有 324 个特征。

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
DLX solver;

int getId(int r, int c, int n) {
    return (r - 1) * 81 + (c - 1) * 9 + n;
}

void insert(int r, int c, int n) {
    int id = getId(r, c, n);
    int room = 3 * ((r - 1) / 3) + (c - 1) / 3 + 1;
    solver.insert(id, 0 * 81 + 9 * (r - 1) + n);
    solver.insert(id, 1 * 81 + 9 * (c - 1) + n);
    solver.insert(id, 2 * 81 + 9 * (room - 1) + n);
    solver.insert(id, 3 * 81 + 9 * (r - 1) + c);
}

int main() {
    int n, m, x, ans[10][10];
    memset(ans, 0, sizeof(ans));
    solver.build(729, 324);
    for (int i = 1; i <= 9; ++i)
        for (int j = 1; j <= 9; ++j) {
            x = read();
            for(int v = 1; v <= 9; ++v)
                if(!x || x == v) insert(i, j, v);
            ans[i][j] = x;
        }
    solver.dance();
    if (solver.crd) {
        for (int i = 1; i <= solver.crd; ++i) {
            int id = solver.ret[i];
            if(!ans[(id - 1) / 81 + 1][(id - 1) % 81 / 9 + 1])
                ans[(id - 1) / 81 + 1][(id - 1) % 81 / 9 + 1] = (id - 1) % 9 + 1;
        }
        for (int i = 1; i <= 9; ++i) {
            for (int j = 1; j <= 9; ++j)
                printf("%d ", ans[i][j]);
            printf("\n");
        }
    } else
        puts("No Solution");
    return 0;
}

Luogu P1074 靶形数独

题目地址 P1074

喜闻乐见的模板题 。

这题跟上一题最大的不同是他需要搜索所有的状态，这就需要我们稍微修改一下 DLX 求解器的部分代码，使得其能够对搜索到的每个解做出相应（我希望写一种能够以给定方式处理每一个搜索到的解的求解器，奈何语言基础太差实在想不出来如何是好）。剩下的部分应该说一模一样，不再详细贴上代码。

N 皇后问题

朴素的搜索可以采用三个状态数组来确认某一个位置是不是可以被放置皇后，我们显然可以使用 DLX 算法来优化这一个过程，只需要将 个皇后的位置全部当成决策，并且将行列、主副对角线总共 个值当成特征就可以了…

并不可以，这个建模有一个显著的问题那就是在 皇后当中对角线并不会被全部占据，只有行列会被完全占据，也就是说这样的建模并不能转化为一个精确覆盖问题。

当然这个问题也是很好解决的，考虑到前 个特征已经足够成为一个精确覆盖问题，我们直接对 DLX 模板内部修改，将算法运行停止条件修改为 R[0] > 2 * n 并且在寻找节点最少的列这一步上面只查找前 列。修改后的 dance 函数如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
bool dance(int dep = 1) {
    if (R[0] > extra) { // extra = 2 * n
        crd = dep - 1;
        return 1;
    }
    int i, j, c = R[0];
    for (i = R[0]; i != 0; i = R[i])
        if (i <= extra && siz[i] < siz[c])
            c = i;
    remove(c);
    for (i = D[c]; i != c; i = D[i]) {
        ret[dep] = row[i];
        for (j = R[i]; j != i; j = R[j]) remove(col[j]);
        if (dance(dep + 1)) return 1;
        for (j = L[i]; j != i; j = L[j]) recover(col[j]);
    }
    recover(c);
    return 0;
}

建模部分很简单就不放了。

Hackergame 2023 小 Z 的谜题

现在回到我学习这个算法的初衷，也就是这个问题的求解。考虑到比赛平台指挥维持几个月的开放，这里以自然语言重新描述一下题目。

在 的立方体内放入如下 个小立方体： 六种规格的方块 分别数量

标答当中直接使用了 Sagemath 内置的模板但是这里我直接用前面的 c++ 模板去写了。

建模部分，考察其中的决策和特征。我们的决策实际上是：在坐标 (x, y, z) 向 (dx, dy, dz) 方向放置编号为 n 的方块，其中方块总共有十六种，决策之间的互斥关系为在空间上不能重叠。那么这里有两种不同的互斥关系：

  • 一个方块放置后，空间内被他占用的格子
  • 一个方块放置后，相同编号的方块不能被放置

分别为二者设计特征（模板内部从 1 开始编号，外部逻辑从 0 开始编号）：

  • 对方块内的每一个坐标 (x, y, z) 的格子，用列 25 * x + 5 * y + z + 1 来编号他
  • 额外开 16 个列代表每一个方块的编号，即 125 + n + 1 列编号

这样我们就得到了一个有 2795 行，141 列的一个待求解矩阵，可以直接使用我们的 DLX 模板来计算。

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
DLX solver;

int num_block = 16;
int count[6] = {3, 4, 2, 2, 2, 3};
int shape[6][3] = {{1, 1, 3}, {1, 2, 2}, {1, 2, 4}, {1, 4, 4}, {2, 2, 2}, {2, 2, 3}};
int dir_[6] = {3, 3, 6, 3, 1, 3};
int dir[16];
int block[16][3];
int idx[17]; 
int id_map[4000][6];

void insertBlock_dir(int x, int y, int z, int cur_idx, int num) {
    int a = 6 - x, b = 6 - y, c = 6 - z;
    for(int i = 0; i < a; ++i) {
        for(int j = 0; j < b; ++j){
            for(int k = 0; k < c; ++k) {
                int id = cur_idx + i * b * c + j * c + k;
                for(int dx = i; dx < i + x; ++dx)
                    for(int dy = j; dy < j + y; ++dy)
                        for(int dz = k; dz < k + z; ++dz)
                            solver.insert(id + 1, dx * 25 + dy * 5 + dz + 1);
                id_map[id + 1][0] = i;
                id_map[id + 1][2] = j;
                id_map[id + 1][4] = k;
                id_map[id + 1][1] = i + x;
                id_map[id + 1][3] = j + y;
                id_map[id + 1][5] = k + z;
                solver.insert(id + 1, 125 + num + 1);
            }
        }
    }
}

void insertBlock(int num) {
    int cur_idx = idx[num], x, y, z;
    for(int i = 0; i < dir[num]; ++i) {
        switch (i) {
            case 0: x = block[num][0]; y = block[num][1]; z = block[num][2]; break;
            case 1: x = block[num][1]; y = block[num][2]; z = block[num][0]; break;
            case 2: x = block[num][2]; y = block[num][0]; z = block[num][1]; break;
            case 3: x = block[num][0]; y = block[num][2]; z = block[num][1]; break;
            case 4: x = block[num][1]; y = block[num][0]; z = block[num][2]; break;
            case 5: x = block[num][2]; y = block[num][1]; z = block[num][0]; break;
        }
        insertBlock_dir(x, y, z, cur_idx, num);
        cur_idx += (6 - block[num][0]) * (6 - block[num][1]) * (6 - block[num][2]);
    }
} 

int main() {
    for(int i = 0, p = 0; i < 6; ++i) {
        for(int j = 0; j < count[i]; ++j) {
            block[p][0] = shape[i][0];
            block[p][1] = shape[i][1];
            block[p][2] = shape[i][2];
            dir[p] = dir_[i];
            ++p;
        }
    } 
    idx[0] = 0;
    for(int i = 0; i < num_block; ++i) 
        idx[i + 1] = idx[i] + dir[i] * (6 - block[i][0])
             * (6 - block[i][1]) * (6 - block[i][2]);
    solver.build(idx[num_block], 125 + num_block);
    for(int i = 0; i < num_block; ++i) 
        insertBlock(i);
    solver.dance();
    for (int i = 1; i <= solver.crd; ++i) {
        int id = solver.ret[i];
        for(int j = 0; j < 6; ++j)
            printf("%d", id_map[id][j]);
    }
    return 0;
}

运行这段代码可以直接得到 flag1 ，程序运行耗时 1.540 s ，使用 O2 编译命令之后来到了 0.768s。

剩下的两个 flag 可以在相同的程序内做出一定的修改来完成求解。

参考资料