Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
sagemath
GitHub Repository: sagemath/sagelib
Path: blob/master/sage/combinat/matrices/dancing_links_c.h
4079 views
1
//*****************************************************************************
2
// Copyright (C) 2008 Carlo Hamalainen <[email protected]>,
3
//
4
// Distributed under the terms of the GNU General Public License (GPL)
5
//
6
// This code is distributed in the hope that it will be useful,
7
// but WITHOUT ANY WARRANTY; without even the implied warranty of
8
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
9
// General Public License for more details.
10
//
11
// The full text of the GPL is available at:
12
//
13
// http://www.gnu.org/licenses/
14
//*****************************************************************************
15
16
// This file contains a C++ port of Knuth's C implementation of the Dancing
17
// Links algorithm. Changes have been minor, including:
18
//
19
// * the main search loop doesn't use goto statements
20
// * there are no hard coded limits
21
// * the main search algorithm can be called multiple times, and a
22
// solution is saved in the public variable 'solution' of the
23
// dancing_links class.
24
25
#include <sstream>
26
#include <set>
27
#include <iostream>
28
#include <vector>
29
#include <string>
30
#include <cassert>
31
32
33
using namespace std;
34
35
template <class T>
36
string to_string(T x)
37
{
38
std::ostringstream stream_out;
39
stream_out << x;
40
return stream_out.str();
41
}
42
43
typedef struct node_struct {
44
int tag;
45
struct node_struct *left, *right;
46
struct node_struct *up, *down;
47
struct column_struct *col;
48
} node;
49
50
typedef struct column_struct {
51
struct node_struct head;
52
int length;
53
int index;
54
55
struct column_struct *prev, *next;
56
} column;
57
58
#define SEARCH_FORWARD 1
59
#define SEARCH_ADVANCE 2
60
#define SEARCH_BACKUP 3
61
#define SEARCH_RECOVER 4
62
#define SEARCH_DONE 5
63
64
class dancing_links {
65
int nr_columns;
66
67
column *root;
68
column * smallest_column()
69
{
70
int minLength = -1;
71
column *minColumn = 0;
72
73
set<column *> seenColumns;
74
75
for(column *p = root->next; p != root; p = p->next) {
76
if (minLength == -1 || p->length < minLength) {
77
minLength = p->length;
78
minColumn = p;
79
}
80
}
81
82
return minColumn;
83
}
84
85
vector<node*> choice;
86
vector<column*> col_array;
87
vector<node*> node_array;
88
89
// For use in search() only.
90
int mode;
91
node *current_node;
92
column *best_col;
93
94
void add_row(vector<int> row, int row_id)
95
{
96
node *rowStart = 0;
97
98
for(vector<int>::iterator i = row.begin(); i != row.end(); i++) {
99
link_entry(*i + 1, row_id, &rowStart);
100
}
101
}
102
103
void cover(column *c)
104
{
105
// Unlink c from the column list.
106
column *left = c->prev;
107
column *right = c->next;
108
left->next = right;
109
right->prev = left;
110
111
// For each row that this column points
112
// to (i.e. has 1's), cover the row too.
113
for(node *row = c->head.down; row != &(c->head); row = row->down) {
114
for(node *n = row->right; n != row; n = n->right) {
115
node *up = n->up;
116
node *down = n->down;
117
118
up->down = down;
119
down->up = up;
120
121
n->col->length--;
122
}
123
}
124
}
125
126
127
void cover_other_columns(node *c)
128
{
129
for(node *p = c->right; p != c; p = p->right)
130
cover(p->col);
131
}
132
133
134
// Links a row, we call this for each column i that contains
135
// a 1. We start things by calling with *rowStart == 0
136
// and this function allocates a new row, otherwise it continues the
137
// row given.
138
void link_entry(int i, int row_id, node **rowStart)
139
{
140
assert(i <= nr_columns);
141
142
// Link this in to column i.
143
node *n = (node *) malloc(sizeof(struct node_struct));
144
node_array.push_back(n);
145
146
n->tag = row_id;
147
n->col = col_array[i];
148
149
if (*rowStart == 0) *rowStart = n;
150
151
n->down = &(col_array[i]->head);
152
153
if (col_array[i]->length == 0) {
154
n->up = &(col_array[i]->head);
155
156
col_array[i]->head.down = n;
157
} else {
158
n->up = col_array[i]->head.up;
159
col_array[i]->head.up->down = n;
160
}
161
162
col_array[i]->head.up = n;
163
164
if (*rowStart != n) {
165
n->left = (*rowStart)->left;
166
n->right = (*rowStart);
167
168
(*rowStart)->left->right = n;
169
(*rowStart)->left = n;
170
} else {
171
n->left = n;
172
n->right = n;
173
}
174
175
col_array[i]->length++;
176
}
177
178
179
// Exact reverse of cover(c).
180
void uncover(column *c)
181
{
182
for(node *row = c->head.up; row != &(c->head); row = row->up) {
183
for(node *n = row->left; n != row; n = n->left) {
184
node *up = n->up;
185
node *down = n->down;
186
187
up->down = down->up = n;
188
189
n->col->length++;
190
}
191
}
192
193
column *left = c->prev;
194
column *right = c->next;
195
left->next = right->prev = c;
196
}
197
198
void uncover_other_columns(node *c)
199
{
200
for(node *p = c->left; p != c; p = p->left)
201
uncover(p->col);
202
}
203
204
205
void save_solution()
206
{
207
solution.clear();
208
209
for(vector<node*>::iterator n = choice.begin(); n != choice.end(); n++) {
210
solution.push_back((*n)->tag);
211
}
212
}
213
214
215
void setup_columns()
216
{
217
assert(nr_columns > 0);
218
219
// The root column, for bookkeeping.
220
root = (column *) malloc(sizeof(struct column_struct));
221
root->index = 0;
222
root->head.tag = 0;
223
col_array.push_back(root);
224
225
for(int i = 0; i < nr_columns; i++) {
226
column *current_col = (column *) malloc(sizeof(struct column_struct));
227
228
current_col->length = 0;
229
current_col->index = i+1;
230
231
current_col->head.tag = -1;
232
233
// The head initially points up/down to itself.
234
current_col->head.up = &(current_col->head);
235
current_col->head.down = &(current_col->head);
236
237
root->prev = current_col;
238
col_array.back()->next = current_col;
239
240
current_col->prev = col_array.back();
241
current_col->next = root;
242
243
col_array.push_back(current_col);
244
}
245
}
246
public:
247
vector<int> solution;
248
249
dancing_links()
250
{
251
nr_columns = -1;
252
current_node = NULL;
253
best_col = NULL;
254
}
255
256
void add_rows(vector<vector<int> > rows) {
257
assert(nr_columns == -1);
258
259
// Calculate the maximum value that appears
260
for(vector<vector<int> >::iterator row = rows.begin(); row != rows.end(); row++) {
261
for(vector<int>::iterator r = row->begin(); r != row->end(); r++) {
262
if (nr_columns < *r) nr_columns = *r;
263
}
264
}
265
266
nr_columns++;
267
268
setup_columns();
269
270
// Add each row
271
int row_id = 0;
272
for(vector<vector<int> >::iterator row = rows.begin(); row != rows.end(); row++) {
273
add_row(*row, row_id);
274
row_id++;
275
}
276
}
277
278
bool search()
279
{
280
assert(nr_columns > 0);
281
282
// If current_node or best_col have changed from being NULL
283
// then we must have already found a solution and we are
284
// entering the search() function again, looking for the
285
// next solution (if any).
286
if (current_node != NULL || best_col != NULL) {
287
mode = SEARCH_RECOVER;
288
} else {
289
mode = SEARCH_FORWARD;
290
}
291
292
while(true) {
293
if (mode == SEARCH_FORWARD) {
294
best_col = smallest_column();
295
cover(best_col);
296
297
current_node = best_col->head.down;
298
choice.push_back(current_node);
299
300
mode = SEARCH_ADVANCE;
301
}
302
303
if (mode == SEARCH_ADVANCE) {
304
if (current_node == &(best_col->head)) {
305
mode = SEARCH_BACKUP; continue;
306
}
307
308
cover_other_columns(current_node);
309
310
if (col_array[0]->next == col_array[0]) {
311
save_solution();
312
return true;
313
}
314
mode = SEARCH_FORWARD; continue;
315
}
316
317
if (mode == SEARCH_BACKUP) {
318
uncover(best_col);
319
if ((int) choice.size() == 1) {
320
mode = SEARCH_DONE; continue;
321
}
322
323
choice.pop_back();
324
325
current_node = choice.back();
326
best_col = current_node->col;
327
328
mode = SEARCH_RECOVER;
329
}
330
331
if (mode == SEARCH_RECOVER) {
332
uncover_other_columns(current_node);
333
334
choice.pop_back();
335
current_node = current_node->down;
336
choice.push_back(current_node);
337
338
mode = SEARCH_ADVANCE; continue;
339
}
340
341
if (mode == SEARCH_DONE) {
342
return false;
343
}
344
}
345
}
346
347
void freemem() {
348
for(vector<column*>::iterator i = col_array.begin(); i != col_array.end(); i++) {
349
free(*i);
350
}
351
352
for(vector<node*>::iterator i = node_array.begin(); i != node_array.end(); i++) {
353
free(*i);
354
}
355
}
356
};
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379