The story so far: I wrote a solver for a generic exact cover problem in Python, and used it to find polyomino tilings. The reason for using Python was because Python’s simplicity, terseness, and flexibility makes it easy to write solvers for other combinatorial problems. So let’s have a look at a second problem.
✴
In Sudoku you are given an \(n×n\) grid of squares, subdivided into \(m×{n\over m}\) blocks, and asked to fill each square with a digit between 1 and \(n\) inclusive such that every digit occurs exactly once in every row, column, and block. (Newspaper problems normally have \(n=9\) and \(m=3\).)
As in the case of polyomino tiling, I’d like to be able to draw a picture of a problem and see a picture of a solution, for example:
PROBLEM = ''' ...84...9 ..1.....5 8...2146. 7.8....9. ......... .5....3.1 .2491...7 9.....5.. 3...84... '''
So I need an encoder and decoder:
from itertools import chain, product from string import ascii_lowercase, ascii_uppercase DIGITS = '123456789' + ascii_uppercase + ascii_lowercase def encode_sudoku(problem=None): """Encode a Sudoku problem. Argument: problem  Given a Sudoku problem in the form of a string of n words, each word consisting of n characters, either a digit or a dot indicating a blank square. (Default: the blank problem.) Result: the pair (n, initial) where n is the size of the problem and initial is a list of tuples (i, j, d) meaning that the digit d is placed at (i, j) in the problem. >>> encode_sudoku('1... .... ...2 ....') (4, [(0, 0, '1'), (2, 3, '2')]) """ rows = problem.split() n = len(rows) if any(len(row) != n for row in rows): raise ValueError("all rows must have length {}".format(n)) initial = [(i, j, d) for i, row in enumerate(rows) for j, d in enumerate(row) if 0 <= DIGITS.find(d) < n] return n, initial def decode_sudoku(n, solution): """Format a Sudoku solution as a string. Arguments: n  The size of the problem. solution  The solution: an iterable of tuples (i, j, d) meaning that the digit d is placed at (i, j) in the solution. """ grid = [['.'] * n for _ in range(n)] for i, j, d in solution: grid[i][j] = d return '\n'.join(''.join(row) for row in grid)
(This choice of encoding comes in handy in make_sudoku
below.) Now the solver itself is simple. Here’s the encoding of the constraints and choices as hashable Python objects:
Constraint  Python representation 

The square at \((i, j)\) must contain a digit.  The tuple (i, j) .

The digit \(d\) must appear in row \(i\).  The tuple (d, 'row', i) .

The digit \(d\) must appear in column \(j\).  The tuple (d, 'col', j) .

The digit \(d\) must appear in block \((k, l)\).  The tuple (d, 'block', k, l) .

Choice  Python representation 
Place the digit \(d\) at \((i, j)\).  The tuple (i, j, d) .

And here’s the solver:
from exactcover import ExactCover from math import floor, sqrt def sudoku(n=9, problem=(), m=None, random=False): """Return an iterator that yields the solutions to a Sudoku problem. Arguments: n  The size of the problem, if no problem is given. (Default: 9.) problem  The problem to solve: an iterable of tuples (i, j, d) meaning that digit d appears at i, j. (Default: blank problem.) m  The height of the blocks. (Default: the square root of n, rounded down.) random  Generate solutions in random order? (Default: False.) """ if m is None: m = int(floor(sqrt(n))) if n <= 0 or len(DIGITS) < n or n % m: raise ValueError("bad dimensions: n={!r} m={!r}".format(n, m)) if not all(0 <= DIGITS.find(d) < n and 0 <= i < n and 0 <= j < n for i, j, d in problem): raise ValueError("bad problem: {!r}".format(problem)) constraints = { (i, j, d): ((i, j), (d, 'row', i), (d, 'col', j), (d, 'block', i // m, j // (n // m))) for i, j, d in product(range(n), range(n), DIGITS[:n]) } return ExactCover(constraints, problem, random)
Here’s a demo:
>>> n, problem = encode_sudoku(PROBLEM) >>> print(decode_sudoku(n, next(sudoku(n, problem)))) 632845179 471369285 895721463 748153692 163492758 259678341 524916837 986237514 317584926
The solver takes the height of the blocks (m
) as an optional argument in order be able to support problem sizes that aren’t perfect squares. For example, 12×12 problems typically have 3×4 blocks, but there is nothing that prevents problems of this size having 2×6 blocks, or even 1×12 blocks. When \(m=1\) the blocks are the same as the rows, and so add no additional contraint, meaning that the solutions are Latin squares. In the example below I use the Sudoku solver to (somewhat inefficiently) count the number of ‘reduced’ \(n×n\) Latin squares (that is, Latin squares where the first row and first column are the numbers from 1 to \(n\) in order):
>>> [sum(1 for _ in sudoku(*encode_sudoku( ... ' '.join([DIGITS[:n]] + [d + '.' * (n1) for d in DIGITS[:n]]))), m=1)) ... for n in range(1, 7)] ... [1, 1, 1, 4, 56, 9408]
This is sequence A000315 in the OEIS.
✴
This Sudoku solver is quick enough that I can use it to set problems. The idea is to create a random filled grid by solving a blank puzzle and passing random=True
to the ExactCover
constructor. Then erase randomly chosen squares (or pairs of squares if the aim is to make a symmetric puzzle) until no more squares can be erased without making the problem ambiguous.
from itertools import islice from random import sample class ImpossibleProblem(Exception): pass class AmbiguousProblem(Exception): pass def solve_sudoku(n=9, problem=(), m=None): """Solve the Sudoku problem and return its unique solution. If the problem is impossible, raise ImpossibleProblem. If the problem has multiple solutions, raise AmbiguousProblem. Optional arguments n, problem, and m are as for sudoku. """ solutions = list(islice(sudoku(n, problem, m), 2)) if len(solutions) == 1: return solutions[0] elif len(solutions) == 0: raise ImpossibleProblem('no solutions') else: raise AmbiguousProblem('two or more solutions') def make_sudoku(n=9, m=None): """Return a random nxn Sudoku problem with 180degree rotational symmetry. The problem returned is minimal in the sense that no symmetric pair of givens can be removed without making the problem ambiguous. The optional arguments are n, the size of the problem (default: 9) and m, the height of the blocks (default: the square root of n, rounded down). """ solution = sorted(next(sudoku(n=n, m=m, random=True))) given = [True] * n**2 k = (n**2 + 1) // 2 for i in sample(range(k), k): given[i] = given[i1] = False p = [solution[j] for j, g in enumerate(given) if g] try: solve_sudoku(n, p, m) problem = p except AmbiguousProblem: given[i] = given[i1] = True return problem
For example:
>>> print(decode_sudoku(9, make_sudoku())) 5..1.7... 6.....4.3 ..9....7. 7..6.5.3. 2.......4 .5.8.4..7 .4....5.. 9.7.....2 ...4.1..6