Riddler: Can You Hunt For The Mysterious Numbers?

2021-01-22

This week’s FiveThirtyEight Riddler Classic offers a sudoku-like puzzle:

There are eight three-digit numbers — each belongs in a row of the table below, with one digit per cell. The products of the three digits of each number are shown in the rightmost column. Meanwhile, the products of the digits in the hundreds, tens, and ones places, respectively, are shown in the bottom row. Can you find all eight three-digit numbers and complete the table? It’s a bit of a mystery, but I’m sure you have it within you to hunt down the answer!

Linear programming

A nice way to solve problems like this is by Integer Linear Programming. Like in my Trees and Tents post, I like to do this in Julia using the JuMP package.

The constraints we have are on the products of the numbers in the rows and the columns. Therefore the problem is not linear if we implement it naively (a variable for every cell denoting the number between 1 and 9 it should contain). We need to do something a little bit more clever. A natural way to look at this problem is via the prime factorizations of the row and the column products. There are four primes between 1 and 9, namely 2, 3, 5 and 7. Given the prime factorization of every row and column, we only need to distribute these primes over the cells. Since \(294 = 2^1\cdot 3^1\cdot 5^0 \cdot 7^2\), we know that in the first row there should be a 2, a 3, no 5 and two 7s. Note that a cell can contain multiple primes, or none at all (meaning it is a 1).

We therefore need a variable for every cell and each of the 4 possible primes, the value of which denotes the multiplicity of this prime in this cell.

In Julia, we first import the required packages, and define the row and column products.

using JuMP, GLPK
using Primes: factor, primes
using DataStructures: DefaultDict

row_prods = [294, 216, 135, 98, 112, 84, 245, 40];
col_prods = [8_890_560, 156_800, 55_566];
prms = primes(1, 9);

n_prms = length(prms);
n_rows = length(row_prods);
n_cols = length(col_prods);

Now, we create lists of the multiplicities of the four primes in the rows and the columns. The [1, 1, 0, 2] denotes the multiplicities of 2, 3, 5 and 7 in the first row.

row_factors = factor.(Dict, row_prods);
col_factors = factor.(Dict, col_prods);

row_factors = [[DefaultDict(0, row_factors[row])[prms[i]] for i in 1:4] for row in 1:n_rows]
col_factors = [[DefaultDict(0, col_factors[col])[prms[i]] for i in 1:4] for col in 1:n_cols]
8-element Array{Array{Int64,1},1}:
 [1, 1, 0, 2]
 [3, 3, 0, 0]
 [0, 3, 1, 0]
 [1, 0, 0, 2]
 [4, 0, 0, 1]
 [2, 1, 0, 1]
 [0, 0, 1, 2]
 [3, 0, 1, 0]
3-element Array{Array{Int64,1},1}:
 [6, 4, 1, 3]
 [7, 0, 2, 2]
 [1, 4, 0, 3]

Now onto the actual problem definition. As stated, we need a variable for every cell and every prime. Every prime can occur 0 or more times, and its multiplicity is an integer.

model = Model(with_optimizer(GLPK.Optimizer))

@variable(model, 0 <= x[1:n_rows, 1:n_cols, 1:n_prms], integer=true)

Rows and column prime multiplicities should add up to the correct numbers.

for i in 1:n_rows, p in 1:n_prms
    @constraint(model, sum(x[i, :, p]) == row_factors[i][p])
end

for j in 1:n_cols, p in 1:n_prms
    @constraint(model, sum(x[:, j, p]) == col_factors[j][p])
end

Finally, we want every cell to be an integer between 1 and 9, which means \(2^{p_2}\cdot 3^{p_3}\cdot 5^{p_5}\cdot 7^{p_7} \le 9\), or \(p_2\log(2) + p_3\log(3) + p_5\log(5)+p_7\log(7) \le \log(9)\), a constraint which is linear in its variables.

for i in 1:n_rows, j in 1:n_cols
    @constraint(model, sum(x[i, j, :] .* log.(prms)) <= log(9))
end

And that is all! We let JuMP do the hard work and come up with a feasible solution. As an example, we print out the number of times a 2 (the first prime) appears in every cell.

JuMP.optimize!(model) 

solution = JuMP.value.(x);

solution[:, :, 1]
8×3 Array{Float64,2}:
 0.0  0.0  1.0
 0.0  3.0  0.0
 0.0  0.0  0.0
 0.0  1.0  0.0
 3.0  1.0  0.0
 0.0  2.0  0.0
 0.0  0.0  0.0
 3.0  0.0  0.0

The final solution to the riddle is then the product of every prime to the power of its multiplicity:

[[Int(prod(prms .^ solution[i, j, :])) for j in 1:n_cols] for i in 1:n_rows]
8-element Array{Array{Int64,1},1}:
 [7, 7, 6]
 [9, 8, 3]
 [9, 5, 3]
 [7, 2, 7]
 [8, 2, 7]
 [7, 4, 3]
 [5, 7, 7]
 [8, 5, 1]

And there’s the solution! I love how easy it is to have JuMP solve such a problem after coming up with a way to phrase it as a linear program.