#!/usr/bin/env python # coding: utf-8 # # # # # # # # #
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 = ''''''.replace('3', str(B)) table = '' + D * ('' + D * '
{}') + '
' display(HTML(CSS + table.format(B, *grid))) # In[53]: show(puzzle0) # 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)) # In[9]: units[20] # In[10]: peers[20] # 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......')) # 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() # # 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 # 10010 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) # That works fine for a 4x4 array. But on a 9x9, with say, 64 empty squares, there will be 964 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](http://www.google.com/search?aq=f&sourceid=chrome&ie=UTF-8&q=10+GHz+*+1024+*+1+million+*+100+*+1+trillion+*+1+billion+*+13+billion+years+%2F+9**64+in+percent) # 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 10010. 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: #

#

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) # In[97]: do1(puzzle1) # In[98]: import cProfile cProfile.run("solve(puzzle1)") # # 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 # # 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) # In[87]: benchmark('3. AC') # # 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') # # 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') # # 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') #

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

    #
  1. 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. #
  2. 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. #
  3. 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. #
  4. 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. #