Solving Every Sudoku Puzzle Quickly

def parse(gridstring):
    """Convert a string into a grid: a list of values. Accepts digits (or '.' or '0'
    for empty); all other characters ignored. In result, '.' used for empty."""
    gridstring = gridstring.replace('0', '.')
    grid = [c for c in gridstring if c in digits or c == '.']
    assert len(grid) == N*N
    return grid

def show(grid):
    "Convert a grid to a gridstring laid out in 2D, with lines between boxes."
    return "INCONSISTENCY" if grid is False else show_template.format(*grid)

show_line = '\n' + '+'.join(['--'*Nbox]*Nbox) + '\n'
show_template = show_line.join(['\n'.join(['|'.join(['{} '*Nbox]*Nbox)]*Nbox)]*Nbox)

>>> puzzle1 = "4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......"
>>> parse(puzzle1)
['4', '.', '.', '.', '.', '.', '8', '.', '5',
 '.', '3', '.', '.', '.', '.', '.', '.', '.',
 '.', '.', '.', '7', '.', '.', '.', '.', '.',
 '.', '2', '.', '.', '.', '.', '.', '6', '.',
 '.', '.', '.', '.', '8', '.', '4', '.', '.',
 '.', '.', '.', '.', '1', '.', '.', '.', '.',
 '.', '.', '.', '6', '.', '3', '.', '7', '.',
 '5', '.', '.', '2', '.', '.', '.', '.', '.',
 '1', '.', '4', '.', '.', '.', '.', '.', '.']

>>> print show(parse(puzzle1))
4 . . |. . . |8 . 5
. 3 . |. . . |. . .
. . . |7 . . |. . .
------+------+------
. 2 . |. . . |. 6 .
. . . |. 8 . |4 . .
. . . |. 1 . |. . .
------+------+------
. . . |6 . 3 |. 7 .
5 . . |2 . . |. . .
1 . 4 |. . . |. . .

>>> ''.join(parse(puzzle1))
'4.....8.5.3..........7......2.....6.....8.4......1.......6.3.7.5..2.....1.4......'

>>> ''.join(parse(puzzle1)) == puzzle1
True

Now let's define what it means to be a solution: as we implied before, a solution is a grid where every unit is filled with a permutation of the digits 1 to 9, and the filled squares of the original puzzle remain unchanged in the solution:

def is_solution(grid, puzzle):
    """A grid is a solution to a puzzle if each unit is a permutation
    of the digits and if the filled squares remain unchanged."""
    filled_squares = [s for s in squares if puzzle[s] in digitset]
    def is_permutation(unit): return set(grid[s] for s in unit) == digitset
    return (grid is not False
            and all(is_permutation(u) for u in allunits)
            and all(grid[s]==puzzle[s] for s in filled_squares))

Searching for solutions

Humans solve Sudoku problems by using logical deduction to infer that a certain digit is (or is not) the value of a square; such deductions are repeated until all squares are filled with a single known digit. If the inference rules are sound and complete and the puzzle is well-formed, the human will only need a pencil--never an eraser. We could write an algorithm to do deductions like that, but the trick is in knowing when you have enough inference rules to solve every possible puzzle.

Therefore we will take a different route. Technically, it is known as model checking rather than deduction, and it consists of searching for a solution that works by trying out possibilities one by one. We will describe a series of search algorithms, starting with the hopelessly inefficient and ending with one that can solve over 5000 of the hardest Sudoku puzzles per second. However, when implemented by a human, this strategy requires an eraser as well as a pencil.

Algorithm 1: generate and test

