Path: blob/master/sage/combinat/matrices/dancing_links_c.h
4079 views
//*****************************************************************************1// Copyright (C) 2008 Carlo Hamalainen <[email protected]>,2//3// Distributed under the terms of the GNU General Public License (GPL)4//5// This code is distributed in the hope that it will be useful,6// but WITHOUT ANY WARRANTY; without even the implied warranty of7// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU8// General Public License for more details.9//10// The full text of the GPL is available at:11//12// http://www.gnu.org/licenses/13//*****************************************************************************1415// This file contains a C++ port of Knuth's C implementation of the Dancing16// Links algorithm. Changes have been minor, including:17//18// * the main search loop doesn't use goto statements19// * there are no hard coded limits20// * the main search algorithm can be called multiple times, and a21// solution is saved in the public variable 'solution' of the22// dancing_links class.2324#include <sstream>25#include <set>26#include <iostream>27#include <vector>28#include <string>29#include <cassert>303132using namespace std;3334template <class T>35string to_string(T x)36{37std::ostringstream stream_out;38stream_out << x;39return stream_out.str();40}4142typedef struct node_struct {43int tag;44struct node_struct *left, *right;45struct node_struct *up, *down;46struct column_struct *col;47} node;4849typedef struct column_struct {50struct node_struct head;51int length;52int index;5354struct column_struct *prev, *next;55} column;5657#define SEARCH_FORWARD 158#define SEARCH_ADVANCE 259#define SEARCH_BACKUP 360#define SEARCH_RECOVER 461#define SEARCH_DONE 56263class dancing_links {64int nr_columns;6566column *root;67column * smallest_column()68{69int minLength = -1;70column *minColumn = 0;7172set<column *> seenColumns;7374for(column *p = root->next; p != root; p = p->next) {75if (minLength == -1 || p->length < minLength) {76minLength = p->length;77minColumn = p;78}79}8081return minColumn;82}8384vector<node*> choice;85vector<column*> col_array;86vector<node*> node_array;8788// For use in search() only.89int mode;90node *current_node;91column *best_col;9293void add_row(vector<int> row, int row_id)94{95node *rowStart = 0;9697for(vector<int>::iterator i = row.begin(); i != row.end(); i++) {98link_entry(*i + 1, row_id, &rowStart);99}100}101102void cover(column *c)103{104// Unlink c from the column list.105column *left = c->prev;106column *right = c->next;107left->next = right;108right->prev = left;109110// For each row that this column points111// to (i.e. has 1's), cover the row too.112for(node *row = c->head.down; row != &(c->head); row = row->down) {113for(node *n = row->right; n != row; n = n->right) {114node *up = n->up;115node *down = n->down;116117up->down = down;118down->up = up;119120n->col->length--;121}122}123}124125126void cover_other_columns(node *c)127{128for(node *p = c->right; p != c; p = p->right)129cover(p->col);130}131132133// Links a row, we call this for each column i that contains134// a 1. We start things by calling with *rowStart == 0135// and this function allocates a new row, otherwise it continues the136// row given.137void link_entry(int i, int row_id, node **rowStart)138{139assert(i <= nr_columns);140141// Link this in to column i.142node *n = (node *) malloc(sizeof(struct node_struct));143node_array.push_back(n);144145n->tag = row_id;146n->col = col_array[i];147148if (*rowStart == 0) *rowStart = n;149150n->down = &(col_array[i]->head);151152if (col_array[i]->length == 0) {153n->up = &(col_array[i]->head);154155col_array[i]->head.down = n;156} else {157n->up = col_array[i]->head.up;158col_array[i]->head.up->down = n;159}160161col_array[i]->head.up = n;162163if (*rowStart != n) {164n->left = (*rowStart)->left;165n->right = (*rowStart);166167(*rowStart)->left->right = n;168(*rowStart)->left = n;169} else {170n->left = n;171n->right = n;172}173174col_array[i]->length++;175}176177178// Exact reverse of cover(c).179void uncover(column *c)180{181for(node *row = c->head.up; row != &(c->head); row = row->up) {182for(node *n = row->left; n != row; n = n->left) {183node *up = n->up;184node *down = n->down;185186up->down = down->up = n;187188n->col->length++;189}190}191192column *left = c->prev;193column *right = c->next;194left->next = right->prev = c;195}196197void uncover_other_columns(node *c)198{199for(node *p = c->left; p != c; p = p->left)200uncover(p->col);201}202203204void save_solution()205{206solution.clear();207208for(vector<node*>::iterator n = choice.begin(); n != choice.end(); n++) {209solution.push_back((*n)->tag);210}211}212213214void setup_columns()215{216assert(nr_columns > 0);217218// The root column, for bookkeeping.219root = (column *) malloc(sizeof(struct column_struct));220root->index = 0;221root->head.tag = 0;222col_array.push_back(root);223224for(int i = 0; i < nr_columns; i++) {225column *current_col = (column *) malloc(sizeof(struct column_struct));226227current_col->length = 0;228current_col->index = i+1;229230current_col->head.tag = -1;231232// The head initially points up/down to itself.233current_col->head.up = &(current_col->head);234current_col->head.down = &(current_col->head);235236root->prev = current_col;237col_array.back()->next = current_col;238239current_col->prev = col_array.back();240current_col->next = root;241242col_array.push_back(current_col);243}244}245public:246vector<int> solution;247248dancing_links()249{250nr_columns = -1;251current_node = NULL;252best_col = NULL;253}254255void add_rows(vector<vector<int> > rows) {256assert(nr_columns == -1);257258// Calculate the maximum value that appears259for(vector<vector<int> >::iterator row = rows.begin(); row != rows.end(); row++) {260for(vector<int>::iterator r = row->begin(); r != row->end(); r++) {261if (nr_columns < *r) nr_columns = *r;262}263}264265nr_columns++;266267setup_columns();268269// Add each row270int row_id = 0;271for(vector<vector<int> >::iterator row = rows.begin(); row != rows.end(); row++) {272add_row(*row, row_id);273row_id++;274}275}276277bool search()278{279assert(nr_columns > 0);280281// If current_node or best_col have changed from being NULL282// then we must have already found a solution and we are283// entering the search() function again, looking for the284// next solution (if any).285if (current_node != NULL || best_col != NULL) {286mode = SEARCH_RECOVER;287} else {288mode = SEARCH_FORWARD;289}290291while(true) {292if (mode == SEARCH_FORWARD) {293best_col = smallest_column();294cover(best_col);295296current_node = best_col->head.down;297choice.push_back(current_node);298299mode = SEARCH_ADVANCE;300}301302if (mode == SEARCH_ADVANCE) {303if (current_node == &(best_col->head)) {304mode = SEARCH_BACKUP; continue;305}306307cover_other_columns(current_node);308309if (col_array[0]->next == col_array[0]) {310save_solution();311return true;312}313mode = SEARCH_FORWARD; continue;314}315316if (mode == SEARCH_BACKUP) {317uncover(best_col);318if ((int) choice.size() == 1) {319mode = SEARCH_DONE; continue;320}321322choice.pop_back();323324current_node = choice.back();325best_col = current_node->col;326327mode = SEARCH_RECOVER;328}329330if (mode == SEARCH_RECOVER) {331uncover_other_columns(current_node);332333choice.pop_back();334current_node = current_node->down;335choice.push_back(current_node);336337mode = SEARCH_ADVANCE; continue;338}339340if (mode == SEARCH_DONE) {341return false;342}343}344}345346void freemem() {347for(vector<column*>::iterator i = col_array.begin(); i != col_array.end(); i++) {348free(*i);349}350351for(vector<node*>::iterator i = node_array.begin(); i != node_array.end(); i++) {352free(*i);353}354}355};356357358359360361362363364365366367368369370371372373374375376377378379