Notebook

Peter Norvig, Nov. 2016

Solving Any Sudoku Puzzle, Quickly¶

Sudoku is a logic puzzle game. Given a grid like the one on below left, fill in the blank squares so that each column, row, and 3x3 box contains all the digits from 1 to 9. The solution is shown below right:

In this notebook I develop and implement an algorithm to solve any Sudoku puzzle, quickly. It turns out to be straightforward using two computer science ideas: constraint satisfaction and search. This notebook is an expansion of my previous Sudoku essay; that one emphasized simplicity; this one emphasizes efficiency and covers everything in more detail.

We will use these definitions:

A grid is a 9×9 array of squares. We can partition the grid

in three ways: - as 9 rows (of size 9×1) stacked on top of each other; - as 9 columns (of size 1×9) placed left-to-right; - as 9 boxes (of size 3×3) arranged in a 3×3 array.

A unit is a row,

column, or box: a collection of 9 squares.

A puzzle is a grid in which some squares have been

filled with a digit and some are empty.

A solution is a grid where every unit contains

all the digits from 1 to 9, and the filled squares of the original puzzle remain unchanged in the solution.

A well-formed puzzle has exactly one solution.
A peer of a square is any other

square in the same row, column, and box.

Sudoku data structures¶

Here are the choices I made to implement these concepts. I represent a grid as an 81-element list, with the upper left square being element 0, then going right-to-left and top-to-bottom with the lower-right square being element 80.

In [41]:

Grid   = list   # an 81-element list, left-to-right, top-to-bottom
Unit   = tuple  # a 9-element sequence
Square = int    # an int from 0 to 80; index into a Grid
Digit  = str    # a character from '1' to '9'
Digits = str    # One or more digits
empty  = '.'    # a square that has not been filled by a digit

puzzle0 = Grid('53..7....6..195....98....6.8...6...34..8.3..17...2...6.6....28....419..5....8..79')
answer0 = Grid('534678912672195348198342567859761423426853791713924856961537284287419635345286179')

def is_solution(solution: Grid, puzzle: Grid) -> bool: 
    """A solution to a puzzle is a grid where every unit contains all the digits 1 to 9
    and the filled squares of the original puzzle remain unchanged in the solution."""
    def unit_contents(unit): return {solution[s] for s in unit}
    return (all(unit_contents(unit) == digits for unit in all_units) and
            all(puzzle[s] == solution[s] for s in filled_squares(puzzle)))

def filled_squares(grid):
    "All the squares in a grid that are already filled with a digit."
    return (s for s in squares if grid[s] in digits)

def empty_squares(grid):
    "All the squares in a grid that are not already filled with a digit."
    return (s for s in squares if grid[s] not in digits)

Utility Functions and Displaying Grids¶

Here is the code to set up all the units and peers and related data structures, alog with the imports and utility functions we will need later on, and code for displaying grids in a pretty fashion. Although we will mostly use 81-square grids, the function use_grid_size(n) can be used to work with n-square grids, where n can be 1, 16, 81, or 256.

In [52]:

import itertools
import time
from IPython.display import HTML, display

def use_grid_size(n=81):
    global    N, D, B, digitstr, digits, squares, all_units, units, peers
    N         = n              # Number of squares in the grid
    D         = isqrt(N)       # Number of different digits
    B         = isqrt(D)       # A boxis of size B x B
    rows      = range(0, N, D) # Squares at left of each row: (0, 9, 18, ... 72)
    cols      = range(D)       # Squares at top of each column: (0, 1, 2, ... 8)
    digitstr  = '123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ_'[:D] # E.g., '123456789' for D=9
    digits    = set(digitstr)   # E.g., {1, 2, 3, 4, 5, 6, 7, 8, 9} for D=9
    squares   = range(N)    # Al1 the squares
    all_rows  = [cross(rows, [c]) for c in cols]
    all_cols  = [cross([r], cols) for r in rows]
    all_boxes = [cross(rs, cs) for rs in chunk(rows, B) for cs in chunk(cols, B)]
    all_units = all_rows + all_cols + all_boxes
    units     = [tuple(unit for unit in all_units if s in unit)
                 for s in squares]
    peers     = [tuple(union(units[s]) - {s})
                 for s in squares]

def cross(rows, cols) -> Unit:
    "A unit of all the ways we can add one of the row numbers to one of the column numbers."
    return Unit(r + c for r in rows for c in cols)

def union(collections) -> set: "Set union"; return set().union(*collections)

def isqrt(n) -> int: "Integer square root."; return int(n ** 0.5)

