In [1] the following problem is described:

- The chessboard shown above has \(8\times 8 = 64\) cells.
- Place numbers 1 through 64 in the cells such that the next number is either directly left, right, below, or above the current number.
- The dark cells must be occupied by prime numbers. Note there are 18 dark cells, and 18 primes in the sequence \(1,\dots,64\). These prime numbers are: 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61.

This problem looks very much like the Numbrix puzzle [3]. We don't have the **given **cells that are provided in Numbrix, but instead, we have some cells that need to be covered by prime numbers.

The path formed by the consecutive values looks like how a rook moves on a chessboard. Which explains the title of this post (and of the original post). If we require a **proper tour** we need values 1 and 64 to be neighbors. Otherwise, if this is not required, I think a better term is: **Hamiltonian path**. I'll discuss both cases, starting with the Hamiltonian path problem.

#### Mathematical Model

As usual for this type of problems, we use the following decision variables: \[\color{darkred}x_{i,j,k} = \begin{cases} 1 & \text{if cell $(i,j)$ has a value $k$} \\ 0 & \text{otherwise}\end{cases}\] With this, our model can look like:

Hamiltonian Path Model |
---|

\[\begin{align} & \sum_k \color{darkred}x_{i,j,k} = 1 && \forall i,j &&\text{One value in a cell}\\ & \sum_{i,j} \color{darkred}x_{i,j,k} = 1 && \forall k &&\text{Unique values} \\ & \color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} && \forall i,j,k\lt 64 && \text{Next value is a neighbor} \\ & \sum_p \color{darkred}x_{i,j,p} = 1 && \text{for dark cells $(i,j)$} &&\text{Prime values} \\ & \color{darkred}x_{i,j,k} \in \{0,1\} \end{align} \] |

Here \(p\) is a subset of \(k\): all the prime values. The **next-value constraint** is explained in more detail in [3]. This is a fascinating constraint. The only thing we are saying here is: the next value is found among the neighbors of \((i,j)\). Or to be more precise: \[\color{darkred}x_{i,j,k}=1\Rightarrow \color{darkred}x_{i-1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i+1,j,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j-1,k+1}=1\>{\bf{or}}\>\color{darkred}x_{i,j+1,k+1}=1\]There is no concept of a path here. Rather, the formation of a path is a natural consequence. I assume in this constraint that when we address outside the board, the corresponding variable is zero. I think this model is quite elegant. It precisely states the conditions we must obey for a valid solution.

Notes:

- This model has no objective.
- We don't require a tour here. There is no restriction that says that values 1 and 64 must be neighbors. That version of the problem is looked at further down.
- This is not a very small model. The number of \(\color{darkred}x_{i,j,k}\) variables is \(8 \times 8 \times 64 = 4096\). The model size (after adding a dummy objective) is:

MODEL STATISTICS BLOCKS OF EQUATIONS 6 SINGLE EQUATIONS 4,179 BLOCKS OF VARIABLES 2 SINGLE VARIABLES 4,097 NON ZERO ELEMENTS 26,661 DISCRETE VARIABLES 4,096

A final question is: is the solution **unique **or are there multiple solutions? To forbid a previously found solution, we can add the cut: \[\sum_{i,j,k} \color{darkblue}x^*_{i,j,k} \cdot \color{darkred}x_{i,j,k} \le {\bf{card}}(k)-1\] where \(\color{darkblue}x^*_{i,j,k}\) is such a previously found solution. This constraint says: at least one cell should change. Note that we never can only change a 0 value to a 1: the total number of ones is given by \({\bf{card}}(k)=64\). So we only need to keep track of the ones. And that is what this constraint is doing: force to flip at least one 1 to a 0. We repeat solving and adding a cut to the model until the model becomes infeasible. An alternative approach would be to use the **solution pool** to multiple optimal solutions. The solution pool is a technology found in some advanced solvers. As the number of solutions is small, I use here the DIY approach using the "no good" cuts.

The feasible paths we find are:

---- 106 PARAMETERresultnegative numbers indicate prime cellsINDEX 1 = solution1 j1 j2 j3 j4 j5 j6 j7 j8 i1 -41 40 39 30 -29 28 21 20 i2 42 -37 38 -31 32 27 22 -19 i3 -43 36 35 34 33 26 -23 18 i4 44 45 46 -3 4 25 24 -17 i5 51 50 -47 -2 -5 8 9 16 i6 52 49 48 1 6 -7 10 15 i7 -53 56 57 58 -59 60 -11 14 i8 54 55 64 63 62 -61 12 -13 INDEX 1 = solution2 j1 j2 j3 j4 j5 j6 j7 j8 i1 -41 40 39 30 -29 28 21 20 i2 42 -37 38 -31 32 27 22 -19 i3 -43 36 35 34 33 26 -23 18 i4 44 45 46 -3 4 25 24 -17 i5 51 50 -47 -2 -5 8 9 16 i6 52 49 48 1 6 -7 10 15 i7 -53 64 63 62 -61 60 -11 14 i8 54 55 56 57 58 -59 12 -13 INDEX 1 = solution3 j1 j2 j3 j4 j5 j6 j7 j8 i1 -41 40 39 30 -29 28 21 20 i2 42 -37 38 -31 32 27 22 -19 i3 -43 36 35 34 33 26 -23 18 i4 44 45 46 -3 4 25 24 -17 i5 49 48 -47 -2 -5 8 9 16 i6 50 51 64 1 6 -7 10 15 i7 -53 52 63 62 -61 60 -11 14 i8 54 55 56 57 58 -59 12 -13 INDEX 1 = solution4 j1 j2 j3 j4 j5 j6 j7 j8 i1 -43 44 45 6 -7 8 9 10 i2 42 -47 46 -5 64 25 24 -11 i3 -41 48 49 4 63 26 -23 12 i4 40 51 50 -3 62 27 22 -13 i5 39 52 -53 -2 -61 28 21 14 i6 38 55 54 1 60 -29 20 15 i7 -37 56 57 58 -59 30 -19 16 i8 36 35 34 33 32 -31 18 -17