First fill in all the empty squares with digits. Then check to see if we have a solution. If not, then fill the empty 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 (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 trying all 100¹⁰ combinations.

Here is the implementation of generate and test for Sudoku:

def solve(gridstring):
    "Generate every combination of digits for the empty squares; test each one."
    from itertools import product
    puzzle = parse(gridstring)
    for candidate in product(*[(v if v in digits else digits) for v in puzzle]):
        if is_solution(candidate, puzzle):
            return candidate
    return None

>>> tiny = '...4 4.3. .4.. .241'

>>> print show(parse(tiny))
. . |. 4 
4 . |3 . 
----+----
. 4 |. . 
. 2 |4 1 

>>> print show(solve(tiny))
1 2 |3 4
3 4 |1 2
----+----
4 3 |2 1
2 1 |4 3

It took 0.8 seconds to solve this tiny 4x4 puzzle. Now consider the 9x9 puzzle defined above as puzzle1. It has 17 filled squares, and thus 64 empty squares. So there are 9⁶⁴ possibile combinations (about 10⁶¹) to consider. How long will it take to run through them?

Suppose we have access to a secret new custom 1024 core 10 GHz CPU that features special hardware instructions that can generate and test a position in one clock cycle, and let's say we could afford to buy a billion of these, 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 10 billion galaxies and in each galaxy convince the inhabitants of 10 billion planets to each buy a billion of these computers, 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 less than 4% done with this one puzzle by now.

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

Algorithm 2: incremental search

Fortunately for us, Sudoku is more 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). 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 start with a 3.

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 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 the incremental search algorithm:

Start filling squares with digits, always making sure that the digit we put in a square is not the same as any peer. 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.

A few things to note about this search tree:

• Each successive level fills in one square, so with nine empty squares, any solution must appear only 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.
• A rectangle with a black X indicates an inconsistency: there is no way to continue because the next unfilled square already has all four values 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 value for the previous square. On the far-left branch of the tree we have to back up three levels before we find a square that has another possibility--we then erase the 1 in the upper-left corner and replace it with a 2.
• The two rightmost rectangles remain unexplored--the search proceeds left-to-right and has not got to them yet. (But since the puzzle is well-formed and we have found the solution, we know the unexplored branches will eventually lead to an inconsistency.)

The code to implement incremental search is straightforward:

def search(grid):
    """Select an unfilled square, try each possible value in order, recursively searching.
    When all squares filled: success; when no more digits to try: return False for failure."""
    if grid is False: return False
    s = select_square(grid)
    if s is None: return grid
    for d in possible_values(grid, s):
        result = search(assign(grid[:], s, d))
        if result: return result
    return False

def select_square(grid):
    "Return the first square that is not filled yet; or None if all squares are filled."
    for s in squares:
        if grid[s] == '.': return s
    return None

def possible_values(grid, s):
    "A square can be any digit that is not already taken by a peer."
    return digitset - set(grid[s] for s in peers[s])

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

def solve(gridstring):
    "Parse gridstring and search for a solution."
    return search(parse(gridstring))

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

>>> print show(solve(puzzle1))
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

However, if you wanted to solve a million puzzles, you'd prefer a faster algorithm. Look again at the search tree above, and consider how we can find the solution faster. Here are some general strategies:

1. We could try to move the solution to the left. 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 be no search--we would arrive directly at the solution.
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--the process of using the rules of Sudoku to conclude that a particular unfilled square must take on a certain value.
3. We could waste less time on the branches that eventually lead to an inconsistency (X). That is, we could notice earlier that a grid is destined to an inconsistency. This is another kind of constraint propagation, but one that concludes that a particular square has no possible value.
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 for any of the three reasons above.

Algorithm 3: constraint satisfaction with forward checking

Consider this branch of the search tree from above:

As soon as we place a 1 in the upper-left corner, there is an inconsistency: there is no value that can go in the square in the second row, second column (because 1 is in the same box, 2 is in the samecolumn, and 3 and 4 are in the same row). But our incremental search does not discover that until three more levels in the tree, when it gets around to selecting the offending square. The idea of forward checking is that we don't wait: every time we place a digit, we check all the unfilled squares and if one of them has no possible values we fail immediately rather than waiting to select the offending square.

