The Dancing Links algorithm from Donald Knuth is a well-known method used to solve sudokus (and various other fun things which can be reduced to an exact cover problem). I recently decided to implement it.

For those who aren't familiar with the algorithm, the main idea is to have a huge matrix of all the candidates where columns and "row dependencies" are taken out the matrix in a recursive manner. Knuth proposed a way of removing and restoring these nodes effectively (without a need of any copy).

I won't go into the details of the algorithm and will assume you have basic understanding of it. A lot of people already manage to popularized the method by writing articles about the naive implementation, so you should be able to document yourself enough to get the idea.

The rest of the article will assume the dimension is `D`

, and `N=D×D`

. For
example, a standard sudoku has `D=3`

and `N=3×3=9`

. Also, the selected language
here is the C.

# A huge matrix

If you look at the initial matrix (empty board) of a `D=3`

sudoku, it looks
like this: Exact Cover Matrix

Each *active* nodes of the matrix is supposed to be double linked horizontally
and vertically to the other *active* nodes.

So in the matrix you will allocate, each node is supposed to have at least 5
pointers: `left`

, `right`

, `up`

, `down`

and the `column`

. Also, they most likely
have a *name* to identify them (in order to print the solution). That name can
be stored in a char, but with the padding it will likely use more.

On a 32 bits machine, each node will occupy something like 24B, and 48B on 64
bits. The matrix has a size of `N×N×N × N×N×4`

, means ```
9×9×9 × 9×9×4 × 24B =
5668704B
```

, or more than `5MB`

. On 64-bits, you need to multiply per 2. And of
course, when you increase the dimension, things get worse. And to be more exact,
you also need a row (`N×N×4`

) for the column headers, so your matrix size would
certainly be `(N×N×N + 1) × N×N×4`

if you do not allocate the header separately.

But this can be greatly reduced with a very simple method.

# Purge the dead nodes

The trick is that the DLX algorithm doesn't require the nodes to be well
aligned in the matrix. What matters is how *active* nodes are **linked**
together (the node selection is based on links and not on the node positions).
The *dead* nodes are never used, so in theory, we can just not allocate them.

Looking at the matrix, you will notice that each column has `N`

actives nodes,
so you can have an array of `N × N×N×4`

nodes. Linking these nodes vertically
is pretty simple: you just have to get down vertically in your array. Getting
the next node horizontally is slightly more tricky; I personally used the size
field in the header to keep track of the last available row (each time a new
node is added in the column, the size is incremented).

In the following schema, the nodes with the same color are double-linked
together horizontally (matrix with `D=2`

):

And with the compact matrix I proposed, it looks like this:

The horizontal links are still illustrated with the color. For the vertical ones, the default matrix layout is kept (you just have to go up and down).

# Compress me the other way around

Another way of doing this is to reduce the matrix not in height, but in width. While this method might get more intuitive (you have one active node per constraint type, so 4 nodes per row), you might have trouble with the column header. Though, this is certainly possible, and I'd be curious to see such implementation.

# Implementation

Here is the code, see the `init()`

function for the matrix initialization. The
other functions mostly follow the naming used in the original paper by Knuth so
you should be able to get a grip on the code easily.

