Exact cover III

,

The story so far: I wrote a solver for a generic exact cover problem in Python, and used it to find polyomino tilings and Sudoku solutions. Now I’m going to look at a couple of problems that don’t translate so easily into exact cover instances.

1. The eight queens problem

The eight queens problem is to place eight queens on a chessboard so that no queen is attacked by any other. (That is, so that there is at most one queen in each row, column, and diagonal.) The problem is well-known to computer programmers because it was used as the main example in Dijkstra’s influential ‘Notes on Structured Programming’ (1970), where he shows how to solve it using a backtracking algorithm, developing the program through a series of stepwise refinements.

Instead of using structured programming, let’s solve it by translating it to an exact cover instance. As in the previous examples, I can specify an encoding of constraints and choices as Python objects:

ConstraintPython representation
Queen on the \(i\)th row.The tuple ('row', i)
Queen on the \(j\)th column.The tuple ('column', j)
Queen on the \(k\)th leading diagonal.The tuple ('leading', k)
Queen on the \(k\)th trailing diagonal.The tuple ('trailing', k)
ChoicePython representation
Place a queen at \((i, j)\).The tuple ('queen', i, j)

Here’s a function that encodes the \(n\)-queens problem as an exact cover instance using these choices and constraints:

from exactcover import ExactCover
from itertools import product

def queens(n=8, random=False):
    """Return an iterator that yields the solutions to the n-queens
    problem.

    Arguments:
    n -- The size of the problem. (Default: 8)
    random -- Generate solutions in random order? (Default: False.)

    """
    constraints = {
        ('queen', i, j): (('row', i), ('column', j),
                          ('leading', i - j + n - 1), ('trailing', i + j))
        for i, j in product(range(n), repeat=2)
    }
    return ExactCover(constraints, random=random)

But this does not work!

>>> next(queens())
Traceback (most recent call last):
  File "", line 1, in 
  File "exactcover.py", line 70, in __next__
    return next(self.iter)
StopIteration

Why is that? Tracing the execution of the exact cover solver, the first thing that happens is that it decides to satisfy the constraint ('leading', 0). The only choice that satisfies this constraint is ('queen', 0, 7). This choice also satisfies the constraint ('row', 0), which causes the choice ('queen', 0, 0) to be eliminated (since it would be in the same row). But now there are no choices that satisfy the constraint ('trailing', 0) and further progress is impossible.

In an exact cover problem, all the constraints must be satisfied. But in this problem there aren’t enough queens to cover all the diagonals: on an 8×8 board there are 15 leading diagonals and 8 queens, but each queen can only cover one leading diagonal. So 7 leading diagonal constraints must go unsatisfied (and similarly for trailing diagonals).

I can get around this difficulty by adding additional choices whose only purpose is to ensure that all the constraints can be satisfied:

ChoicePython representation
Cover the \(k\)th leading diagonal.The tuple ('leading', k)
Cover the \(k\)th trailing diagonal.The tuple ('trailing', k)

The idea is that after placing the eight queens, the remaining 14 diagonal constraints can be satisfied using 14 of these extra choices. They will appear in the solution, but they are redundant and can be ignored. Artificial devices like these that help to translate one problem into another are called ‘gadgets’. With these gadgets in hand, it’s easy to construct the exact cover instances:

from exactcover import ExactCover
from itertools import product

def queens(n=8, random=False):
    """Return an iterator that yields the solutions to the n-queens
    problem.

    Arguments:
    n -- The size of the problem. (Default: 8)
    random -- Generate solutions in random order? (Default: False.)

    """
    def constraints():
        for i, j in product(range(n), repeat=2):
            yield ('queen', i, j), (('row', i),
                                    ('column', j),
                                    ('leading', i - j + n - 1),
                                    ('trailing', i + j))
        for k in range(2 * n - 1):
            yield ('leading', k), (('leading', k),)
            yield ('trailing', k), (('trailing', k),)
    return ExactCover(dict(constraints()), random=random)