In this case incremental search is not too bad; it generates two nodes in the search tree that could have been pruned. But in other cases this could be millions of nodes--forward checking is definitely worth it for Sudoku.

However, there is a cost to forward checking: we need to figure out if any non-selected square has no possible values left. To do that efficiently we will keep a collection of remaining possible values for every square; we can store the collection right in the grid. That is, grid[s] will hold a collection of all the remaining possible digits for square s. The collection could be represented by a Python set or list, but I chose instead to use a string of digits. So the value for the upper-left square in the tiny puzzle (before it is assigned a value) would be '123'.

We'll change solve so that it calls initialize to set all the unfilled squares to their initial value according to this convention. We'll also change assign(grid, s, d) so that it eliminates d from all the peers of s and fails if any peer has no remaining possible values. The eliminate function removes a digit from the possibilities for a square, and also does forward checking---checks if a peer is reduced to no values. (We add an extra level of indirection--the function check calls forward_check--because we will soon introduce other kinds of checks.) With these added complications one function becomes easier: possible_values can now just look up the collection of possible values rather than having to compute it. The functions search and select_square remain unchanged.

def solve(gridstring):
    "Parse gridstring, handle initial constraints, and search for a solution."
    return search(initialize(parse(gridstring)))

def initialize(puzzle):
    """First set each square in the result grid to have the set of all possible values.
    Then assign each filled square to the specified value."""
    grid = N * N * [digits]
    for s in squares:
        if puzzle[s] in digitset: grid = assign(grid, s, puzzle[s])
    return grid

def possible_values(grid, s): return grid[s]

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

def eliminate(grid, s, d):
    "Remove d from possibilities for grid[s]. If checking finds an inconsistency return False."
    if d not in grid[s]: return grid # Already eliminated d
    grid[s] = grid[s].replace(d, '')
    return check(grid, s, d)

def check(grid, s, d): return forward_check(grid, s)

def forward_check(grid, s):
    "Fail if s has no possible value, or one value which is the same as a peer."
    nvals = len(grid[s])
    return not (nvals==0 or (nvals==1 and any(grid[p]==grid[s] for p in peers[s])))

How much does forward checking help? For the easy puzzles, not much. One reason they are called "easy" is that they don't require forward checking. For the harder puzzle1 puzzle the solving speed is improved by a factor of 5. But that's still 43 seconds. Let's keep improving.

Algorithm 4: arc consistency

This idea is called arc consistency; it gets its name from the practice of describing a constraint satisfaction problem as a graph where the nodes are the variables (in Sudoku, the squares) and the arcs are the constraints. We make the arcs (the peer constraints) consistent one by one, and we can end up traversing multiple links in the graph (whereas the forward checking algorithm only looks ahead one arc in the graph).

The implementation of arc consistency for Sudoku is a simple change to what we already have:

def check(grid, s, d): return arc_consistent(grid, s)

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

Compared to forward checking, arc consistency reduces the time to solve puzzle1 from 43 seconds down to 0.16 seconds (a 266-fold speedup). The 50 easy problems are speeded up from a rate of 5 per second to 282 per second, a 56-fold speedup. In fact, 70% of the easy problems have their search tree pruned to a single node with this approach. Clearly, arc consistency is very effective on Sudoku problems (as it is on many constraint satisfaction problems), and one might be tempted to stop here. But let's try just a few more of the strategies that have proven popular for Sudoku (they're easy to implement).

Algorithm 5: dual consistency

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 False
    return True

Algorithm 6: naked pairs