```
#include <limits.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define D 3
#define N (D*D)
struct node {
union {
char size; // headers exclusively use size (to find most interesting column)
uint32_t name; // cells exclusively use name (to print solution)
};
struct node *left, *up, *right, *down;
struct node *column;
};
enum {CST_CELL, CST_ROW, CST_COL, CST_BOX, CST_N};
#define NODES_HEIGHT (N+1) /* +1 for header */
#define NODES_WIDTH (CST_N * N*N)
struct sudansu {
uint32_t solution[N*N];
struct node h;
struct node nodes[NODES_HEIGHT * NODES_WIDTH];
struct node *rows[N*N*N];
};
static unsigned get_node_position(unsigned type, unsigned r, unsigned c, unsigned z)
{
switch (type) {
case CST_CELL: return c + r*N;
case CST_ROW: return z + r*N;
case CST_COL: return z + c*N;
case CST_BOX: return z + c/D*N + r/D*D*N;
default: abort();
}
}
#define CELL_NAME(row, col, numchr) ((row)<<16 | (col)<<8 | (numchr))
static void init(struct sudansu *ctx)
{
unsigned i, j, r, c, z, cst_type;
struct node *x, *row_nodes[CST_N]; // one node per line per constraint type
memset(ctx, 0, sizeof(*ctx));
/* vertical links */
x = ctx->nodes;
for (j = 0; j < NODES_HEIGHT; j++) {
for (i = 0; i < NODES_WIDTH; i++) {
x->up = j == 0 ? &ctx->nodes[(NODES_HEIGHT-1)*NODES_WIDTH + i] : x-NODES_WIDTH;
x->down = j == NODES_HEIGHT-1 ? &ctx->nodes[ 0*NODES_WIDTH + i] : x+NODES_WIDTH;
x->column = &ctx->nodes[i];
x++;
}
}
/* header horizontal links */
ctx->h.right = ctx->nodes;
ctx->h.left = ctx->nodes + NODES_WIDTH-1;
for (i = 0; i < NODES_WIDTH; i++) {
ctx->nodes[i].left = i == 0 ? &ctx->h : &ctx->nodes[i - 1];
ctx->nodes[i].right = i == NODES_WIDTH-1 ? &ctx->h : &ctx->nodes[i + 1];
}
/* nodes horizontal links */
i = 0;
for (r = 0; r < N; r++) {
for (c = 0; c < N; c++) {
for (z = 0; z < N; z++) {
/* locate the CST_N nodes per row */
for (cst_type = 0; cst_type < CST_N; cst_type++) {
unsigned col = get_node_position(cst_type, r, c, z);
struct node *colhead = &ctx->nodes[cst_type*N*N + col];
// locate first unused node of the colum
x = colhead->down + NODES_WIDTH * colhead->size++;
x->name = CELL_NAME(r, c, z+'1');
row_nodes[cst_type] = x;
}
/* link them together */
for (cst_type = 0; cst_type < CST_N; cst_type++) {
row_nodes[cst_type]->left = row_nodes[cst_type == 0 ? CST_N-1 : cst_type-1];
row_nodes[cst_type]->right = row_nodes[cst_type == CST_N-1 ? 0 : cst_type+1];
}
/* index the row since the row positions in ctx->nodes are not linear */
ctx->rows[i++] = row_nodes[0];
}
}
}
}
static void cover(struct node *c)
{
struct node *i, *j;
c->right->left = c->left;
c->left->right = c->right;
for (i = c->down; i != c; i = i->down) {
for (j = i->right; j != i; j = j->right) {
j->down->up = j->up;
j->up->down = j->down;
j->column->size--;
}
}
}
static void uncover(struct node *c)
{
struct node *i, *j;
for (i = c->up; i != c; i = i->up) {
for (j = i->left; j != i; j = j->left) {
j->column->size++;
j->down->up = j->up->down = j;
}
}
c->right->left = c->left->right = c;
}
static int print_solution(const uint32_t *solution)
{
unsigned i, j, k;
char grid[N*N] = {0};
const char *p = grid;
for (i = 0; i < N*N; i++) {
const char row = solution[i]>>16 & 0xff;
const char col = solution[i]>> 8 & 0xff;
const char num = solution[i] & 0xff;
grid[row*N + col] = num;
}
for (j = 0; j < N; j++) {
for (i = 0; i < N; i++) {
printf(" %c", *p ? *p : '?');
if (i+1 != N && (i+1) % D == 0)
printf(" |");
p++;
}
if (j+1 != N && (j+1) % D == 0) {
printf("\n ");
for (k = 0; k < (N+D-1)*2 - 1; k++)
printf("-");
}
printf("\n");
}
printf("\n");
return 1;
}
static struct node *choose_a_column(struct node *h)
{
int size = INT_MAX;
struct node *j, *c = h->right;
if (c == h)
return NULL;
for (j = h->right; j != h; j = j->right)
if (j->size < size)
c = j, size = j->size;
return c;
}
static int search(struct sudansu *ctx, unsigned k)
{
struct node *c, *r, *j;
c = choose_a_column(&ctx->h);
if (!c)
return print_solution(ctx->solution);
cover(c);
for (r = c->down; r != c; r = r->down) {
ctx->solution[k] = r->name;
for (j = r->right; j != r; j = j->right)
cover(j->column);
if (search(ctx, k + 1))
return 1;
ctx->solution[k] = 0;
for (j = r->left; j != r; j = j->left)
uncover(j->column);
}
uncover(c);
return 0;
}
static int parse_grid(struct sudansu *ctx, FILE *stream)
{
char buf[16];
unsigned i, j, k = 0;
for (j = 0; j < N; j++) {
if (!fgets(buf, sizeof buf, stream))
return -1;
for (i = 0; i < N; i++) {
if (buf[i] >= '1' && buf[i] <= '1'+N) {
struct node *x, *r = ctx->rows[j*N*N + i*N + buf[i]-'1'];
ctx->solution[k++] = CELL_NAME(j, i, buf[i]);
cover(r->column);
for (x = r->right; x != r; x = x->right)
cover(x->column);
}
}
}
return k;
}
int main(int ac, char **av)
{
int k = 0;
struct sudansu ctx;
init(&ctx);
if (ac == 2 && strcmp(av[1], "-") == 0) {
k = parse_grid(&ctx, stdin);
if (k < 0)
k = 0;
}
search(&ctx, k);
return 0;
}
```