Exact cover II


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:


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.

    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.

    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:

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

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

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 + '.' * (n-1) 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):

class AmbiguousProblem(Exception):

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')
        raise AmbiguousProblem('two or more solutions')

def make_sudoku(n=9, m=None):
    """Return a random nxn Sudoku problem with 180-degree rotational
    symmetry. The problem returned is minimal in the sense that no
    symmetric pair of givens can be removed without making the problem

    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[-i-1] = False
        p = [solution[j] for j, g in enumerate(given) if g]
            solve_sudoku(n, p, m)
            problem = p
        except AmbiguousProblem:
            given[i] = given[-i-1] = True
    return problem

For example:

>>> print(decode_sudoku(9, make_sudoku()))