Exact cover I


1. Introduction

Karp (1972) introduces the name ‘exact cover’ for the following decision problem:

INPUT: Family \(\{S_j\}\) of subsets of a set \(\{u_i, i=1,2,\ldots,t\}\)

PROPERTY: There is a subfamily \(\{T_h\} ⊆ \{S_j\}\) such that the sets \(T_h\) are disjoint and \(∪T_h = ∪S_j = \{u_i, i=1,2,\ldots,t\}\)

(The problem being to determine, for an arbitrary input, whether it possesses the property.) Karp goes on to prove that this problem is NP-complete, using a reduction from chromatic number, which has a reduction from satisfiability with at most three literals per clause (3SAT for short), which in turn has a reduction from satisfiability, which he proves is NP-complete directly.

I’m going to call the subsets \(S_j\) ‘choices’ and the \(u_i\) ‘constraints’, and to say that \(S_j\) ‘satisfies’ \(u_i\) iff \(u_i ∈ S_j\). The problem is then to determine if there is a set of choices that satisfies all the constraints exactly once each.1

Instances of this problem often arise in the solution of combinatorial problems, except that in this context we are interested in computing or counting the solutions, not just deciding if the problem has a solution. For example, consider the problem of tiling a region of 60 squares with the twelve pentominoes:2

In this problem the constraints are the 60 squares to be tiled (each square must be covered exactly once), and the 12 pentominoes to be placed (each must be placed exactly once). The choices are the different locations and orientations at which each pentomino can be placed.

Pentomino tiling problems were among the first combinatorial problems to be solved using a backtracking algorithm. Scott (1958) describes how he programmed the MANIAC II to find the 65 essentially different tilings of an 8×8 square (minus the central 2×2 square) by the twelve pentominoes, two of which are shown here.3

2. ‘Algorithm X’

How might this be solved? Well, Knuth (2000) gives the following algorithm, which he calls ‘Algorithm X’:4

  1. If the set of unsatisfied constraints is empty, the problem is solved. Yield the solution and return.
  2. Otherwise, choose an unsatisfied constraint c.
  3. If there are no choices that satisfy this constraint, the problem cannot be solved from this position. Return.
  4. Otherwise, for each choice i that satisfies c:
    1. For each constraint j satisfied by i:
      1. Remove j from the set of unsatisfied constraints.
      2. For each choice k that satisfies j, remove k from the set of choices.
    2. Recursively run this algorithm on the reduced set of choices and constraints.
    3. For each constraint j satisfied by i:
      1. Restore j to the set of unsatisfied constraints.
      2. For each choice k that satisfies j, restore k to the set of choices.

Knuth comments on this:

Algorithm X is simply a statement of the obvious trial-and-error approach. (Indeed, I can’t think of any other reasonable way to do the job, in general.)

3. Python implementation

An obvious question is, why Python, since combinatorial problems are computationally intensive, and Python is pretty slow? Well, Python’s combination of simplicity, terseness, and flexibility makes it a good language for experimenting with and modifying code. And Python’s native sets can be used to implement the sets of constraints and choices, meaning that I will be able to represent the constraints and choices themselves by any Python objects, avoiding the necessity for complex encoding and decoding. The algorithm needs:

from collections import defaultdict
from collections.abc import Iterator
from random import shuffle

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

    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.

    A solution is a list of choices such that each constraint is
    satisfied by exactly one of the choices in the solution.

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

    def __init__(self, constraints, initial=(), random=False):
        """Construct an ExactCover solver,

        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.)

        For example:

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

        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]:

        # The set of constraints which are currently unsatisfied.
        self.unsatisfied = set(self.choices)

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

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

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

    def _solve(self):
        if not self.unsatisfied:
            # No remaining unsatisfied constraints.
            yield list(self.solution)

        # Pick the constraint with the fewest remaining choices
        # (Knuth's "S heuristic").
        best = min(self.unsatisfied, key=lambda j:len(self.choices[j]))
        choices = list(self.choices[best])
        if self.random:

        # Try each choice in turn and recurse.
        for i in choices:
            yield from self._solve()

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

        for j in self.constraints[i]:
            for k in self.choices[j]:
                for l in self.constraints[k]:
                    if l != j:

    def _unchoose(self, i):
        """Unmake choice i; restore constraints and choices."""
        last = self.solution.pop()
        assert i == last
        for j in self.constraints[i]:
            for k in self.choices[j]:
                for l in self.constraints[k]:
                    if l != j:

This code is pretty compact, but it has a small flaw, namely the selection of the constraint with the fewest remaining choices (Knuth’s ‘S heuristic’). This needs to take time \(O(\log n)\) but here it takes \(O(n)\) and so will become a bottleneck as the number of constraints increases. However, it’s not easy to make a data structure that supports the necessary operations efficiently.5 I’ll try to revisit this in a later article in this series.

4. Pentomino tiling

First, input and output. It’s a good idea to arrange for these to be easily readable and checkable. I don’t want to have to type in lists of coordinates of tiles in pentominoes; I’d much rather draw a picture:

      L    OO   T    ZZ

I     L      N                      Y
I  FF L  PP  N TTT       V   W  X  YY ZZ
I FF  L  PP NN  T  U U   V  WW XXX  Y  Z
I  F  LL P  N   T  UUU VVV WW   X   Y  ZZ

If some other form is more convenient for processing, let the computer do the work:

from collections import defaultdict

_EMPTY = ' '

def encode_polyominoes(picture):
    """Given a "picture" of some polyominoes (in ASCII-art form), return a
    list of tuples (name, tiles) for each polyomino, where tiles are
    the coordinates (i, j) of its tiles.

    >>> sorted(encode_polyominoes(TETROMINOES))
    ... # doctest: +NORMALIZE_WHITESPACE
    [('I', ((1, 0), (1, 1), (1, 2), (1, 3))),
     ('L', ((1, 6), (1, 7), (1, 8), (2, 6))),
     ('O', ((1, 11), (1, 12), (2, 11), (2, 12))),
     ('T', ((1, 15), (1, 16), (1, 17), (2, 16))),
     ('Z', ((1, 20), (1, 21), (2, 21), (2, 22)))]

    pieces = defaultdict(list)
    for i, row in enumerate(picture.split('\n')):
        for j, c in enumerate(row):
            if c != _EMPTY:
                pieces[c].append((i, j))
    return [(c, tuple(tiles)) for c, tiles in pieces.items()]