Or plotted using R/ggplot:

All feasible paths |

#### Finding feasible tours

Looking at these solutions, we see that just one solution forms a proper tour. The obvious question is: can we add a constraint to the model that just delivers this tour. Actually, this is rather simple. We need to change the equation next. The mathematical formulation looks like: \[\begin{cases}\color{darkred}x_{i-1,j,k+1}+\color{darkred}x_{i+1,j,k+1}+\color{darkred}x_{i,j-1,k+1}+\color{darkred}x_{i,j+1,k+1} \ge \color{darkred}x_{i,j,k} & \forall i,j,k\lt 64 \\ \color{darkred}x_{i-1,j,1}+\color{darkred}x_{i+1,j,1}+\color{darkred}x_{i,j-1,1}+\color{darkred}x_{i,j+1,1} \ge \color{darkred}x_{i,j,64} & \forall i,j\end{cases}\]

Here is how we can do this a bit more easily in GAMS. The original version had (see the appendix below for the complete model):

next(i,j,k+1)..

x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k);

next(i,j,k)..

x(i-1,j,k++1)+x(i+1,j,k++1)+x(i,j-1,k++1)+x(i,j+1,k++1) =g= x(i,j,k);

Here the ++ operator is used to implement a

**circular lead**. Also, note that the equation next is now indexed by all k. The net effect is that we add 64 rows to the model (one for each cell). If we run the model with this equation we get a single solution:

---- 106 PARAMETERresultnegative numbers indicate prime cellsINDEX 1 = solution1 j1 j2 j3 j4 j5 j6 j7 j8 i1 -41 40 39 30 -29 28 21 20 i2 42 -37 38 -31 32 27 22 -19 i3 -43 36 35 34 33 26 -23 18 i4 44 45 46 -3 4 25 24 -17 i5 49 48 -47 -2 -5 8 9 16 i6 50 51 64 1 6 -7 10 15 i7 -53 52 63 62 -61 60 -11 14 i8 54 55 56 57 58 -59 12 -13

#### Conclusions

**path**while obeying the prime-value-only cells, is a fairly close variant of the Numbrix model discussed in [3]. Finding a

**tour**can be handled by changing a single equation (but in an interesting way).

#### References

- Another Rook's Tour of the Chessboard, https://puzzling.stackexchange.com/questions/111818/another-rooks-tour-of-the-chessboard
- Chessboard, https://en.wikipedia.org/wiki/Chessboard
- Numbrix puzzle via MIP, https://yetanothermathprogrammingconsultant.blogspot.com/2020/12/numbrix-puzzle-via-mip.html
- Prime numbers in GAMS, https://www.gams.com/34/gamslib_ml/libhtml/gamslib_prime.html. Unfortunately, it skips the prime number 2.

#### Appendix: GAMS model

