Peter Norvig, Feb 2021

CrossProduct Puzzle¶

The 538 Riddler poses a type of puzzle called CrossProduct, which works like this:

Replace each "?" in the table with a single digit so that the product of the digits in each row equals the number to the right of the row, and the product of the digits in each column equals the number above the column.

Sample puzzle:

6615	15552	420	[6×3]
?	?	?	210
?	?	?	144
?	?	?	54
?	?	?	135
?	?	?	4
?	?	?	49

Solution:

6615	15552	420	[6×3]
7	6	5	210
9	8	2	144
3	9	2	54
5	9	3	135
1	4	1	4
7	1	7	49

We could solve CrossProduct puzzles by hand, but why not write a program to do it?

Data type definitions¶

Here are the data types we will use in trying to solve CrossProduct puzzles:

Digit: a single digit, from 1 to 9 (but not 0).
Row: a tuple of digits that forms a row in the table, e.g. (7, 6, 5).
Table: a table of digits that fill in all the "?"s; a list of rows, e.g. [(7, 6, 5), (9, 8, 2), ...].
Products: a list of the numbers that corresponding digits must multiply to, e.g. in the puzzle above:
- [6615, 15552, 420] for the column products;
- [210, 144, 54, 135, 4, 49] for the row products.
Puzzle: a puzzle to be solved, as defined by the row products and column products.

In [1]:

from typing      import Tuple, List, Set, Iterable, Optional
from numpy       import divide, prod, transpose
from collections import namedtuple
import random

Digit    = int
Row      = Tuple[Digit, ...] 
Table    = List[Row]   
Product  = int
Products = List[Product] 
Puzzle   = namedtuple('Puzzle', 'row_prods, col_prods')

The puzzles¶