def solution_tiles(solution):
    """Generate all the tiles in a polyomino solution."""
    return (t for _, tiles in solution for t in tiles)

def solution_bounds(solution):
    """Return the bounds of the region covered by a polyomino solution."""
    h, w = map(max, zip(*solution_tiles(solution)))
    return h + 1, w + 1

def decode_solution(solution):
    """Format a polyomino solution as a string."""
    h, w = solution_bounds(solution)
    grid = [[_EMPTY] * w for _ in range(h)]
    for c, tiles in solution:
        for i, j in tiles:
            grid[i][j] = c
    return '\n'.join(''.join(row) for row in grid)

Now for the encoding of a tiling problem as an exact cover instance:

ConstraintPython representation
The square with coordinates \((i, j)\) must be covered.The tuple (i, j).
The piece named ‘X’ must be placed somewhere.The string 'X'.
ChoicePython representation
Place the piece named ‘X’ so that it covers the squares \((i_0, j_0)\), \((i_1, j_1)\), …The tuple ('X', ((i0, j0), (i1, j1), ...)).

This encoding has the following useful properties:

  1. A solution (a list of choices) is completely self-descriptive.
  2. A solution is in the same format as the output of encode_polyominoes, which means I can pass the latter to any function taking a solution. This will come in handy later on when I have code for drawing solutions.
from exactcover import ExactCover
from itertools import product