def naked_pairs(grid, s):
    """If two squares s and s2 in a unit have the same two values,
    eliminate those values from every other square s3 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, s2)
               for s3 in u:
                    if s != s3 != s2:
                        if not all(eliminate(grid, s3, v) for v in vals):
                            return False
    return True

Benchmarks

Filename	N puzzles	Description
puzzle1.txt	1	hard puzzle, called puzzle1 above
easy50.txt	50	easy puzzles
hardest.txt	11	puzzles from around the web advertised as "hardest"
top95.txt	95	hard puzzles
top870.txt	870	hard puzzles
top2365.txt	2365	hard puzzles
subig20.txt	17,445	hard puzzles
min17.txt	49,151	hard puzzles with 17 filled squares
ALL.txt	69,091	all of the above (duplicates eliminated)

The code to gather statistics on the runtime of our algorithms on these puzzle sets can be seen in sudoku.py.

I ran the benchmark code on the five serious algorithms: incremental search (incr), forward checking (FC), and arc consistency by itself (AC), with dual consistency (AC+Dual), and with dual plus naked pairs (AC+D+NP). For each benchmark file we can see the name of the file, the number of puzzles, the average speed of solving (expressed two ways, as a frequency in Hertz and as the average number of seconds), the standard deviation of the results, and then, to give an idea of the distribution of times, the minimum, the first percentile, 25th percentile, median, 75th percentile, 99th percentile, and maximum times. Here are the results:

version     name        N       Hz      avg   stddev      min       1%      25%   median      75%      99%      max
   incr  puzzle1        1        0  225.729    0.000  225.729  225.729  225.729  225.729  225.729  225.729  225.729

     FC  puzzle1        1        0   42.132    0.000   42.132   42.132   42.132   42.132   42.132   42.132   42.132
     FC   easy50       50        5    0.172    0.272    0.004    0.004    0.010    0.050    0.161    0.987    0.987

     AC  puzzle1        1        6    0.159    0.000    0.159    0.159    0.159    0.159    0.159    0.159    0.159
     AC   easy50       50      285    0.004    0.005    0.001    0.001    0.001    0.002    0.003    0.025    0.025
     AC  hardest       11       22    0.044    0.113    0.001    0.001    0.002    0.010    0.014    0.401    0.401
     AC    top95       95        0    1.030    2.668    0.001    0.001    0.024    0.072    0.329   16.211   16.211
     AC   top870      870        4    0.221    1.111    0.001    0.001    0.006    0.020    0.060    6.943   16.197

AC+Dual  puzzle1        1       62    0.016    0.000    0.016    0.016    0.016    0.016    0.016    0.016    0.016
AC+Dual   easy50       50      188    0.005    0.001    0.005    0.005    0.005    0.005    0.005    0.009    0.009
AC+Dual  hardest       11       67    0.015    0.010    0.005    0.005    0.006    0.011    0.018    0.035    0.035
AC+Dual    top95       95       15    0.067    0.128    0.005    0.005    0.010    0.021    0.055    0.890    0.890
AC+Dual   top870      870       33    0.030    0.074    0.005    0.005    0.008    0.014    0.028    0.298    1.273

AC+D+NP  puzzle1        1      147    0.007    0.000    0.007    0.007    0.007    0.007    0.007    0.007    0.007
AC+D+NP   easy50       50      183    0.005    0.001    0.005    0.005    0.005    0.005    0.005    0.007    0.007
AC+D+NP  hardest       11       83    0.012    0.007    0.005    0.005    0.007    0.009    0.015    0.031    0.031
AC+D+NP    top95       95       30    0.033    0.063    0.005    0.005    0.007    0.012    0.025    0.489    0.489
AC+D+NP   top870      870       58    0.017    0.033    0.005    0.005    0.006    0.010    0.017    0.137    0.596

Algorithms 6-9: variable ordering

Here are the results; the benchmarks uphold the idea that MinRV is a good strategy; it performs the best on all the benchmarks except easy50, where it falls just slightly behind firstV. This is true on the average time to solution, and even more so on the maximum time.

version     name        N       Hz      avg   stddev      min       1%      25%   median      75%      99%      max
 firstV   easy50       50      183    0.005    0.001    0.005    0.005    0.005    0.005    0.005    0.007    0.007
 firstV  hardest       11       83    0.012    0.007    0.005    0.005    0.007    0.009    0.015    0.031    0.031
 firstV    top95       95       30    0.033    0.063    0.005    0.005    0.007    0.012    0.025    0.489    0.489
 firstV   top870      870       58    0.017    0.033    0.005    0.005    0.006    0.010    0.017    0.137    0.596

randomV   easy50       50      165    0.006    0.003    0.005    0.005    0.005    0.005    0.006    0.019    0.019
randomV  hardest       11      113    0.009    0.004    0.005    0.005    0.006    0.006    0.010    0.018    0.018
randomV    top95       95        7    0.127    0.203    0.005    0.005    0.021    0.047    0.120    1.090    1.090
randomV   top870      870       25    0.039    0.101    0.005    0.005    0.010    0.018    0.032    0.490    1.762

  MinRV   easy50       50      182    0.005    0.001    0.005    0.005    0.005    0.005    0.005    0.008    0.008
  MinRV  hardest       11      142    0.007    0.002    0.005    0.005    0.005    0.006    0.010    0.011    0.011
  MinRV    top95       95       57    0.017    0.017    0.005    0.005    0.007    0.012    0.021    0.132    0.132
  MinRV   top870      870       73    0.014    0.012    0.005    0.005    0.007    0.010    0.016    0.055    0.185

  MaxRV   easy50       50      164    0.006    0.002    0.005    0.005    0.005    0.005    0.005    0.017    0.017
  MaxRV  hardest       11       25    0.038    0.058    0.005    0.005    0.007    0.013    0.058    0.210    0.210
  MaxRV    top95       95        0    1.161    2.647    0.005    0.005    0.075    0.172    0.913   19.633   19.633
  MaxRV   top870      870        4    0.223    0.983    0.005    0.005    0.025    0.050    0.109    5.129   19.664

import random

def select_square(values):
    "Return an unfilled square with the fewest possible values; or None if no unfilled squares."
    def numvalues(s): return len(values[s])
    unfilled = [s for s in squares if numvalues(s) > 1]
    return None if unfilled==[] else random.choice(unfilled)

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

def select_square(values):
    "Return an unfilled square with the fewest possible values; or None if no unfilled squares."
    def numvalues(s): return len(values[s])
    unfilled = [s for s in squares if numvalues(s) > 1]
    return None if unfilled==[] else max(unfilled, key=numvalues)

Algorithm 10-12: value ordering

def constraining_value(values, s, d):
    "Count the number of peers of s that would have 0 or 1 value left without d."
    return sum((d in values[s2] and len(values[s2]) <= 2) for s2 in peers[s])

def possible(values, s):
    return sorted(values[s], key=lambda d: constraining_value(values, s, d))


def possible(values, s):
    return sorted(values[s], key=lambda d: -constraining_value(values, s, d))

def possible(values, s):
    vals = list(values[s])
    random.shuffle(vals)
    return vals

benchmark([1,2,3,4], "RandVal")

 version    name        N       Hz      avg   stddev      min       1%      25%   median      75%      99%      max
MostCVl   easy50       50      184    0.005    0.000    0.005    0.005    0.005    0.005    0.005    0.008    0.008
MostCVl  hardest       11      113    0.009    0.006    0.005    0.005    0.006    0.006    0.009    0.029    0.029
MostCVl    top95       95       57    0.018    0.018    0.005    0.005    0.008    0.012    0.021    0.137    0.137
MostCVl   top870      870       72    0.014    0.013    0.005    0.005    0.007    0.011    0.016    0.057    0.193

LeastCV   easy50       50      181    0.006    0.001    0.005    0.005    0.005    0.005    0.005    0.009    0.009
LeastCV  hardest       11      124    0.008    0.004    0.005    0.005    0.006    0.006    0.010    0.018    0.018
LeastCV    top95       95       54    0.018    0.018    0.005    0.005    0.008    0.012    0.022    0.135    0.135
LeastCV   top870      870       72    0.014    0.013    0.005    0.005    0.008    0.011    0.016    0.058    0.194

RandVal   easy50       50      185    0.005    0.000    0.005    0.005    0.005    0.005    0.005    0.008    0.008
RandVal  hardest       11       91    0.011    0.007    0.005    0.005    0.006    0.007    0.014    0.029    0.029
RandVal    top95       95       52    0.019    0.018    0.005    0.005    0.008    0.013    0.023    0.123    0.123
RandVal   top870      870       75    0.013    0.012    0.005    0.005    0.007    0.010    0.015    0.056    0.209

Robustness

version     name        N       Hz      avg   stddev      min       1%      25%   median      75%      99%      max
 Normal    top95       95       38    0.026    0.025    0.008    0.008    0.011    0.018    0.032    0.196    0.196
 Normal   top870      870       49    0.020    0.018    0.008    0.008    0.011    0.016    0.023    0.082    0.275
 Normal  top2365     2365       54    0.018    0.018    0.008    0.008    0.010    0.014    0.021    0.076    0.435
 Normal  subig20    17445       83    0.012    0.022    0.007    0.007    0.008    0.008    0.010    0.066    1.583
 Normal    min17    49151       49    0.020    0.052    0.007    0.007    0.008    0.010    0.016    0.179    3.385
 Normal     200K   200000       47    0.021    0.051    0.007    0.007    0.008    0.011    0.017    0.178    3.732

1 . . |. . . |. 6 4        1 3 2 |8 5 7 |9 6 4        Solution in 3.732 seconds
. . . |3 2 . |. . .        9 6 4 |3 2 1 |5 7 8        
. 5 . |. . . |. . .        7 5 8 |9 6 4 |1 3 2        
------+------+------       ------+------+------        
. 2 . |6 4 . |. . .        5 2 1 |6 4 8 |7 9 3        
6 . . |. . . |. . 1        6 8 9 |5 7 3 |2 4 1        
. . . |. . . |8 . .        3 4 7 |2 1 9 |8 5 6        
------+------+------       ------+------+------        
. . . |. . 5 |3 2 .        8 1 6 |4 9 5 |3 2 7        
4 . . |7 . . |. . .        4 9 3 |7 8 2 |6 1 5        
. . . |. . . |. . .        2 7 5 |1 3 6 |4 8 9    

version     name        N     max  N-ratio max-ratio
 Normal    top95       95   0.196
 Normal   top870      870   0.275  9.2      1.4
 Normal  top2365     2365   0.435  2.7      1.6
 Normal  subig20    17445   1.583  7.4      3.6
 Normal    min17    49151   3.385  2.8      2.1
 Normal     200K   200000   3.732  4.0      1.1

Algorithm 13: Concurrency

Therefore I decided to port the program to Java, which has good facilities for concurrency. The translation is fairly straightforward, but there are some changes; some because of the nature of Java, and some because I wanted to take advantage of this rewrite opportunity to make some choices that will enhance efficiency. Here are the main differences between the Python version and the Java version:

• The Java version uses concurrency. I decided to do this in the simplest way possible: Given a benchmark file of 100,000 puzzles, and 10 threads, just divide the list of puzzles into 10 sublists of 10,000 puzzles each, and give one to each thread. If the variance of solving puzzles were larger we would want a more complicated scheme that allows a thread that finishes early to take some work from a thread that is not yet finished. But that would complicate the program, and it seems that the variance is low enough that it is not necessary.
• Java's compiler should generate much more efficient code than Python's interpreter.
• Java doesn't have Python's concise generator expressions and list comprehensions, so the code will be more verbose, and will rely more on for loops and less on expressions.
• Java has good facilities at the bit-twiddling level, and we will take advantage of that in the representation of sets of possible digits within a grid square. For example, the set of values {1, 4, 9}, which was represented in Python by the string '149', will be represented in Java by the binary bit pattern 100001001 (which is the int 265, or 1+8+256).
• Java is a less dynamic language than Python. In Python a list can be extended in length, and can hold objects of any type. A Java array cannot change size, and must be declared to hold a specific type; for example with "int[81] grid". Python's flexibility comes with a performance penalty. Checking to see if a digit is a possibility for some grid square with the code d in grid[s] means at run time Python will have to fetch the grid object, check its type to see how to do indexing, fetch the underlying array from the list object, fetch the s object, get the numeric value of the s object, fetch grid[s], check its type (str) to see how to do the membership test, and then test if the string d is a substring of grid[s]. In Java, the corresponding code is (d & grid[s]) > 0), and at run time Java just has to check that s is less than 81, fetch the memory location of grid plus s, do a bitwise and with the value od d, and compare to 0.
• In the Python version, we made a fresh copy of the grid every time we recursively called search, by executing grid[:]. Over the course of a benchmark Python needs to allocate (and reclaim) storage for millions of these grids. In Java we will save the time of managing memory by pre-allocating 81 grids (one for each level in the search tree) in an array called gridpool, and then instead of creating a new grid, we can just copy the current grid into the pre-allocated grid at the current level. We still need to do the same number of copies, but we save the time of managing storage.

When we start adding threads, the performance increases almost linearly with the number of threads up to 6 (we reach a spedup of 5.66 and 5.36 in two runs with 6 threads). That is what I would expect with 6 CPUs. However, as we add more threads the gains increase at about half the slope up to 10 threads (where we have a speedup of 7.20 or 7.26 in two runs). From there we have a long flat plateau where the speedup is between 7 and 7.5 no matter how many threads we add. (This is 700 times faster than the Python program.) The top measured performance happens to be at 40 threads, with 43,873 puzzles solved per second, but this is more due to random variation than to any special property of the number 40. (if we repeated the experiment, the peak might come at 20, or 12. There is variation because the computer is doing other things (I didn't shut down all other processes) and because the Java JIT compiler performs slightly differently each run. I did let the JIT "warm up" by solving a few hundred thousand puzzles before I started the measurements.)

Alternative approaches

• Heuristic search: traditional heuristic search, as exemplified by the A* Algorithm>/a> is a "best-first" approach that skips around in the search tree. Our search was depth-first, meaning we always went as deep as possible before backtracking. A* keeps a list of all states that have not been proven inconsistent and looks first at the one that is estimated to be on the shortest path to the goal. This works well for things like finding a route to a destination on a map, but one reason it does not seem promising for Sudoku is that all Sudoku solutions have the same path length, so it is difficult to come up with a useful heuristic.
• Local search: so far we have considered incremental search, in which we start with a partial assignment of digits to squares and incrementally add new assignments, never violating the rules of Sudoku. In a local search approach we start with a complete (random) assignment of digits to all squares and expect that it violates the rules (for example, there might be two instances of a '3' in the same row). Then we make small (local) changes to the assignment so that it violates fewer rules. Local search is good when it is sufficient to find an answer that is an approximation to the best answer, and for problems that have a large number of solutions (with a very large number of non-solutions). Since these conditions do not apply to Sudoku, my guess is that local search would not perform well. Local search techniques include hill climbing, simulated annealing, and genetic algorithms.
• Dancing links: this algorithm has proven to be effectivce for Sudoku. It is a constraint-satisfaction approach that considers both the squares and the locations as variables. In other words, it implements our dual consistency rule by keeping track of variables, rather than by computing places for a digit on the fly. It might well be more efficient than what I have implemented, but it didn't seem that interesting to me, since it is doing roughly the same thing.

Don't let my reluctance stop you: it would be a great project to implement any of these techniques for Sudoku.

Why?

Translations