def chunk(sequence, B) -> list:
    "Chunk sequence into subsequences of length B."
    return [sequence[i:i+B] for i in range(0, len(sequence), B)]

def first(iterable) -> object: 
    "Return the first item from iterable, or None."
    return next(iterable, None)

def show(grid):
    "Display the Sudoku grid."
    use_grid_size(len(grid))
    CSS = '''<style>
    .s3 td {width:2em; height:2em; border:1px solid grey; text-align:center }
    .s3 td:first-child  {border-left:solid}
    .s3 td:nth-child(3n) {border-right:solid}
    .s3 tr:first-child  {border-top:solid}
    .s3 tr:nth-child(3n) {border-bottom:solid}
    </style>'''.replace('3', str(B))
    table = '<table class="s{}">' + D * ('<tr>' + D * '<td>{}') + '</table>'
    display(HTML(CSS + table.format(B, *grid)))

In [53]:

show(puzzle0)

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

We can how the numbering scheme for squares works; I'll display a grid that is filled, not with digits, but with square numbers. Then I'll show the units and peers for square 20:

In [44]:

show(range(81))

0	1	2	3	4	5	6	7	8
9	10	11	12	13	14	15	16	17
18	19	20	21	22	23	24	25	26
27	28	29	30	31	32	33	34	35
36	37	38	39	40	41	42	43	44
45	46	47	48	49	50	51	52	53
54	55	56	57	58	59	60	61	62
63	64	65	66	67	68	69	70	71
72	73	74	75	76	77	78	79	80

In [9]:

units[20]

Out[9]:

((2, 11, 20, 29, 38, 47, 56, 65, 74),
 (18, 19, 20, 21, 22, 23, 24, 25, 26),
 (0, 1, 2, 9, 10, 11, 18, 19, 20))

In [10]:

peers[20]

Out[10]:

(0, 65, 2, 1, 9, 74, 11, 10, 18, 19, 21, 22, 23, 24, 25, 26, 29, 38, 47, 56)

In [45]:

# A grid where 'X' marks square 20, 'r' marks squares in same row, 'c' marks squares in same column,
# and capital letters, 'C', 'R', and 'B', mark squares in the same box.
show(Grid('BBC......BBC......RRXrrrrrr..c........c........c........c........c........c......'))

B	B	C	.	.	.	.	.	.
B	B	C	.	.	.	.	.	.
R	R	X	r	r	r	r	r	r
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.

We need a test suite:

In [46]:

puzzle0 = Grid('53..7....6..195....98....6.8...6...34..8.3..17...2...6.6....28....419..5....8..79')
answer0 = Grid('534678912672195348198342567859761423426853791713924856961537284287419635345286179')
wrong0  = Grid('532678912672195348198342567859761423426853791713924856961537284287419635345286179')
puzzle1 = Grid('4.7.6.8.5.3....9.....7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......')
answer1 = Grid('417369825632158947958724316825437169791586432346912758289643571573291684164875293')
wrong1  = Grid('427369825632158947958724316825437169791586432346912758289643571573291684164875293')

def tests():
    "A suite of unit tests."
    use_grid_size(81)
    assert N == 81 and D == 9 and B == 3
    assert len(squares) == 81
    assert len(all_units) == 27
    for s in squares:
        assert len(units[s]) == 3 
        assert len(peers[s]) == 20
    assert union(['feed', 'beef', 'cafe', 'face']) == {'a', 'b', 'c', 'd', 'e', 'f'}
    assert isqrt(16) == 4 and isqrt(81) == 9 and isqrt(256) == 16
    assert cross((0, 9, 18), (0, 1, 2)) == (0, 1, 2, 9, 10, 11, 18, 19, 20)
    assert chunk('abcdefghi', 3) == ['abc', 'def', 'ghi']
    assert units[20] == ((2, 11, 20, 29, 38, 47, 56, 65, 74),
                         (18, 19, 20, 21, 22, 23, 24, 25, 26),
                         (0, 1, 2, 9, 10, 11, 18, 19, 20))
    assert set(peers[20]) == {0, 1, 2, 9, 10, 11, 18, 19, 21, 22, 
                              23, 24, 25, 26, 29, 38, 47, 56, 65, 74}
    assert is_solution(answer0, puzzle0)
    assert not is_solution(wrong0, puzzle0)
    assert is_solution(answer1, puzzle1)
    assert not is_solution(wrong1, puzzle1)
    return 'pass'

tests()

Out[46]:

'pass'