def polyomino(pieces, region, random=False):
    """Return an iterator that yields the solutions to a polyomino tiling

    Required arguments:
    pieces -- an iterable of pieces, each being a tuple whose first
        element is the name of the piece and whose second element is a
        sequence of coordinates of the tiles in the piece (in an
        arbitrary coordinate system).
    region -- the region to fill (an iterable of coordinates x, y).

    Optional argument:
    random -- Generate solutions in random order? (Default: False.)

    >>> region = set(product(range(3), range(7))) - {(1, 2)}
    >>> sol = min(polyomino(encode_polyominoes(TETROMINOES), region))
    >>> print(decode_solution(sol))

    region = list(region)
    region_set = set(region)
    constraints = {}
    for c, tiles in pieces:
        # Origin of piece at the first tile.
        oi, oj = tiles[0]
        # Adjust tiles so that they are relative to the origin of the piece.
        tiles = [(ti - oi, tj - oj) for ti, tj in tiles]
        for (oi, oj), ri, rj, rk in product(region, (-1, +1), (-1, +1), (0, 1)):
            # Place origin of piece at oi, oj, reflecting vertically
            # if ri == -1, horizontally if rj == -1, and diagonally if
            # rk == 1.
            placing = []
            for t in tiles:
                p = ri * t[rk] + oi, rj * t[1-rk] + oj
                if p not in region_set:
                placing = tuple(sorted(placing))
                choice = c, placing
                if choice not in constraints:
                    constraints[choice] = (c,) + placing
    return ExactCover(constraints=constraints, random=random)

Now we’re all set to find tiling solutions:

>>> region = set(product(range(8), repeat=2)) - set(product((3, 4), repeat=2))
>>> solution = next(polyomino(encode_polyominoes(PENTOMINOES), region))
>>> print(decode_solution(solution))

5. Drawing solutions

Let’s draw that solution in Scalable Vector Graphics format, using the svgwrite Python package:

from svgwrite import Drawing

def solutions_svg(solutions, filename, columns=1, size=25, padding=5,
                  colour=lambda _: "white", stroke_colour="black",
    """Format polyomino tilings as an SVG image.

    Required arguments:
    solutions -- iterable of solutions to the tiling problem, each of
        which is a sequence of piece placements, each of which is a
        tuple whose first element is the name of the piece, and whose
        second element is a sequence of pairs (i, j) giving the
        locations of the tiles in the piece.
    filename -- where to save the SVG drawing.

    Optional arguments:
    columns -- number of solutions per row (default: 1).
    size -- width and height of each tile (default: 25).
    padding -- padding around the image (default: 5)
    colour -- function taking a piece name and returning its colour
        (default: a function returning white for each piece).
    stroke -- stroke colour (default: black).
    stroke_width -- width of strokes between pieces (default: 3).

    solutions = list(solutions)
    h, w = solution_bounds(solutions[0])
    rows = (len(solutions) + columns - 1) // columns
    drawing_size = (2 * padding + (columns * (w + 1) - 1) * size,
                    2 * padding + (rows    * (h + 1) - 1) * size)
    drawing = Drawing(debug=False, filename=filename, size=drawing_size)
    for i, solution in enumerate(solutions):
        y, x = divmod(i, columns)
        oj = padding + (x * (w + 1)) * size
        oi = padding + (y * (h + 1)) * size
        group = drawing.g(stroke=stroke_colour, stroke_linecap="round",
        grid = [[_EMPTY] * w for _ in range(h)]
        for c, placing in solution:
            piece = drawing.g(fill=colour(c))
            for i, j in placing:
                grid[i][j] = c
                piece.add(drawing.rect((j * size + oj, i * size + oi),
                                       (size, size)))
        edges = drawing.path(stroke_width=stroke_width)
        for i, j in product(range(h + 1), range(w)):
            if ((_EMPTY if i == 0 else grid[i-1][j])
                != (_EMPTY if i == h else grid[i][j])):
                edges.push(['M', j * size + oj, i * size + oi, 'l', size, 0])
        for i, j in product(range(h), range(w + 1)):
            if ((_EMPTY if j == 0 else grid[i][j-1])
                != (_EMPTY if j == w else grid[i][j])):
                edges.push(['M', j * size + oj, i * size + oi, 'l', 0, size])

COLOURS = dict(I="#EEAAAA", F="#DDBB99", L="#CCCC88",
               P="#BBDD99", N="#AAEEAA", T="#99DDBB",
               U="#88CCCC", V="#99BBDD", W="#AAAAEE",
               X="#BB99DD", Y="#CC88CC", Z="#DD99BB")

Then, using:

>>> solutions_svg([solution], 'solution-8x8.svg', colour=COLOURS.get)

produces the image:

How many solutions are there in all?

>>> all_solutions = list(polyomino(encode_polyominoes(PENTOMINOES), region))
>>> len(all_solutions)

But these include solutions that are reflections and rotations of each other. To reduce this collection to the essentially different solutions, one approach is to map each solution to a canonical form (by sorting the pieces and the tiles in each piece) and then discard duplicates:

from itertools import starmap

def canonical(solution):
    """Return a canonical version of a polyomino tiling solution."""
    h, w = solution_bounds(solution)
    all_tiles = sorted(solution_tiles(solution))
    orientations = [(False, True)] * 3
    if sorted((h - i - 1, j) for i, j in all_tiles) != all_tiles:
        # Can't reflect vertically.
        orientations[0] = (False,)
    if sorted((i, w - j - 1) for i, j in all_tiles) != all_tiles:
        # Can't reflect horizontally.
        orientations[1] = (False,)
    if sorted((j, i) for i, j in all_tiles) != all_tiles:
        # Can't reflect diagonally.
        orientations[2] = (False,)
    solution = sorted(solution)
    def oriented(ri, rj, rk):
        oriented_solution = []
        for c, tiles in solution:
            oriented_tiles = []
            for i, j in tiles:
                if ri: i = h - i - 1 # reflect vertically
                if rj: j = w - j - 1 # reflect horizontally
                if rk: i, j = j, i   # reflect diagonally
                oriented_tiles.append((i, j))
            oriented_solution.append((c, tuple(sorted(oriented_tiles))))
        return tuple(oriented_solution)
    return min(starmap(oriented, product(*orientations)))

This reduces the 520 solutions to the 65 essentially different solutions:

>>> unique_solutions = set(map(canonical, all_solutions))
>>> len(unique_solutions)
>>> solutions_svg(unique_solutions, 'all-8x8.svg', columns=13, colour=COLOURS.get)

A. References

  1.  I find these terms help me understand the problem and the algorithm. In Karp’s setting there are ‘sets’ and ‘subsets’, and in Knuth’s there are ‘rows’ and ‘columns’, but I find it hard to remember what these mean, especially in the context of a puzzle which may itself have rows and columns.

  2.  The colours are from this image by R. A. Nonenmacher, licensed under CC-BY-SA.

  3.  I haven’t been able to find a copy of this technical report online; all I know is from the summary in Knuth (2000). If you have a copy, or know where I can find one, please let me know.

  4.  Knuth presents this as a nondeterministic algorithm, but I think it’s slightly clearer to present it as a deterministic backtracking algorithm that yields all the solutions.

  5.  If you’re saying, “it’s easy, dummy, just drop in the SortedSet class from the sortedcontainers package,” it’s because you haven’t tried it yet.