Example: Magic Squares

The model will eliminate symmetries by imposing order constraints. Propagation will be drastically improved by exploiting a redundant constraint.

Problem Specification

The Magic Square Puzzle consists in finding for given N an N x N - matrix such that:

• Every field of the matrix is an integer between 1 and N².
• The fields of the matrix are pairwise distinct.
• The sums of the rows, columns, and the two main diagonals are all equal.

A magic square for N = 3 looks as follows:

2 7 6
9 5 1
4 3 8

The Magic Square Puzzle is extremely hard for large N. Even for N = 4, our script will have to explore more than 10000 nodes to find a solution.

Viewpoint and Constraints

We model the problem by having a variable F_ij for every field (i,j) of the matrix. Moreover, we have one additional variable S and require that the sum of every row, column, and main diagonal equals S.

Without loss of generality, we can impose the following order constraints eliminating symmetries:

F₁₁ < F_NN F_N1 < F_1N F₁₁ < F_N1

Since the sum of the sums of the rows must equal the sum of all fields, we have the redundant constraint

${\frac{{N^2}}{{2}}}$ * (N² + 1) = N * S

To see this, note that sum of all fields equals

1 + 2 + 3 + ... + N² = ${\frac{{N^2}}{{2}}}$ * (N² + 1)

and that the sum of each of the N rows must be S.

Branching Strategy

We branch on the variables F₁₁,..., F_NN with a first-fail strategy splitting the domain of the selected variable and trying the lower half first.

Script

val maxValue = (valOf (Int.maxInt))div 2

fun magicsquares n space =
    let
        val nn = n * n
        val matrix = Vector.tabulate(n,fn x => 
                              FD.rangeVec(space,n,(1,nn)))
        val matrixv =Vector.concat(Vector.toList(matrix))
        val sum = FD.range(space,(1,maxValue))
        val sumn = FD.range(space,(1,maxValue))
        val diagonal1 = Vector.tabulate(n,fn x =>(x,x))
        val diagonal2 = Vector.tabulate(n,fn x =>(n-1-x,x))
        val diagonal1vars = Vector.map(fn(x,y) => 
               Vector.sub(Vector.sub(matrix,x),y))diagonal1  
        val diagonal2vars = Vector.map(fn(x,y) => 
               Vector.sub(Vector.sub(matrix,x),y))diagonal2
        val field11 = Vector.sub(Vector.sub(matrix,0),0) 
        val field1N = Vector.sub(Vector.sub(matrix,0),n-1)
        val fieldN1 = Vector.sub(Vector.sub(matrix,n-1),0)
        val fieldNN = Vector.sub(Vector.sub(matrix,n-1),n-1)                 
    in
        FD.distinct(space,matrixv,FD.DOM);
        (* diagonals *)
        post(space,SUMV(diagonal1vars) `= FD(sum),FD.DOM);
        post(space,SUMV(diagonal2vars) `= FD(sum),FD.DOM);
        (* columns *)
        Vector.appi(fn(y,z)=>
         let 
            val colmn = Vector.tabulate(n,fn x =>
                         Vector.sub(Vector.sub(matrix,x),y)) 
         in
            post(space,SUMV(colmn) `= FD(sum),FD.DOM)
         end)matrix;
        (* rows *)
        Vector.appi(fn(x,y) =>
                     post(space,SUMV(y) `= FD(sum),FD.DOM))matrix;
        (* Eliminate symmetries *)
        FD.rel(space,field11,FD.LE,fieldNN);
        FD.rel(space,fieldN1,FD.LE,field1N);
        FD.rel(space,field11,FD.LE,fieldN1);
        (*redundant constraints *)
        post(space,FD(sum) `* `n `= FD(sumn),FD.DOM);
        FD.relI(space,sumn,FD.EQ,(nn *(nn + 1))div 2);
        
        FD.branch(space,matrixv,FD.B_SIZE_MIN,FD.B_SPLIT_MIN);
        {matrix}  
    end

With the Explorer you can find out that for N = 3 there is exactly one magic square satisfying the ordering constraints of our model. Without the ordering constraints there are 8 different solutions. Omitting the propagator for the redundant constraint will increase the search tree by one order of magnitude.