Solver 1: Generate and Test¶

Here is a simple (but very inefficient) algorithm to solve Sudoku puzzles:

Generate and test solver:

First fill in all the squares with digits. Then check to see if we have a solution. If not, then fill the squares with a different combination of digits. Repeat with different combinations until a solution is found.

Sometimes generate and test is the best you can do. Herb Simon, in his classic book [The Sciences of the Artificial](http://books.google.com/books?id=k5Sr0nFw7psC) (page 194), describes the task of opening a safe that has ten dials, each with 100 numbers on it. If the safe is well-designed there is no better approach than to generate-and-test with all 100¹⁰ combinations. The generate-and-test solver is simple:

In [49]:

def solve(puzzle): 
    "Find a solution to the puzzle."
    return first(grid for grid in generate_all_grids(puzzle) 
                 if is_solution(grid, puzzle))

def generate_all_grids(grid):
    "An iterable of all possible ways to fill in grid."
    values = [(digitstr if d is empty else d) for d in grid]
    return itertools.product(*values)

In [54]:

def do1(puzzle):
    "Do one puzzle; showing puzzle and solution and printing elapsed time."
    show(puzzle)
    t0 = time.clock()
    solution = solve(Grid(puzzle))
    t1 = time.clock()
    assert is_solution(solution, puzzle)
    show(solution)
    print('{:.3f} seconds'.format(t1 - t0))

use_grid_size(16)
tiny = Grid('...44.3..4...241')
do1(tiny)

.	.	.	4
4	.	3	.
.	4	.	.
.	2	4	1

2	3	1	4
4	1	3	2
1	4	2	3
3	2	4	1

0.210 seconds

That works fine for a 4x4 array. But on a 9x9, with say, 64 empty squares, there will be 9⁶⁴ combinations of digits. Suppose we have access to a secret new custom 10 GHz GPU capable of doing 1024 generate and test operations in a single clock cycle, and let's say we can fit a million of these units in a data center, and we could afford a hundred data centers, and while we're shopping, let's say we also pick up a time machine and go back 13 billion years to the early days of the universe, and a fleet of starships with which we visit a trillion galaxies and in each galaxy convince the inhabitants of a billion planets to each buy a similar setup, and we all start our programs running, managing to distribute the cases perfectly so that no computation is duplicated or wasted. Then we can estimate that we'd be only 3% of the way through with the first puzzle.

Generate and test is not the algorithm you're looking for. Move along.

Solver 2: Combinatorial Search¶

In Simon's safe example, he next supposes a defective safe, in which a faint click is heard whenever a dial is set to the right position. With this safe it only takes 100 × 10 trials to discover the correct combination, not 100¹⁰. The moral is that if we have some selective knowledge of components of the safe (or puzzle) that tell us if a partial solution is on track, our search will be much easier.

Fortunately for us, Sudoku is like the defective safe. To be precise, it is like a safe with 81 dials, each with 9 numbers, and sometimes if you put one of the dials in the wrong position you hear a sound, but sometimes not (depending on the positions of the other dials and on how carefully you listen). For example, in the lower-left corner of the tiny puzzle above, if we try 1, 2, or 4 we would immediately get feedback that those digits won't work, because they already appear elsewhere in the bottom row. Therefore, the lower-left corner must be filled by a 3. There are 262,144 ways to fill in the 9 empty squares, but right away we can eliminate 196,608 of them; we need only consider the ones that have a 3 in that position.

Unfortunately for us, Sudoku does not give up all its secrets so easily. Sometimes, when we consider a square, we can't immediately tell what digit to put there. For example, in the upper-left of the tiny puzzle, only 4 can be eliminated as a possibility. That's ok--we're allowing ourselves the use of an eraser, so just guess one of the remaining three possibilities; if it doesn't pan out, erase it and try one of the others. This gives us:

Combinatorial search algorithm: Start filling squares with digits, one at a time, always making sure that the digit we put in a square is not the same as any peer. If there are several possible digits, pick one, but remember the others. If we reach a square where every digit conflicts with a peer, back up and consider a different digit for the previous square. Stop when we successfully fill every square, or give up if we tried all possibilities without finding a solution.

The progress of this algorithm can be described as a search tree, where each node represents a partially-filled (incremental) state of the puzzle, and a branch corresponds to filling in a square with a digit. (We display the new digits in red.) Here is a search tree for the tiny puzzle:

A few things to note about this particular search tree:

Each successive level fills in one square, so with nine empty squares, the solution must appear on the ninth level of the tree (counting the root as the zeroth level).
This is a well-formed puzzle, so there will be exactly one solution, and thus exactly one node at the ninth level.
A rectangle with a black X indicates an inconsistency: there is no way to continue because the next unfilled square already has all four digits taken by its peers.
When we hit a X we have to pull out our eraser, or in other words we have to backtrack--go up a level in the tree and try the next possible digit for the previous square. On the far-left branch of the tree we have to back up three levels---erasing first the 2, then the 3, then the 1--before we can pencil in another possible digit (the 2 in the upper-left corner which we see in the middle branch of the tree).
The rightmost branch of the tree remains unexplored--the search proceeds left-to-right and has not got to that brnach yet. (But since the puzzle is well-formed and we have found the solution, we know the unexplored branches cannot contain another solution.)

The code to implement combinatorial search is straightforward:

In [99]:

def search(grid) -> Grid:
    """Select an unfilled square, try each possible digit there, recursively searching.
    When all squares filled: success; when no more digits to try: return None for failure."""
    square = select_empty_square(grid)
    if square is None: # No empty square means the grid is a solution
        return grid
    for digit in possible_digits(grid, square):
        result = search(assign(Grid(grid), square, digit))
        if result: 
            return result

solve = search

def select_empty_square(grid) -> Square:
    "Return the first square that is not filled with a digit; or None if all squares are filled."
    return (None if grid is None else first(empty_squares(grid)))

def assign(grid, s, d) -> Grid:
    "For now, simply assign grid[s] = d."
    grid[s] = d
    return grid

def possible_digits(grid, s):
    "A square can be filled with any digit that is not already taken by a peer."
    peer_digits = {grid[s] for s in peers[s]}
    return digits - peer_digits

The key function is search. It obeys the convention that if it is passed Inconsistent (to indicate an invalid grid) it returns Inconsistent, meaning that no solution can be found. Otherwise it selects some unfilled square to work on. If select_square returns None that is actually good news: it means that all the squares are already filled, and we are done: we have a solution, so we just return it. Otherwise, for each possible digit that could fill square s we try to assign that digit to s and search for a solution from there. If some call to search succeeds, return it; but if a digit assignment causes the search to fail, just keep going with the next digit assignment. If all digit assignments fail, then the whole call to search fails, and we back up to the previous recursive call.

We call this type of combinatorial search a constraint satisfaction search: think of each square as being a variable that can take on a value (a digit), and the rules of Sudoku as being constraints on the possible values. A state in the search tree is then an assignment of values to some subset of variables in a way that does not violate any constraint. The constraint satisfaction approach has found many fun applications in puzzles and serious applications in problem solving.

Here we see solve (and therefore search) in action:

That's all there is to it! The only reason we don't stop this article now is that the program is still too slow for some purposes. We're not talking billions-of-years slow, but it did take 3 minutes and 45 seconds to solve this puzzle (a rather hard one). On a benchmark of 50 easy puzzles the program was fast enough, taking a total of 9.5 seconds (a rate of 5 puzzles per second, or 5 Hz).

In [96]:

do1(puzzle0)

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

5	3	4	6	7	8	9	1	2
6	7	2	1	9	5	3	4	8
1	9	8	3	4	2	5	6	7
8	5	9	7	6	1	4	2	3
4	2	6	8	5	3	7	9	1
7	1	3	9	2	4	8	5	6
9	6	1	5	3	7	2	8	4
2	8	7	4	1	9	6	3	5
3	4	5	2	8	6	1	7	9

0.024 seconds

In [97]:

do1(puzzle1)

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

4	1	7	3	6	9	8	2	5
6	3	2	1	5	8	9	4	7
9	5	8	7	2	4	3	1	6
8	2	5	4	3	7	1	6	9
7	9	1	5	8	6	4	3	2
3	4	6	9	1	2	7	5	8
2	8	9	6	4	3	5	7	1
5	7	3	2	9	1	6	8	4
1	6	4	8	7	5	2	9	3

17.216 seconds

In [98]:

import cProfile
cProfile.run("solve(puzzle1)")

         63904089 function calls (62617098 primitive calls) in 30.046 seconds

   Ordered by: standard name

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
 52321162    8.946    0.000    8.946    0.000 <ipython-input-1-9492bb1f8314>:17(is_filled)
  1286992    0.444    0.000   20.602    0.000 <ipython-input-52-b1db65fde4cb>:36(first)
1286992/1    3.380    0.000   30.046   30.046 <ipython-input-95-3dd63506ce25>:1(search)
  1286992    1.496    0.000   22.374    0.000 <ipython-input-95-3dd63506ce25>:14(select_empty_square)
  2573983   11.162    0.000   20.108    0.000 <ipython-input-95-3dd63506ce25>:17(<genexpr>)
  1286991    0.238    0.000    0.238    0.000 <ipython-input-95-3dd63506ce25>:19(assign)
  1286991    1.272    0.000    4.055    0.000 <ipython-input-95-3dd63506ce25>:24(possible_digits)
  1286991    2.783    0.000    2.783    0.000 <ipython-input-95-3dd63506ce25>:26(<setcomp>)
        1    0.000    0.000   30.046   30.046 <string>:1(<module>)
        1    0.000    0.000   30.046   30.046 {built-in method builtins.exec}
  1286992    0.326    0.000   20.159    0.000 {built-in method builtins.next}
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}