$ontext Place
the numbers 1..64 on an 8x8 board such that:1. next
values are neighbors. A neighbor is defined asbelow, left, right, or above.2.
certain cells need to contain a prime numberhttps://puzzling.stackexchange.com/questions/111818/another-rooks-tour-of-the-chessboard$offtext *----------------------------------------------------------------------* data*----------------------------------------------------------------------setsi 'rows' /i1*i8/j 'columns' /j1*j8/k 'values' /1*64/kprime(k) 'primes' /2,3,5,7,11,13,17,19,23,29, 31,37,41,43,47,53,59,61 / ; acronym p;table board(i,j) 'locations with a p should
get a prime value'j1 j2 j3 j4 j5 j6 j7 j8 i1 p p i2 p p p i3 p p i4 p p i5 p p p i6 p i7 p p p i8 p p ; *----------------------------------------------------------------------* model*----------------------------------------------------------------------set sol /solution1*solution10/;set solfound(sol) 'solutions found';parameter
solution(sol,i,j,k) 'solutions
found so far';solfound(sol) = no;solution(sol,i,j,k) = 0; binary variable x(i,j,k);variable dummy;equationsonevalue(i,j) 'each cell has one value
(between 1 and 81)'unique(k) 'each value must be in exactly one cell'next(i,j,k) 'x(i,j,k)=1 => a neighbor is k+1'isprime(i,j) 'we need a prime on (i,j)'cut(sol) 'cut off previous solutions'obj 'dummy'; obj.. dummy =e= 0; onevalue(i,j).. sum(k, x(i,j,k)) =e= 1;unique(k).. sum((i,j),x(i,j,k)) =e= 1;next(i,j,k+1).. x(i-1,j,k+1)+x(i+1,j,k+1)+x(i,j-1,k+1)+x(i,j+1,k+1) =g= x(i,j,k); isprime(i,j)$(board(i,j)=p).. sum(kprime, x(i,j,kprime)) =e= 1;cut(solfound).. sum((i,j,k),solution(solfound,i,j,k)*x(i,j,k)) =l= card(k)-1;model m /all/;option
mip=cplex,threads=8;*----------------------------------------------------------------------* solve loop*----------------------------------------------------------------------loop(sol,solve m minimizing dummy using mip;break$(m.modelstat <> %modelStat.optimal%and m.modelstat <>
%modelStat.integerSolution%);solfound(sol) = yes;solution(sol,i,j,k) = round(x.l(i,j,k)); display x.l;); *----------------------------------------------------------------------* reporting*----------------------------------------------------------------------parameter
result(sol,i,j) 'negative
numbers indicate prime cells';result(sol,i,j) = sum(k, k.val*solution(sol,i,j,k));result(sol,i,j)$(board(i,j)=p) = -result(sol,i,j); option
result:0:1:1;display result; |

- An acronym allows a string to be used in a numerical table (like inf, eps etc.).
- The equation next is indexed by \(k+1\). This means the equation is generated for all but the last k.
- If \(i-1\),\(i+1\),\(j-1\),\(j+1\) is outside the domain (i.e. outside the board), GAMS will return zero for the variable (or in other words: ignore that entry). That is exactly what we need here. So we don't need to add additional protection to remain inside the board.
- This version of the model generates all feasible Hamiltonian paths.
- We could have reduced the number of variables a bit by only generating valid \(\color{darkred}x_{i,j,k}\): for prime cells, non-prime \(k\)'s are not allowed. That would also eliminate the isprime equation.

https://gist.github.com/zayenz/46fb486a4454aac72127376436b823ef is a MiniZinc model of the proper tour version of the problem that uses the global circuit constraint to model the tour-part.

ReplyDeleteHere is the output from finding all the solutions using MiniZinc 2.5.4 and Gecode 6.3.0 as the solver

```

Running prime-position-tour.mzn

Board:

[| 41, 40, 39, 30, 29, 28, 21, 20 |

42, 37, 38, 31, 32, 27, 22, 19 |

43, 36, 35, 34, 33, 26, 23, 18 |

44, 45, 46, 3, 4, 25, 24, 17 |

49, 48, 47, 2, 5, 8, 9, 16 |

50, 51, 64, 1, 6, 7, 10, 15 |

53, 52, 63, 62, 61, 60, 11, 14 |

54, 55, 56, 57, 58, 59, 12, 13 |]

Next:

[| 9, 1, 2, 12, 4, 5, 15, 7 |

17, 11, 3, 13, 21, 6, 23, 8 |

25, 10, 18, 19, 20, 14, 31, 16 |

26, 27, 35, 29, 37, 22, 30, 24 |

41, 33, 34, 28, 45, 39, 47, 32 |

42, 50, 44, 36, 46, 38, 55, 40 |

57, 49, 43, 51, 52, 53, 63, 48 |

58, 59, 60, 61, 62, 54, 64, 56 |]

----------

==========

%%%mzn-stat: initTime=0.003434

%%%mzn-stat: solveTime=0.29611

%%%mzn-stat: solutions=1

%%%mzn-stat: variables=256

%%%mzn-stat: propagators=236

%%%mzn-stat: propagations=1076445

%%%mzn-stat: nodes=3713

%%%mzn-stat: failures=1856

%%%mzn-stat: restarts=0

%%%mzn-stat: peakDepth=22

%%%mzn-stat-end

Finished in 656msec

```

That is very fast. Of course the model is smaller (say in terms of number of variables) and you are using a special purpose constraint (circuit). It demonstrates that a good CP model is quite different than a MIP model.

DeleteHow long did it take for your model to solve the problem? This is probably a case that really benefits from having the right global constraint available.

DeleteThe model uses a mod-constraint, which is not the best in terms of propagation. It might be beneficial to replace it with an extensional constraint using a table, just listing the valid pairs of values ([<1,2>,<2,3>,...,<64,1>]) instead. It was fast enough that it didn't matter here.

Note that it is very important to run the search on the next-variables in the model. If the board-variables (with the numbers for each position) are used for searhc, the problem is much harder. I would guess that an even faster solution would be to use Warnsdorffs heuristic for the search (follow the current path in the next-variables, and go to the next possibility with the fewest remaining possibilities), since that works well for Kinghts tour problems. For that one would have to implement it in a solver like Gecode. Might not work for this case, but could be interesting to try some time.

Cplex took 6 seconds to find the tour (default settings on a slow laptop).

Delete