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 NPcomplete, 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 NPcomplete 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}
How might this be solved? Well, Knuth (2000) gives the following algorithm, which he calls ‘Algorithm X’:^{4}
Knuth comments on this:
Algorithm X is simply a statement of the obvious trialanderror approach. (Indeed, I can’t think of any other reasonable way to do the job, in general.)
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:
unsatisfied
.
choices
, mapping each constraint to the set of choices satisfying that constraint.
constraints
, mapping each choice to the set of constraints satisfied by that choice.
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]: self.choices[j].add(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. try: for i in initial: self._choose(i) 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) return # 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: shuffle(choices) # Try each choice in turn and recurse. for i in choices: self._choose(i) yield from self._solve() 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]: self.unsatisfied.remove(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 choice i; restore constraints and choices.""" last = self.solution.pop() assert i == last for j in self.constraints[i]: self.unsatisfied.add(j) for k in self.choices[j]: for l in self.constraints[k]: if l != j: self.choices[l].add(k)
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.
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:
TETROMINOES = ''' IIII LLL OO TTT ZZ L OO T ZZ ''' PENTOMINOES = ''' I 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 ASCIIart 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:
Constraint  Python 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' .

Choice  Python 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:
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 problem. 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)) OOTTTZZ OO TZZL IIIILLL """ 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[1rk] + oj if p not in region_set: break placing.append(p) else: 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)) YYYYWZZP LNYWWZPP LNWWZZPP LNN TTT LLN FTV UUXFFFTV UXXXFVVV UUXIIIII
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", stroke_width=3): """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", stroke_width=0.25) drawing.add(group) grid = [[_EMPTY] * w for _ in range(h)] for c, placing in solution: piece = drawing.g(fill=colour(c)) group.add(piece) 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) group.add(edges) for i, j in product(range(h + 1), range(w)): if ((_EMPTY if i == 0 else grid[i1][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][j1]) != (_EMPTY if j == w else grid[i][j])): edges.push(['M', j * size + oj, i * size + oi, 'l', 0, size]) drawing.save() 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], 'solution8x8.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) 520
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) 65 >>> solutions_svg(unique_solutions, 'all8x8.svg', columns=13, colour=COLOURS.get)
↩ 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.
↩ The colours are from this image by R. A. Nonenmacher, licensed under CCBYSA.
↩ 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.
↩ 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.
↩ 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.