Solver 3. Arc Consistency¶

In [93]:

def solve(puzzle: Grid) -> Grid:
    grid = Grid(digitstr if d is empty else d
                for d in puzzle)
    return search(grid)

def possible_digits(grid, s) -> Digits: return grid[s]

def assign(grid, s, d) -> Grid:
    """Assign grid[s] = d and eliminate d from the peers of s.
    Return the updated grid, or Inconsistent if inconsistency detected."""
    if d not in grid[s]: return Inconsistent # d is not among the possibilities
    grid[s] = d
    if not all(eliminate(grid, p, d) for p in peers[s]): 
        return None
    return grid

def eliminate(grid, s, d) -> bool:
    "Remove d from possibilities for grid[s]. If checking finds an inconsistency return None."
    if d not in grid[s]: 
        return True # Already eliminated d
    grid[s] = grid[s].replace(d, '')
    return all(constraint(grid, s, d) 
               for constraint in constraints)

def arc_consistent(grid, s, d):
    "Return true if s has multiple digits left, or one that we can consistently assign."
    ndigits = len(grid[s])
    return ndigits >= 2 or (ndigits == 1 and assign(grid, s, grid[s]))

constraints = [arc_consistent]

In [94]:

do1(puzzle1) # Should be 0.2

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-94-ed46d3590e0c> in <module>()
----> 1 do1(puzzle1) # Should be 0.2

