4
$\begingroup$

I think there is exactly one magic square of order 3 (though I could be wrong about that). When I run this, it does not come up with any solutions. Can you help? The following code is meant to return all order 3 magic squares.

n = 3;
c = (n (n^2 + 1))/2;
equations = {
   Element[a11, Integers],
   Element[a21, Integers],
   Element[a31, Integers],
   Element[a12, Integers],
   Element[a22, Integers],
   Element[a32, Integers],
   Element[a13, Integers],
   Element[a23, Integers],
   Element[a33, Integers],
   1 <= a11 <= n,
   1 <= a21 <= n,
   1 <= a31 <= n,
   1 <= a12 <= n,
   1 <= a22 <= n,
   1 <= a13 <= n,
   1 <= a23 <= n,
   1 <= a33 <= n,
   a11 + a21 + a31 == c,
   a12 + a22 + a32 == c,
   a13 + a23 + a33 == c,
   a11 + a12 + a13 == c,
   a21 + a22 + a23 == c,
   a11 + a22 + a33 == c,
   a13 + a22 + a31 == c,
   a31 + a32 + a33 == c
 };

Reduce[equations, {a11, a21, a31, a12, a22, a32, a13, a23, a33}]
Improve this question
$\endgroup$
1
  • $\begingroup$ You have imposed the criteria that all numbers should be less or equal to three. For a traditional magic square, this is not a requirement, however. It also makes it impossible to find a solution. The sum of a row of three element cannot be fifteen if each element must be at most three. $\endgroup$
    – C. E.
    Jun 13, 2020 at 23:51

1 Answer 1

Reset to default
8
$\begingroup$

There's a much more compact way to represent this problem:

magicSquareConstraints[n_, k_, c_] := 
 Module[{sq = Table[a[i, j], {i, n}, {j, n}]},
  Join[
   (1 <= # <= k) & /@ Flatten[sq],
   (Total[#] == c) & /@ sq,
   (Total[#] == c) & /@ Transpose[sq],
   {
    Total[#] == c &@Diagonal[sq],
    Total[#] == c &@Diagonal[Reverse /@ sq]
   }
   (* uncomment this if you want unique entries 
    , Not@*Equal @@@ Subsets[Flatten@sq, {2}]
   *)
  ]]

With[{n = 3, k = 9, c = 15, s = 2},
 mtx = Table[a[i, j], {i, n}, {j, n}];
 mtx /. FindInstance[magicSquareConstraints[n, k, c],
   Flatten[mtx], Integers, s]
 ]

(* results: {{{3, 7, 5}, {7, 5, 3}, {5, 3, 7}},
            {{8, 1, 6}, {3, 5, 7}, {4, 9, 2}}} *)

n is the dimensions of the matrix, k is the maximum number (from 1 to k) to put in the boxes, c is the desired total along any row, column or diagonal, and s is the desired number of solutions to find. In the above s=2 and I have two solutions that are magic squares.

Also uncomment the constraint Not@*Equal @@@ Subsets[Flatten[sq], {2}] if you want all entries to be unique numbers. For example, with {n = 4, k = 16, c = 34, s = 1} and the unique constraint, it took about 10 minutes to find this $4\times4$ magic square: $$ \left( \begin{array}{cccc} 1 & 2 & 15 & 16 \\ 13 & 14 & 3 & 4 \\ 12 & 7 & 10 & 5 \\ 8 & 11 & 6 & 9 \\ \end{array} \right) $$

Update: There also appears to be a magic square ResourceFunction on the function repository. It's very fast though it only handles odd dimensions, and it only seems to give me one solution:

msq = ResourceFunction["MagicSquare"]
msq[5] // MatrixForm

I have also been looking into a linear programming solution after reading this on Wolfram blog. LinearOptimization is a lot faster in the 4x4 case than using FindInstance but Mathematica crashes when I add in the diagonals constraints and I don't know why and I have now 'fixed' the crash by moving my diagonal constraints:

n = 4; c = 34;
ones = ConstantArray[1, n^2];
vars = Table[v[i, j], {i, n}, {j, n}];
fvars = Flatten[vars];
auxvars = Map[Range[n^2].# &, vars, {2}];
diag1 = Diagonal[auxvars];
diag2 = Diagonal[Reverse@auxvars];
constraints = Join[
   (* sum of the rows == c*)
   Total[#] == c & /@ auxvars,
   (* sum of the columns == c *)
   Total[#] == c & /@ Transpose[auxvars],
   (* sum of the diagonals == c *)
   Total[#] == c & /@ {diag1},
   Total[#] == c & /@ {diag2},
   {
    (* all numbers used *)
    Total[fvars] == ConstantArray[1, n^2],
    (* each variable has a single 1 *)
    Total[#] == 1 & /@ fvars,
    (* vector integer constraints for linear programming *)
    0 \[VectorLessEqual] # \[VectorLessEqual] 1 & /@ fvars,
    # \[Element] Vectors[n^2, Integers] & /@ fvars
    }];
Map[FirstPosition[#, 1][[1]] &,
 vars /. LinearOptimization[0, constraints, Flatten[vars]], {2}] // MatrixForm

$$ \left( \begin{array}{cccc} 1 & 12 & 15 & 6 \\ 7 & 9 & 4 & 14 \\ 10 & 8 & 13 & 3 \\ 16 & 5 & 2 & 11 \\ \end{array} \right) $$ Finally, as I'm becoming a bit obsessed with this answer, I'll just leave you with the following weird property I noticed with these two $4\times4$ magic squares: if you apply Mod[# + 7, 16] + 1 to all elements of a $4\times4$ square with sums of $34$, you get another magic square with permuted elements except the diagonals sum to $18$ and $50$.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Nice! Thanks... $\endgroup$
    – Veritas Lux
    Jun 16, 2020 at 16:51
  • $\begingroup$ @VeritasLux Please accept if you feel this is a satisfactory answer, unless there is something else I can help with. $\endgroup$
    – flinty
    Jun 16, 2020 at 16:53
  • $\begingroup$ for some reason, it doesn't give me the 4x4 magic square. I'm think it is my computer or perhaps how mathematica is configured, such as some depth bailout measure by which it kills the search for a 4x4 magic square. It works for about 20 minutes and looked like something ended the kernel with no magic square. Also, your code doesn't seem to want to work for 2x2 magic squares. I did get all 8 3x3 magic squares! Any idea why this behavior is? I have 32GB of ram but pretty sure mathematica is using a fraction of that. TIA $\endgroup$
    – Veritas Lux
    Jun 17, 2020 at 19:52
  • $\begingroup$ That's a shame about 4x4 - it took a while on my machine too but it did complete. As for 2x2 - the smallest magic square is 3x3 and 2x2 is impossible. If you remove both the diagonal constraints and the unique constraint then 2x2 works and you'll get solutions with a distinct number on each diagonal, or all four entries the same. $\endgroup$
    – flinty
    Jun 17, 2020 at 20:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.