Here’s a simple decoder for the solutions. Note that the decoder ignores the gadgets (the choices representing unattacked diagonals) and only looks at the positions of the queeens.

def decode_queens(n, solution):
    grid = [['+'] * n for _ in range(n)]
    for choice in solution:
        if choice[0] == 'queen':
            _, i, j = choice[:3]
            grid[i][j] = '♛'
    return '\n'.join(map(''.join, grid))

Here’s the solver in operation:

>>> print(decode_queens(12, next(queens(12))))
++++++♛+++++
++++++++♛+++
+++++♛++++++
++♛+++++++++
+++++++++♛++
+++♛++++++++
++++++++++♛+
+++++++♛++++
++++♛+++++++
+♛++++++++++
+++++++++++♛
♛+++++++++++

The solver can count the solutions for small board sizes (sequence A000170 in the OEIS):

>>> [sum(1 for _ in queens(n)) for n in range(1, 11)]
[1, 0, 0, 2, 10, 4, 40, 92, 352, 724]

Here’s code for drawing solutions as SVG images, using the svgwrite package:

from svgwrite import Drawing

def solutions_svg(n, solutions, filename, columns=10, size=20, margin=10,
                  fg_colour='black', bg_colour='#eee'):
    """Format n-queens solutions as an SVG image.

    Required arguments:
    n -- size of board.
    solutions -- iterable of solutions to the n-queens problem, each
        of which is a sequence of choices, some of which are tuples
        ('queen', i, j) giving the locations of the queens.
    filename -- where to save the SVG drawing.

    Optional arguments:
    columns -- number of solutions per row (default: 1).
    size -- width and height of each board square (default: 20).
    margin -- margin between boards (default: 10).
    fg_colour -- colour of queens and even squares (default: black).
    bg_colour -- colour of odd squares (default: #eee).

    """
    solutions = list(solutions)
    rows = (len(solutions) + columns - 1) // columns
    drawing_size = (columns * (n * size + margin) - margin,
                    rows    * (n * size + margin) - margin)
    drawing = Drawing(debug=False, filename=filename, size=drawing_size)
    square_symbol = drawing.symbol(id='square')
    drawing.add(square_symbol)
    square_path = drawing.path(stroke=fg_colour, stroke_width=.25)
    square_symbol.add(square_path)
    for i in range(0, size, 3):
        square_path.push(['M', size - i, 0, 'L', 0, size - i])
        if i: square_path.push(['M', i, size, 'L', size, i])
    board_symbol = drawing.symbol(id='board')
    drawing.add(board_symbol)
    board_symbol.add(drawing.rect((0, 0), (n * size, n * size), fill=bg_colour))
    for i, j in product(range(n), repeat=2):
        if (i + j) % 2:
            board_symbol.add(drawing.use('#square', x=i * size, y=j * size))
    queen_symbol = drawing.symbol(id='q')
    drawing.add(queen_symbol)
    queen_symbol.add(drawing.text('♛', x=[size // 10], y=[size * 3 // 4],
                                  style='color:{} font-size:{}'.format(
                                      fg_colour, (size * 4 // 5))))
    for i, solution in enumerate(solutions):
        y, x = divmod(i, columns)
        x = x * (size * n + margin)
        y = y * (size * n + margin)
        drawing.add(drawing.use("#board", x=x, y=y))
        for choice in solution:
            if choice[0] == 'queen':
                _, i, j = choice[:3]
                drawing.add(drawing.use('#q', x=x + i * size, y=y + j * size))
    drawing.save()

And here are the 92 solutions on the 8×8 board:

2. The sixteen queens problem

The sixteen queens problem is to place sixteen queens on a chessboard so that there are no more than two queens in each row, column and diagonal.

(This problem is a simplification of puzzle 317 in Dudeney’s Amusements in Mathematics (1917), which asks:

Place two pawns in the middle of the chessboard, one at Q4 and the other at K5. Now, place the remaining fourteen pawns (sixteen in all) so that no three shall be in a straight line in any possible direction. Note that I purposely do not say queens, because by the words "any possible direction" I go beyond attacks on diagonals. The pawns must be regarded as mere points in space—at the centres of the squares.

This is known as the no-three-in-line problem, but the extra oblique lines make the problem too complicated to use as an example here.)

How can I represent this as an exact cover instance? Consider the row constraints: there must be two queens in each row. So I can have a constraint for the first queen in a row and another constraint for the second queen in a row. But when choosing to place a queen at \((i, j)\), it might be the first queen or the second queen in the row. So there had better be two choices: “place a queen at \((i, j)\) that’s the first queen in row \(i\)” and “place a queen at \((i, j)\) that’s the second queen in row \(i\)”. A similar analysis applies to the columns, to the leading diagonals, and to the trailing diagonals. So each constraint splits into 2, and each choice splits into 2×2×2×2 = 16:

ConstraintPython representation
\(a\)th queen on the \(i\)th row.The tuple ('row', i, a)
\(b\)th queen on the \(j\)th column.The tuple ('column', j, b)
\(c\)th queen on the \(k\)th leading diagonal.The tuple ('leading', k, c)
\(d\)th queen on the \(k\)th trailing diagonal.The tuple ('trailing', k, d)
ChoicePython representation
Place a queen at \((i, j)\) that’s the \(a\)th queen in its row, the \(b\)th in its column, the \(c\)th in its leading diagonal, and the \(d\)th in its trailing diagonal.The tuple ('queen', i, j, a, b, c, d)
Omit the \(a\)th queen on the \(k\)th leading diagonal.The tuple ('leading', k, a)
Omit the \(a\)th queen on the \(k\)th trailing diagonal.The tuple ('trailing', k, a)
def queens(n=8, random=False):
    """Return an iterator that yields the solutions to the 2n-queen problem.

    Arguments:
    n -- The size of the problem. (Default: 8)
    random -- Generate solutions in random order? (Default: False.)

    """
    def constraints():
        for i, j, a, b, c, d in product(*map(range, (n, n, 2, 2, 2, 2))):
            yield ('queen', i, j, a, b, c, d), (
                ('row', i, a),
                ('column', j, b),
                ('leading', i - j + n - 1, c),
                ('trailing', i + j, d))
        for k, a in product(range(2 * n - 1), range(2)):
            yield ('leading', k, a), (('leading', k, a),)
            yield ('trailing', k, a), (('trailing', k, a),)
    return ExactCover(dict(constraints()), random=random)

This works, but it’s slower than a sleepy sloth, taking minutes to emit the first solution. The reason for the poor performance is the combinatorial explosion caused by the gadgets. It’s one thing to use a gadget in order to show that one problem can be reduced to another in a proof that a problem belongs to a complexity class. It’s quite another when you actually want to compute some solutions. It would make much more sense to modify the exact cover solver.

3. Multiset cover

There are various ways that the solver might be modified, but here I’m going to make two changes: first, represent the constraints as a multiset; and second, only require some of the constraints to be satisfied.

from collections import Counter, Iterator, defaultdict
from random import shuffle

class MultiCover(Iterator):
    """An iterator that yields solutions to an EXACT COVER problem with
    multiplicity.

    An EXACT COVER problem consists of a set of "choices". Each choice
    satisfies one or more "constraints". Choices and constraints may
    be represented by any hashable Python objects, with the proviso
    that all choices must be distinct, as must all constraints.

    Some of the constraints are required (must be satisfied) and some
    are optional (may be satisfied). Each constraint c has a
    multiplicity counts[c].

    A solution is a list of choices such that each of the
    "must-satisfy" constraints is satisfied by exactly counts[c] of
    the choices in the solution.

    Required argument:
    constraints -- A map from each choice to an iterable of
        constraints satisfied by that choice.

    Optional arguments:
    initial -- An iterable of choices that must appear in every
        solution. (Default: no choices.)
    random -- Generate solutions in random order? (Default: False.)
    counts -- Map from constraint to the number of times that
        constraint must be satisfied. (By default all constraints must
        be satisfied exactly once.)
    satisfy -- Iterable of contraints that must be satisfied. (Default:
        all constraints.)

    For example:

        >>> data = dict(A = [1, 4, 7],
        ...             B = [1, 4],
        ...             C = [4, 5, 7],
        ...             D = [3, 5, 6],
        ...             E = [2, 3, 6, 7],
        ...             F = [2, 7])
        >>> sorted(next(MultiCover(data)))
        ['B', 'D', 'F']

    By passing a counts dictionary you can require that some
    constraints are satisfied multiple times. This example
    requires that constraint 7 is satisfied exactly twice:

        >>> sorted(next(MultiCover(data, counts={7: 2})))
        ['A', 'D', 'F']

    If only some constraints must be satisfied, these can be specified
    by the satisfy argument:

        >>> sorted(map(sorted, MultiCover(data, satisfy=[1, 2, 3, 4, 6, 7])))
        [['B', 'D', 'F'], ['B', 'E']]

    """
    # This implements Donald Knuth's "Algorithm X"
    # http://lanl.arxiv.org/pdf/cs/0011047.pdf

    def __init__(self, constraints, initial=(), random=False, counts=None,
                 satisfy=None):
        self.random = random
        self.constraints = constraints

        # A map from constraint to the set of choices that satisfy
        # that constraint.
        self.choices = defaultdict(set)
        for i in self.constraints:
            for j in self.constraints[i]:
                self.choices[j].add(i)

        # The multi-sets of constraints that must/may be satisfied
        # (and currently are not).
        self.must_satisfy = Counter()
        self.may_satisfy = Counter()
        if satisfy:
            self.must_satisfy.update(satisfy)
            self.may_satisfy.update(i for i in self.choices.keys()
                                    if i not in self.must_satisfy)
            if counts is not None:
                for k, v in counts.items():
                    if k in self.must_satisfy:
                        self.must_satisfy[k] = v
                    else:
                        self.may_satisfy[k] = v
        else:
            # All constraints must be satisfied.
            self.must_satisfy.update(self.choices.keys())
            if counts is not None:
                for k, v in counts.items():
                    self.must_satisfy[k] = v

        # Dictionary mapping constraint to the multiset it belongs to.
        self.multiset = dict()
        self.multiset.update((i, self.must_satisfy) for i in self.must_satisfy)
        self.multiset.update((i, self.may_satisfy) for i in self.may_satisfy)

        # The partial solution currently under consideration,
        # implemented as a stack of choices.
        self.solution = []

        # Make all the initial choices.
        try:
            for i in initial:
                self._choose(i)
            self.it = self._solve()
        except KeyError:
            # Initial choices were contradictory, so there are no solutions.
            self.it = iter(())

    def __next__(self):
        return next(self.it)

    next = __next__             # for compatibility with Python 2

    def _solve(self, best=None, last=None):
        if not self.must_satisfy:
            # No remaining constraints need to be satisfied.
            yield list(self.solution)
            return

        if best in self.must_satisfy:
            # Only consider choices that are "later" than the last choice (in
            # the arbitrary order induced by the id function). This avoids
            # transpostions by enforcing an order on the choices for each
            # constraint.
            choices = [i for i in self.choices[best] if id(i) > id(last)]
        else:
            # Find the constraint with the fewest remaining choices.
            best = min(self.must_satisfy, key=lambda j:len(self.choices[j]))
            choices = list(self.choices[best])
        if self.random:
            shuffle(choices)

        # Try each choice in turn and recurse.
        for i in choices:
            self._choose(i)
            for solution in self._solve(best, i):
                yield solution
            self._unchoose(i)

    def _choose(self, i):
        """Make choice i; mark constraints satisfied; and remove any
        choices that clash with it.

        """
        self.solution.append(i)
        for j in self.constraints[i]:
            multiset = self.multiset[j]
            multiset[j] -= 1
            if not multiset[j]:
                del multiset[j]
                for k in self.choices[j]:
                    for l in self.constraints[k]:
                        if l != j:
                            self.choices[l].remove(k)

    def _unchoose(self, i):
        """Unmake the most recent choice; restore constraints and choices."""
        assert i == self.solution.pop()
        for j in self.constraints[i]:
            multiset = self.multiset[j]            
            if j not in multiset:
                for k in self.choices[j]:
                    for l in self.constraints[k]:
                        if l != j:
                            self.choices[l].add(k)
            multiset[j] += 1

Using the multiset cover solver, I can represent the eight queens problem more straightforwardly, with no need for gadgets:

ConstraintPython representationMultiplicityRequired?
Queen on the \(i\)th row.The tuple ('row', i)1Yes
Queen on the \(j\)th column.The tuple ('column', j)1Yes
Queen on the \(k\)th leading diagonal.The tuple ('leading', k)1No
Queen on the \(k\)th trailing diagonal.The tuple ('trailing', k)1No
ChoicePython representation
Place a queen at \((i, j)\).The tuple ('queen', i, j)
from exactcover import MultiCover
from itertools import product

def queens(n=8, random=False):
    """Return an iterator that yields the solutions to the n-queens
    problem.

    Arguments:
    n -- The size of the problem. (Default: 8)
    random -- Generate solutions in random order? (Default: False.)

    """
    constraints = {
        ('queen', i, j): (('row', i), ('column', j),
                          ('leading', i - j + n - 1), ('trailing', i + j))
        for i, j in product(range(n), repeat=2)
    }
    return MultiCover(constraints, random=random,
                      satisfy=product(('row', 'column'), range(n)))

This is several times faster than the version of this function using gadgets. Here’s the version using gadgets finding all the 11×11 solutions:

>>> from timeit import timeit
>>> timeit(lambda:print(sum(1 for _ in queens(11))), number=1)
2680
8.097410938004032

Compare the above with the version using the multiset cover solver:

>>> timeit(lambda:print(sum(1 for _ in queens(11))), number=1)
2680
1.597617062041536

This can count the number of solutions for small board sizes (sequence A000170 in the OEIS):

>>> [sum(1 for _ in queens(i)) for i in range(4, 13)]
[2, 10, 4, 40, 92, 352, 724, 2680, 14200, 73712]

Now the sixteen queens solver is straightforward:

ConstraintPython representationMultiplicityRequired?
Queen on the \(i\)th row.The tuple ('row', i)2Yes
Queen on the \(j\)th column.The tuple ('column', j)2Yes
Queen on the \(k\)th leading diagonal.The tuple ('leading', k)2No
Queen on the \(k\)th trailing diagonal.The tuple ('trailing', k)2No
ChoicePython representation
Place a queen at \((i, j)\).The tuple ('queen', i, j)
from exactcover import MultiCover
from itertools import chain, product

def queens(n=8, random=False):
    """Return an iterator that yields the solutions to the 2n-queens
    problem.

    Arguments:
    n -- The size of the problem. (Default: 8)
    random -- Generate solutions in random order? (Default: False.)

    """
    constraints = {
        ('queen', i, j): (('row', i), ('column', j),
                          ('leading', i - j + n - 1), ('trailing', i + j))
        for i, j in product(range(n), repeat=2)
    }
    counts = dict.fromkeys(
        chain(product(('row', 'column'), range(n)),
              product(('leading', 'trailing'), range(2 * n + 1))), 2)
    return MultiCover(constraints, random=random, counts=counts,
                      satisfy=product(('row', 'column'), range(n)))

This quickly produces solutions:

>>> print(decode_queens(8, next(queens(8))))
♛++++♛++
♛+++++♛+
+++♛♛+++
+♛++++♛+
+++♛♛+++
+♛+++++♛
++♛++++♛
++♛++♛++

And it’s fast enough to count the number of solutions for small board sizes (sequence A225623 in the OEIS):

>>> [sum(1 for _ in queens(i)) for i in range(2, 9)]
[1, 2, 11, 92, 1097, 19448, 477136]

There are 92 solutions to the 5×5 instance of this problem, exactly as many as there were to the eight queens problem:

A. References