<ipython-input-54-65f8133e36ab> in do1(puzzle)
      3     show(puzzle)
      4     t0 = time.clock()
----> 5     solution = solve(Grid(puzzle))
      6     t1 = time.clock()
      7     assert is_solution(solution, puzzle)

<ipython-input-93-3fd012a56996> in solve(puzzle)
      2     grid = Grid(digitstr if d is empty else d
      3                 for d in puzzle)
----> 4     return search(grid)
      5 
      6 def possible_digits(grid, s) -> Digits: return grid[s]

<ipython-input-56-3dd63506ce25> in search(grid)
      6         return grid
      7     for digit in possible_digits(grid, square):
----> 8         result = search(assign(Grid(grid), square, digit))
      9         if result:
     10             return result

<ipython-input-56-3dd63506ce25> in search(grid)
      6         return grid
      7     for digit in possible_digits(grid, square):
----> 8         result = search(assign(Grid(grid), square, digit))
      9         if result:
     10             return result

<ipython-input-56-3dd63506ce25> in search(grid)
      2     """Select an unfilled square, try each possible digit there, recursively searching.
      3     When all squares filled: success; when no more digits to try: return None for failure."""
----> 4     square = select_empty_square(grid)
      5     if square is None: # No empty square means the grid is a solution
      6         return grid

<ipython-input-91-5e0e53165fa8> in select_empty_square(grid)
      1 def select_empty_square(grid):
      2     "Return an unfilled square with the fewest possible digits; or None if no unfilled squares."
----> 3     unfilled = [s for s in squares if len(grid[s]) > 1]
      4     return min(unfilled, key=lambda s: len(grid[s])) if unfilled else None

<ipython-input-91-5e0e53165fa8> in <listcomp>(.0)
      1 def select_empty_square(grid):
      2     "Return an unfilled square with the fewest possible digits; or None if no unfilled squares."
----> 3     unfilled = [s for s in squares if len(grid[s]) > 1]
      4     return min(unfilled, key=lambda s: len(grid[s])) if unfilled else None

TypeError: 'NoneType' object is not subscriptable

Benchmarking¶

In [86]:

