Crossing the CrossProduct off my list

OK, so I know that I may be jumping the gun here since there is another week of CrossProduct yet to come, so I am sorry. SPOILER ALERT!!!! So, you know avert your eyes, if you want to keep solving CrossProduct the old fashioned way with pen and paper, which I definitely did for the first two weeks and even for a first stab at this week's. However, as fun as this logic puzzle has been, I figured that I could outsource the heavy lifting to a much better logician with more brain power ... so I built a generic solver.

The general gist of the solver is not too dissimilar from the way that I typically attacked the problem with pen and paper: use the prime factorizations of the inputs to recast the puzzle as an integer constraint satisfaction problem. Assume we have a CrossProduct puzzle with $I$ rows and $J$ columns.

For each element $n_i$ in the input column, we can write the unique prime factorization as $$n_i = 2^{A_i} 3^{B_i} 5^{C_i} 7^{D_i}, \,\, i=1, \dots, I.$$ We know that there can be no higher factors by design since the puzzle is to write $n_i$ as the product of three single digit numbers. Similarly, for each element $m_j$ in the input row, we can write the unique prime factorization as $$m_j = 2^{\alpha_j} 3^{\beta_j} 5^{\gamma_j} 7^{\delta_j}, \,\, j= 1, \dots, J.$$

Let $u_{ij}, v_{ij}, w_{ij}, x_{ij}$ be nonnegative integers for $i \in \{1, 2, \dots, I\}, j \in \{1, 2, \dots, J\}.$ We want to model the entry in the $i$th row and $j$th column as $$f_{ij} = 2^{u_{ij}} 3^{v_{ij}} 5^{w_{ij}} 7^{x_{ij}},$$ such that following product conditions hold $$\prod_{j=1}^J f_{ij} = n_i, \prod_{i=1}^I f_{ij} = m_j, \,\, i=1,\dots,I, \, j=1, \dots, J.$$

The condition that $f_{ij}$ is a single digit number can be translated to the linear condition $$\ln(2) u_{ij} + \ln(3) v_{ij} + \ln(5) w_{ij} + \ln(7) x_{ij} \leq \ln(9), \,\, i = 1, \dots, I, j = 1, \dots, J.$$ Similarly, the conditions that the products across rows and columns are equal to the values in the input column and row, respectively, are given by \begin{align*} \sum_{j=1}^J u_{ij} &= A_i,\,\, i=1, \dots, I, \,\, &\sum_{i=1}^I u_{ij} = \alpha_j, \,\, j =1, \dots, J\\ \sum_{j=1}^J v_{ij} &= B_i,\,\, i=1, \dots, I, \,\, &\sum_{i=1}^I v_{ij} = \beta_j, \,\, j =1, \dots, J\\ \sum_{j=1}^J w_{ij} &= C_i,\,\, i=1, \dots, I, \,\, &\sum_{i=1}^I w_{ij} = \gamma_j, \,\, j =1, \dots, J\\ \sum_{j=1}^J x_{ij} &= D_i,\,\, i=1, \dots, I, \,\, &\sum_{i=1}^I x_{ij} = \delta_j, \,\, j =1, \dots, J\end{align*}

Using the mip solver for Python, I set up in a small Jupyter notebook, which is pasted below. The solver is very quick, at least for the type of problem sizes that one could hope to solve by pen and paper. I haven't fully tested the bounds of its perfomance with either $I$ or $J$ get big, but presumably they would still be better than solving by hand... Again, last chance to avert your eyes!!!!

In [1]:

import numpy as np
from mip import Model, xsum, INTEGER

def get_prime_factors(number, primes=[2, 3, 5, 7]):
    factors = ()
    for p in primes:
        k = 0
        while number % p == 0:
            k += 1
            number /= p
        factors += (k,)
    return factors

def solve_cross_product(column, row):
    n_rows = len(column)
    n_cols = len(row)
        
    crossProduct = Model()
    u = [[crossProduct.add_var('u({},{})'.format(i,j), var_type=INTEGER) for j in range(n_cols)] for i in range(n_rows)]
    v = [[crossProduct.add_var('v({},{})'.format(i,j), var_type=INTEGER) for j in range(n_cols)] for i in range(n_rows)]
    w = [[crossProduct.add_var('w({},{})'.format(i,j), var_type=INTEGER) for j in range(n_cols)] for i in range(n_rows)]
    x = [[crossProduct.add_var('x({},{})'.format(i,j), var_type=INTEGER) for j in range(n_cols)] for i in range(n_rows)]
    
    
    for i in range(n_rows):
        for j in range(n_cols):
            ## all values must be nonnegative integers
            crossProduct += u[i][j] >= 0
            crossProduct += v[i][j] >= 0
            crossProduct += w[i][j] >= 0
            crossProduct += x[i][j] >= 0
    
            ## all entries must yield single digit entries
            crossProduct += np.log(2.0) * u[i][j] + np.log(3.0) * v[i][j] + np.log(5.0) * w[i][j] + np.log(7.0) * x[i][j] <= np.log(9.0)
            
    ## each row must have the appropriate factors in order to multiply to the values in the input column
    for i, c_ in enumerate(column):
        factors = get_prime_factors(c_)
        crossProduct += xsum(u[i][j] for j in range(n_cols)) == factors[0]
        crossProduct += xsum(v[i][j] for j in range(n_cols)) == factors[1]
        crossProduct += xsum(w[i][j] for j in range(n_cols)) == factors[2]
        crossProduct += xsum(x[i][j] for j in range(n_cols)) == factors[3]
        
    ## each column must have the appropriate factors in order to multiply to the values in the input row
    for j, r_ in enumerate(row):
        factors = get_prime_factors(r_)
        crossProduct += xsum(u[i][j] for i in range(n_rows)) == factors[0]
        crossProduct += xsum(v[i][j] for i in range(n_rows)) == factors[1]
        crossProduct += xsum(w[i][j] for i in range(n_rows)) == factors[2]
        crossProduct += xsum(x[i][j] for i in range(n_rows)) == factors[3]
    
    status = crossProduct.optimize()
    if crossProduct.num_solutions:
        y = np.ones((n_rows, n_cols))
        for v in crossProduct.vars:
            name = v.name
            row_, col_ = v.name.split('(')[-1].split(')')[0].split(',')
            val = v.x
            y[int(row_),int(col_)] *= 2**val if name[0] == 'u' else 3**val if name[0] == 'v' else 5**val if name[0] == 'w' else 7**val
        print(y)

In [2]:

column = np.array([280, 168, 162, 360, 60, 256, 126])
row = np.array([183708, 245760, 117600])
solve_cross_product(column, row)

[[7. 8. 5.]
 [3. 8. 7.]
 [9. 6. 3.]
 [9. 8. 5.]
 [3. 5. 4.]
 [4. 8. 8.]
 [9. 2. 7.]]