Here are the puzzles given by 538 Riddler (they promised one a week for four weeks; so far we've seen three):

In [2]:

puzzles = (
    Puzzle([135,  45,  64, 280, 70],             [3000,   3969,    640]),
    Puzzle([210, 144,  54, 135,  4,  49],        [6615,  15552,    420]),
    Puzzle([280, 168, 162, 360, 60, 256, 126], [183708, 245760, 117600]))

Strategy¶

Here's my strategy:

To solve a puzzle, first find all ways to fill the first row, and for each way, solve the rest of the puzzle.
To fill a row, first find all ways to fill the first digit, and for each way, fill the rest of the row.

So the first step is to define fill_one_row(row_prod, col_prods) to return a set of digit-tuples that can legally fill a row that has the given row product in a puzzle with the given column products.

If col_prods is empty, then there is one solution (the 0-length tuple) if row_prod is 1, and no solution otherwise.
Otherwise, try each digit d that divides both the row_prod and the first number in col_prods, and then try all ways to fill the rest of the row.

In [3]:

def fill_one_row(row_prod: Product, col_prods: Products) -> Set[Row]:
    "All permutations of digits that multiply to `row_prod` and evenly divide `col_prods`."
    if not col_prods:
        return {()} if row_prod == 1 else set()
    else:
        return {(d, *rest) for d in range(1, 10)
                if (row_prod / d).is_integer() and (col_prods[0] / d).is_integer()
                for rest in fill_one_row(row_prod // d, col_prods[1:])}

Some examples:

In [4]:

fill_one_row(210, [6615, 15552, 420]) # There are 2 ways to fill this row

Out[4]:

{(5, 6, 7), (7, 6, 5)}

In [5]:

fill_one_row(54, [6615, 15552, 420]) # There are 8 ways to fill this row

Out[5]:

{(1, 9, 6),
 (3, 3, 6),
 (3, 6, 3),
 (3, 9, 2),
 (9, 1, 6),
 (9, 2, 3),
 (9, 3, 2),
 (9, 6, 1)}

Now we can solve the rest of a puzzle:

solve(puzzle) finds the first solution. (A well-formed puzzle has exactly one solution, but some might have more, or none.)
solutions(puzzle) yields all possible solutions to a puzzle. There are three main cases to consider:
- A puzzle with no rows has the empty table, [], as a solution, as long as the column products are all 1.
- A puzzle with rows might have solutions, as long as the column products are all integers. Call fill_row to get all possible ways to fill the first row, and for each one recursively call solutions to get all the possible ways of filling the rest of the rows (making sure to pass in an altered col_prods where each element is divided by the corresponding element in the first row).
- Otherwise there are no solutions.

In [6]:

def solve(puzzle) -> Optional[Table]: return next(solutions(puzzle), None)

def solutions(puzzle) -> Iterable[Table]:
    """Yield all tables that solve the puzzle.
    The product of the digits in row r must equal row_prods[r], for all r.
    The product of the digits in column c must equal col_prods[c], for all c."""
    row_prods, col_prods = puzzle
    if not row_prods and all(c == 1 for c in col_prods):
        yield []
    elif row_prods and all(c == int(c) for c in col_prods):
        yield from ([row1, *rows]
                    for row1 in fill_one_row(row_prods[0], col_prods)
                    for rows in solutions(Puzzle(row_prods[1:], list(divide(col_prods, row1)))))

Solutions¶

Here are solutions to the three puzzles posed by The Riddler:

In [7]:

[solve(p) for p in puzzles]

Out[7]:

[[(3, 9, 5), (5, 9, 1), (8, 1, 8), (5, 7, 8), (5, 7, 2)],
 [(7, 6, 5), (9, 8, 2), (3, 9, 2), (5, 9, 3), (1, 4, 1), (7, 1, 7)],
 [(7, 8, 5), (3, 8, 7), (9, 6, 3), (9, 8, 5), (3, 5, 4), (4, 8, 8), (9, 2, 7)]]

Those are the correct solutions. However, we could make them look nicer.

Prettier solutions¶

In [8]:

from IPython.display import Markdown, display

def pretty(puzzle) -> Markdown:
    """A puzzle and its solution in pretty Markdown format."""
    row_prods, col_prods = puzzle
    head = row(col_prods + [f'[{len(row_prods)}×{len(col_prods)}]'])
    dash = row(['---'] * (1 + len(col_prods)))
    rest = [row(r + (f'**{rp}**',))
            for r, rp in zip(solve(puzzle), row_prods)]
    return Markdown('\n'.join([head, dash, *rest]))

def row(items) -> str: 
    """Make a markdown table row."""
    return '|' + '|'.join(map(str, items)) + '|'

In [9]:

for p in puzzles:
    display(pretty(p))

3000	3969	640	[5×3]
3	9	5	135
5	9	1	45
8	1	8	64
5	7	8	280
5	7	2	70

6615	15552	420	[6×3]
7	6	5	210
9	8	2	144
3	9	2	54
5	9	3	135
1	4	1	4
7	1	7	49

183708	245760	117600	[7×3]
7	8	5	280
3	8	7	168
9	6	3	162
9	8	5	360
3	5	4	60
4	8	8	256
9	2	7	126

Making new puzzles¶

Can we make new puzzles? Can we make well-formed ones (those with exactly one solution)? Here is an approach:

Make a table filled with random digits (random_table).
Make a puzzle from the row and column products of the table (table_puzzle).
Repeat N times (random_puzzles).
Optionally, check if puzzles are well-formed.

In [10]:

def random_table(nrows, ncols) -> Table:
    "Make a table of random digits of the given size."
    return [tuple(random.randint(1, 9) for c in range(ncols))
            for r in range(nrows)]

def table_puzzle(table) -> Puzzle:
    "Given a table, compute the puzzle it is a solution for."
    return Puzzle([prod(row) for row in table], 
                  [prod(col) for col in transpose(table)])

def random_puzzles(N, nrows, ncols, seed=42) -> List[Puzzle]: 
    "Return a list of `N` random puzzles."
    random.seed(seed) # For reproducability
    return [table_puzzle(random_table(nrows, ncols)) for _ in range(N)]

def well_formed(puzzle) -> bool: 
    "Does the puzzle have exactly one solution?"
    S = solutions(puzzle)
    first, second = next(S, None), next(S, None)
    return first is not None and second is None

In [11]:

random_table(nrows=5, ncols=3)

Out[11]:

[(6, 1, 7), (9, 7, 3), (9, 3, 5), (5, 6, 9), (6, 6, 4)]

In [12]:

puz = table_puzzle(_)

In [13]:

well_formed(puz)

Out[13]:

False

In [14]:

len(list(solutions(puz)))

Out[14]:

In [15]:

pretty(puz)

Out[15]:

14580	756	3780	[5×3]
6	1	7	42
9	7	3	189
5	3	9	135
9	6	5	270
6	6	4	144

How likely are random puzzles (of various sizes) to be well-formed?

In [16]:

N = 200
for r, c in [(3, 3), (3, 4), (4, 3), (3, 5), (5, 3), (4, 4), (6, 3)]:
    w = sum(map(well_formed, random_puzzles(N, r, c))) / N
    print(f'{w:3.0%} of random puzzles with {r} rows and {c} cols ({r * c:2} cells) are well-formed')

33% of random puzzles with 3 rows and 3 cols ( 9 cells) are well-formed
18% of random puzzles with 3 rows and 4 cols (12 cells) are well-formed
14% of random puzzles with 4 rows and 3 cols (12 cells) are well-formed
 8% of random puzzles with 3 rows and 5 cols (15 cells) are well-formed
 2% of random puzzles with 5 rows and 3 cols (15 cells) are well-formed
 4% of random puzzles with 4 rows and 4 cols (16 cells) are well-formed
 0% of random puzzles with 6 rows and 3 cols (18 cells) are well-formed

We see that most puzzles are not well-formed. Smaller sizes are more likely to yield well-formed puzzles.

Speed¶

How long does it take to solve random puzzles? We can do a thousand small (5x3) puzzles in about two seconds:

In [17]:

%time all(solve(p) for p in random_puzzles(1000, 5, 3))

CPU times: user 1.56 s, sys: 48.6 ms, total: 1.6 s
Wall time: 1.56 s

Out[17]:

True

Puzzles that are even a little bit larger can be a lot slower, and there is huge variability in the time to solve. For example, a single 10 x 6 puzzle can take from a few milliseconds to tens of seconds:

In [18]:

[p10x6] = random_puzzles(1, 10, 6)
%time pretty(p10x6)

CPU times: user 3.03 s, sys: 10.3 ms, total: 3.04 s
Wall time: 3.03 s

Out[18]:

24576	979776	274400	2177280	1792000	524880	[10×6]
4	1	5	6	2	2	480
1	7	1	3	4	3	252
8	9	2	9	2	1	2592
8	3	7	8	5	6	40320
1	6	7	3	5	3	1890
1	8	1	7	4	6	1344
3	9	2	8	5	6	12960
4	6	5	1	7	9	7560
1	1	8	2	4	5	320
8	2	7	5	8	3	13440

In general, the time to solve a puzzle can grow exponentially in the number of cells. Consider a row in a six-column puzzle, where the products are all 5040. There are 3,960 ways to fill this row:

In [19]:

n = 5040
len(fill_one_row(n, [n] * 6))

Out[19]:

If four rows all had a similar number of possibilities and didn't constrain each other, that would be hundreds of trillions of combinations to try—an infeasible number. We will need a faster algorithm for larger puzzles.

Faster Speed¶

To speed things up, we could encode the puzzle as a constraint satisfaction problem (CSP), and use a highly-optimized CSP solver. But even without going to a professional-grade CSP solver, we could borrow the heuristics they use. There are four main considerations in CSP solving:

Variable definition: In solutions, we are treating each row as a variable, and asking "which of the possible values returned by fill_one_row will work as the value of this row? An alternative would be to treat each cell as a variable, and fill in the puzzle one cell at a time rather than one row at a time. This has the advantage that each variable has only 9 possible values, not thousands of possibilities.
Variable ordering: In solutions, we consider the variables (the rows) in strict top-to-bottom order. It is usually more efficient to reorder the variables, filling in first the variable with the minimum number of possible values. The reasoning is that if you have a variable with only 2 possibilities, you have a 50% chance of guessing right the first time, whereas if there were 100 possibilities, you have only a 1% chance of guessing right.
Value ordering: The function fill_one_row returns values in sorted lexicographic order, lowest first. We could reorder the values to pick the one that imposes the least constraints first (that is, the value that allows the most possibilities for the other variables).
Domain-specific heuristics: CSP solvers are general, but sometimes knowledge that is specific to a problem can be helpful. One fact about CrossProduct is that the digits 5 and 7 are special in the sense that if a row (or column) product is divisible by 5 (or 7), then the digit 5 (or 7) must appear in the row (or column). That is not true for the other digits (for example, if a row product is divisible by 8, then an 8 may appear in the row, or it might be a 2 and a 4, or three 6s, etc.).

Usually variable ordering is the most productive heuristic. Let's try it. The function reorder takes a puzzle and returns a version of the puzzle with the row products permuted so that the rows with the fewest possible fillers come first:

In [20]:

def reorder(puzzle) -> Puzzle:
    """Create a version of puzzle with the rows reordered so the rows with the fewest
    number of possible fillers come first."""
    def fillers(r): return len(fill_one_row(r, puzzle.col_prods))
    rows = sorted(puzzle.row_prods, key=fillers)
    return Puzzle(rows, puzzle.col_prods)

In [21]:

p2 = puzzles[2]
p2, reorder(p2)

Out[21]:

(Puzzle(row_prods=[280, 168, 162, 360, 60, 256, 126], col_prods=[183708, 245760, 117600]),
 Puzzle(row_prods=[256, 280, 162, 360, 126, 168, 60], col_prods=[183708, 245760, 117600]))

In [22]:

# How many ways are there to fill each row?
{r: len(fill_one_row(r, p2.col_prods)) 
 for r in reorder(p2).row_prods}

Out[22]:

{256: 1, 280: 2, 162: 2, 360: 2, 126: 5, 168: 7, 60: 8}

Now I'll define a set of test puzzles and see how long it takes to solve them all, and compare that to the time to solve the reordered versions:

In [23]:

test_puzzles = random_puzzles(20, 10, 3)

In [24]:

%time all(solve(p) for p in test_puzzles)

CPU times: user 20.6 s, sys: 453 ms, total: 21.1 s
Wall time: 20.7 s

Out[24]:

True

In [25]:

%time all(solve(reorder(p)) for p in test_puzzles)

CPU times: user 134 ms, sys: 3.41 ms, total: 137 ms
Wall time: 134 ms

Out[25]:

True

That's a nice improvement—150 times faster on this small test set! I'm curious whether we would get even more speedup by treating each cell as a separate variable, or by considering value ordering, but I'll leave that as an exercise for the reader.

Tests¶

A suite of unit tests:

In [26]:

def test():
    "Test suite for CrossProduct functions."
    assert fill_one_row(1, [])  == {()}
    assert fill_one_row(2, [])  == set()
    assert fill_one_row(9, [9]) == {(9,)}
    assert fill_one_row(10, [10])  == set()
    assert fill_one_row(73, [360, 360, 360])  == set()
    
    assert solve(Puzzle([], []))   == []
    assert solve(Puzzle([], [1]))  == []
    assert solve(Puzzle([], [2]))  == None
    assert solve(Puzzle([5], [5])) == [(5,)]
    assert solve(Puzzle([0], [0])) == None # Maybe should allow zero as a digit?
    assert solve(Puzzle([2, 12], [3, 8])) == [(1, 2), (3, 4)]

    assert fill_one_row(729, [90, 126, 81]) == {(9, 9, 9)} # Unique fill
    
    assert fill_one_row(729, [90, 126, 81, 30]) == {
     (3, 9, 9, 3), (9, 3, 9, 3), (9, 9, 3, 3), (9, 9, 9, 1)}
    
    # 72 has the most ways to fill a 3-digit row
    assert max(range(1, 100), key=lambda n: len(fill_one_row(n, [5*7*8*9]*3))) == 72
    assert fill_one_row(72, [72, 72, 72])  == { 
        (1, 8, 9),
        (1, 9, 8),
        (2, 4, 9),
        (2, 6, 6),
        (2, 9, 4),
        (3, 3, 8),
        (3, 4, 6),
        (3, 6, 4),
        (3, 8, 3),
        (4, 2, 9),
        (4, 3, 6),
        (4, 6, 3),
        (4, 9, 2),
        (6, 2, 6),
        (6, 3, 4),
        (6, 4, 3),
        (6, 6, 2),
        (8, 1, 9),
        (8, 3, 3),
        (8, 9, 1),
        (9, 1, 8),
        (9, 2, 4),
        (9, 4, 2),
        (9, 8, 1)}
    
    assert fill_one_row(7**8, [7]*9) == {
        (1, 7, 7, 7, 7, 7, 7, 7, 7),
        (7, 1, 7, 7, 7, 7, 7, 7, 7),
        (7, 7, 1, 7, 7, 7, 7, 7, 7),
        (7, 7, 7, 1, 7, 7, 7, 7, 7),
        (7, 7, 7, 7, 1, 7, 7, 7, 7),
        (7, 7, 7, 7, 7, 1, 7, 7, 7),
        (7, 7, 7, 7, 7, 7, 1, 7, 7),
        (7, 7, 7, 7, 7, 7, 7, 1, 7),
        (7, 7, 7, 7, 7, 7, 7, 7, 1)}
    
    assert solve(Puzzle([210, 144, 54, 135, 4, 49], [6615, 15552, 420])) == [
        (7, 6, 5), 
        (9, 8, 2), 
        (3, 9, 2), 
        (5, 9, 3), 
        (1, 4, 1), 
        (7, 1, 7)]
    
    assert sorted(solutions(Puzzle([8, 8, 1], [8, 8, 1]))) == [ # Multi-solution puzzle
        [(1, 8, 1), 
         (8, 1, 1), 
         (1, 1, 1)],
        [(2, 4, 1), 
         (4, 2, 1), 
         (1, 1, 1)],
        [(4, 2, 1), 
         (2, 4, 1), 
         (1, 1, 1)],
        [(8, 1, 1), 
         (1, 8, 1), 
         (1, 1, 1)]]
    
    assert not list(solutions(Puzzle([8, 8, 1], [8, 8, 2]))) # Unsolvable puzzle
    
    assert solve(Puzzle([1470, 720, 270, 945, 12, 343], 
                        [6615, 15552, 420, 25725])) == [ # 4 column puzzle
        (7, 6, 5, 7),
        (9, 8, 2, 5),
        (3, 9, 2, 5),
        (5, 9, 3, 7),
        (1, 4, 1, 3),
        (7, 1, 7, 7)]
    
    puzz  = Puzzle([6, 120, 504], [28, 80, 162])
    table = [(1, 2, 3), 
             (4, 5, 6), 
             (7, 8, 9)]
    assert solve(puzz) == table
    assert table_puzzle(table) == puzz
    assert well_formed(puzz)
    
    assert not well_formed(Puzzle([7, 7], [7, 7]))
    assert well_formed(Puzzle([64, 224, 189, 270, 405, 144, 105], 
                              [308700, 12960, 1119744]))
    
    assert row((1, 2, 3)) == '|1|2|3|'
    
    col_prods = [193536, 155520, 793800]
    assert (reorder(Puzzle([10, 48,  36,  7, 32,  81, 252, 160, 21, 90], col_prods)) == 
                    Puzzle([ 7, 10, 160, 21, 81, 252,  90,  32, 48, 36], col_prods))
    return True
    
test()

Out[26]:

True