puzzles = """
4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......
52...6.........7.13...........4..8..6......5...........418.........3..2...87.....
6.....8.3.4.7.................5.4.7.3..2.....1.6.......2.....5.....8.6......1....
48.3............71.2.......7.5....6....2..8.............1.76...3.....4......5....
....14....3....2...7..........9...3.6.1.............8.2.....1.4....5.6.....7.8...
......52..8.4......3...9...5.1...6..2..7........3.....6...1..........7.4.......3.
6.2.5.........3.4..........43...8....1....2........7..5..27...........81...6.....
.524.........7.1..............8.2...3.....6...9.5.....1.6.3...........897........
6.2.5.........4.3..........43...8....1....2........7..5..27...........81...6.....
.923.........8.1...........1.7.4...........658.........6.5.2...4.....7.....9.....
6..3.2....5.....1..........7.26............543.........8.15........4.2........7..
.6.5.1.9.1...9..539....7....4.8...7.......5.8.817.5.3.....5.2............76..8...
..5...987.4..5...1..7......2...48....9.1.....6..2.....3..6..2.......9.7.......5..
3.6.7...........518.........1.4.5...7.....6.....2......2.....4.....8.3.....5.....
1.....3.8.7.4..............2.3.1...........958.........5.6...7.....8.2...4.......
6..3.2....4.....1..........7.26............543.........8.15........4.2........7..
....3..9....2....1.5.9..............1.2.8.4.6.8.5...2..75......4.1..6..3.....4.6.
45.....3....8.1....9...........5..9.2..7.....8.........1..4..........7.2...6..8..
.237....68...6.59.9.....7......4.97.3.7.96..2.........5..47.........2....8.......
..84...3....3.....9....157479...8........7..514.....2...9.6...2.5....4......9..56
""".split()[5:6]

benchmarks = {}

def benchmark(label, puzzles=puzzles):
    "Run `solve` on these puzzles; record and verify results for this label; print all results."
    n = len(puzzles)
    t0 = time.clock()
    results = [solve(Grid(p)) for p in puzzles]
    avg = (time.clock() - t0) / len(puzzles)
    for (r, p) in zip(results, puzzles):
        assert is_solution(r, p) 
    benchmarks[label] = '{:.3f} sec/puzzle ({:.1f} Hz)'.format(avg, 1/avg)
    for L in sorted(benchmarks):
        print('{:10s} {}'.format(L, benchmarks[L]))

In [84]:

len(puzzles)

Out[84]:

In [87]:

benchmark('3. AC')

3. AC      11.219 sec/puzzle (0.1 Hz)

4. Dual Consistency¶

In [72]:

def dual_consistent(grid, s, d):
    """After eliminating d from grid[s], check each unit of s and make sure there is some
    position in the unit for d. If only one possible place left for d, assign it."""
    for u in units[s]:
        places_for_d = [s2 for s2 in u if d in grid[s2]]
        nplaces = len(places_for_d)
        if nplaces==0 or (nplaces==1 and not assign(grid, places_for_d[0], d)):
            return None
    return True

In [88]:

constraints = [arc_consistent, dual_consistent]
benchmark('4. AC+Dual')

3. AC      11.219 sec/puzzle (0.1 Hz)
4. AC+Dual 26.392 sec/puzzle (0.0 Hz)

5. Naked pairs¶

In [89]:

def naked_pairs(grid, s, ignore):
    """Look for two squares in a unit with the same two possible digits. 
    For example, if s and s2 both have the value '35', then we know that 3 and 5
    must go in those two squares. We don't know which is which, but we can eliminate 
    3 and 5 from any other square s3 that is in the unit."""
    vals = grid[s]
    if len(vals) != 2: return True
    for u in units[s]:
        for s2 in u:
            if s2 != s and grid[s2] == vals:
                # Found naked pair: s and s2; remove their two vals from others in unit
                for s3 in u:
                    if s != s3 != s2:
                        if not all(eliminate(grid, s3, v) for v in vals):
                            return None
    return True

In [90]:

constraints = [arc_consistent, dual_consistent, naked_pairs]
benchmark('5. AC+D+NP')

3. AC      11.219 sec/puzzle (0.1 Hz)
4. AC+Dual 26.392 sec/puzzle (0.0 Hz)
5. AC+D+NP 29.821 sec/puzzle (0.0 Hz)

6. Variable Ordering¶

In [91]:

def select_empty_square(grid):
    "Return an unfilled square with the fewest possible digits; or None if no unfilled squares."
    unfilled = [s for s in squares if len(grid[s]) > 1]
    return min(unfilled, key=lambda s: len(grid[s])) if unfilled else None

In [92]:

benchmark('6. VarOrd')

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-92-3052b7d99217> in <module>()
----> 1 benchmark('6. VarOrd')

<ipython-input-86-a0d9359377e8> in benchmark(label, puzzles)
     28     n = len(puzzles)
     29     t0 = time.clock()
---> 30     results = [solve(Grid(p)) for p in puzzles]
     31     avg = (time.clock() - t0) / len(puzzles)
     32     for (r, p) in zip(results, puzzles):

<ipython-input-86-a0d9359377e8> in <listcomp>(.0)
     28     n = len(puzzles)
     29     t0 = time.clock()
---> 30     results = [solve(Grid(p)) for p in puzzles]
     31     avg = (time.clock() - t0) / len(puzzles)
     32     for (r, p) in zip(results, puzzles):

<ipython-input-85-3fd012a56996> in solve(puzzle)
      2     grid = Grid(digitstr if d is empty else d
      3                 for d in puzzle)
----> 4     return search(grid)
      5 
      6 def possible_digits(grid, s) -> Digits: return grid[s]

<ipython-input-56-3dd63506ce25> in search(grid)
      6         return grid
      7     for digit in possible_digits(grid, square):
----> 8         result = search(assign(Grid(grid), square, digit))
      9         if result:
     10             return result

<ipython-input-56-3dd63506ce25> in search(grid)
      6         return grid
      7     for digit in possible_digits(grid, square):
----> 8         result = search(assign(Grid(grid), square, digit))
      9         if result:
     10             return result

<ipython-input-56-3dd63506ce25> in search(grid)
      2     """Select an unfilled square, try each possible digit there, recursively searching.
      3     When all squares filled: success; when no more digits to try: return None for failure."""
----> 4     square = select_empty_square(grid)
      5     if square is None: # No empty square means the grid is a solution
      6         return grid

<ipython-input-91-5e0e53165fa8> in select_empty_square(grid)
      1 def select_empty_square(grid):
      2     "Return an unfilled square with the fewest possible digits; or None if no unfilled squares."
----> 3     unfilled = [s for s in squares if len(grid[s]) > 1]
      4     return min(unfilled, key=lambda s: len(grid[s])) if unfilled else None

<ipython-input-91-5e0e53165fa8> in <listcomp>(.0)
      1 def select_empty_square(grid):
      2     "Return an unfilled square with the fewest possible digits; or None if no unfilled squares."
----> 3     unfilled = [s for s in squares if len(grid[s]) > 1]
      4     return min(unfilled, key=lambda s: len(grid[s])) if unfilled else None

TypeError: 'NoneType' object is not subscriptable

However, if you want to solve a million hard puzzles, you probably don't want to wait a few years, and you would prefer a faster algorithm. Look again at the search tree above, and consider how we can find the solution faster. Here are four general strategies:

We could try to move the solution to the left. Since we go top-down, left-to-right, this means we would find it sooner. This can potentially be done with better value ordering, the order in which the function possible_values considers which digit to try next. If we were able to choose right every time, there would in effect be no search--the solution would be in the lower-left corner of the tree and we would arrive at it directly with no backtracking.
We could try to move the solution closer to the top. If we only fill in one square at a time then the solution will always be at the same level (the number of empty squares in the original puzzle). But if we could fill in several squares at a time, the solution gets closer to the top and we can find it faster. Filling in multiple squares at a time can be accomplished with constraint propagation--using the constraints imposed by the rules of Sudoku to conclude that an unfilled square must take on a certain digit, thus sanctioning us to fill in several squares at once.
We could prune branches (mark them with an X and backtrack) earlier. This can also be done with constraint propagation, but in this case rather than filling in squares, we would notice that some square has no possible digit, and thus that we can backtrack immediately.
We could consider the squares in a better order, a process called variable ordering, which in our program is implemented by the function select_square. At each node in the tree we pick one square to branch on. The current implementation of select_square just considers the squares in order from 0 to 80 and chooses the first unfilled one, but it could be more clever. A different ordering could be better (or worse) for any of the three reasons above.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

0	1	2	3	4	5	6	7	8
9	10	11	12	13	14	15	16	17
18	19	20	21	22	23	24	25	26
27	28	29	30	31	32	33	34	35
36	37	38	39	40	41	42	43	44
45	46	47	48	49	50	51	52	53
54	55	56	57	58	59	60	61	62
63	64	65	66	67	68	69	70	71
72	73	74	75	76	77	78	79	80

B	B	C	.	.	.	.	.	.
B	B	C	.	.	.	.	.	.
R	R	X	r	r	r	r	r	r
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

5	3	4	6	7	8	9	1	2
6	7	2	1	9	5	3	4	8
1	9	8	3	4	2	5	6	7
8	5	9	7	6	1	4	2	3
4	2	6	8	5	3	7	9	1
7	1	3	9	2	4	8	5	6
9	6	1	5	3	7	2	8	4
2	8	7	4	1	9	6	3	5
3	4	5	2	8	6	1	7	9

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

4	1	7	3	6	9	8	2	5
6	3	2	1	5	8	9	4	7
9	5	8	7	2	4	3	1	6
8	2	5	4	3	7	1	6	9
7	9	1	5	8	6	4	3	2
3	4	6	9	1	2	7	5	8
2	8	9	6	4	3	5	7	1
5	7	3	2	9	1	6	8	4
1	6	4	8	7	5	2	9	3

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

0	1	2	3	4	5	6	7	8
9	10	11	12	13	14	15	16	17
18	19	20	21	22	23	24	25	26
27	28	29	30	31	32	33	34	35
36	37	38	39	40	41	42	43	44
45	46	47	48	49	50	51	52	53
54	55	56	57	58	59	60	61	62
63	64	65	66	67	68	69	70	71
72	73	74	75	76	77	78	79	80

B	B	C	.	.	.	.	.	.
B	B	C	.	.	.	.	.	.
R	R	X	r	r	r	r	r	r
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

5	3	4	6	7	8	9	1	2
6	7	2	1	9	5	3	4	8
1	9	8	3	4	2	5	6	7
8	5	9	7	6	1	4	2	3
4	2	6	8	5	3	7	9	1
7	1	3	9	2	4	8	5	6
9	6	1	5	3	7	2	8	4
2	8	7	4	1	9	6	3	5
3	4	5	2	8	6	1	7	9

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

4	1	7	3	6	9	8	2	5
6	3	2	1	5	8	9	4	7
9	5	8	7	2	4	3	1	6
8	2	5	4	3	7	1	6	9
7	9	1	5	8	6	4	3	2
3	4	6	9	1	2	7	5	8
2	8	9	6	4	3	5	7	1
5	7	3	2	9	1	6	8	4
1	6	4	8	7	5	2	9	3

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

0	1	2	3	4	5	6	7	8
9	10	11	12	13	14	15	16	17
18	19	20	21	22	23	24	25	26
27	28	29	30	31	32	33	34	35
36	37	38	39	40	41	42	43	44
45	46	47	48	49	50	51	52	53
54	55	56	57	58	59	60	61	62
63	64	65	66	67	68	69	70	71
72	73	74	75	76	77	78	79	80

B	B	C	.	.	.	.	.	.
B	B	C	.	.	.	.	.	.
R	R	X	r	r	r	r	r	r
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.
.	.	c	.	.	.	.	.	.

5	3	.	.	7	.	.	.	.
6	.	.	1	9	5	.	.	.
.	9	8	.	.	.	.	6	.
8	.	.	.	6	.	.	.	3
4	.	.	8	.	3	.	.	1
7	.	.	.	2	.	.	.	6
.	6	.	.	.	.	2	8	.
.	.	.	4	1	9	.	.	5
.	.	.	.	8	.	.	7	9

5	3	4	6	7	8	9	1	2
6	7	2	1	9	5	3	4	8
1	9	8	3	4	2	5	6	7
8	5	9	7	6	1	4	2	3
4	2	6	8	5	3	7	9	1
7	1	3	9	2	4	8	5	6
9	6	1	5	3	7	2	8	4
2	8	7	4	1	9	6	3	5
3	4	5	2	8	6	1	7	9

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.

4	1	7	3	6	9	8	2	5
6	3	2	1	5	8	9	4	7
9	5	8	7	2	4	3	1	6
8	2	5	4	3	7	1	6	9
7	9	1	5	8	6	4	3	2
3	4	6	9	1	2	7	5	8
2	8	9	6	4	3	5	7	1
5	7	3	2	9	1	6	8	4
1	6	4	8	7	5	2	9	3

4	.	7	.	6	.	8	.	5
.	3	.	.	.	.	9	.	.
.	.	.	7	.	.	.	.	.
.	2	.	.	.	.	.	6	.
.	.	.	.	8	.	4	.	.
.	.	.	.	1	.	.	.	.
.	.	.	6	.	3	.	7	.
5	.	.	2	.	.	.	.	.
1	.	4	.	.	.	.	.	.