# The magic square function

There are some matrices such that the sum of columns, the sum of rows, and the sum of diagonals are the same value. Here is an example:

 8     1     6
3     5     7
4     9     2


The sum of rows is 15, and so is sum of columns and the sum of diagonals.

In MATLAB, refer to Cleve Moler's book "Experiments with MATLAB". I can generate this kind of matrix using the magicfunction, with an argument specifying the size of the matrix:

4-by-4:

>> magic(4)
ans =
16     2     3    13
5    11    10     8
9     7     6    12
4    14    15     1


5-by-5:

>> magic(5)
ans =
17    24     1     8    15
23     5     7    14    16
4     6    13    20    22
10    12    19    21     3
11    18    25     2     9


Is there a similar function in Mathematica? Or, maybe there's some way to build this kind of matrix?

• @Amzoti No, I haven't. Thank you for pointing it out, great source!
– Nick
Feb 3, 2015 at 7:33
• @Amzoti Unfortunately that Notebook relies on a package that is no longer available as far as I can tell. Here is an image of the download page in the Internet Archive, but the MagicSquares link does not work: web.archive.org/web/20110215143835/http://library.wolfram.com/… Feb 3, 2015 at 7:51

Translated Cleve Moler's magic() function from Matlab code to Mathematica.

Grid[Partition[MatrixForm@magic[#] & /@ {3, 4, 5, 6, 7, 8, 9, 10}, 4],
Frame -> All, FrameStyle -> LightGray] code:

magic[n_Integer /; (n > 0 && n != 2)] := Module[{m, j, k, p, i},
(*Translation of Cleve Moler's magic magic() function to Mathematica*)
Which[
Mod[n, 2] == 1, m = oddOrderMagicSquare[n],
Mod[n, 4] == 0,
j = Floor  @ Abs [ Mod[Range[n], 4]/2];
k = Outer[Equal, j, j] /. {True -> 1, False -> 0};
m = Outer[Plus, Range[1, n*n, n], Range[0, n - 1]];
p = Position[k, 1];
(m[[Sequence @@ #]] = n*n + 1 - m[[Sequence @@ #]]) & /@ p,
True,
p = n/2;
m = oddOrderMagicSquare[p];
m = ArrayFlatten@{{m, m + 2*p^2}, {m + 3*p^2, m + p^2}};
If[n != 2,
i = Range[p];
k = (n - 2)/4;
j = {Range[k], Range[n - k + 2, n]};
j = Flatten@DeleteCases[j, {}];
m[[Join[i, i + p], j]] = m[[Join[i + p, i], j]]
]
];
m
];
oddOrderMagicSquare[n_] := Module[{p},
p = Range[n];
Transpose[n*Mod[Map[p + # &, p - (n + 3)/2], n] +
Mod[Map[p + # &, 2*p - 2], n] + 1]
];

• Put some conditions on n == 2. No magic square of order 2.
– Hans
Feb 3, 2015 at 16:21
• @Hans thanks. I changed the signature to check for n!=2 Feb 3, 2015 at 17:52
• Thank you very much. The output is also beautiful.
– Nick
Feb 6, 2015 at 12:41

@Nasser's answer is nice, but slowly when Mod[n,4]==0. Here is a faster code, efficiency is close to Matlab :

ClearAll[magic]
magic[n_?OddQ] := oddOrderMagicSquare[n];

magic[n_ /; n~Mod~4 == 0] :=
Module[{J, K1, M},
J = Floor[(Range[n]~Mod~4)/2.0];
K1 = Abs@Outer[Plus, J, -J]~BitXor~1;
M = Outer[Plus, Range[1, n^2, n], Range[0, n - 1]];
M + K1 (n*n + 1 - 2 M)
] // ExperimentalCompileEvaluate;

magic[n_?EvenQ] :=
Module[{p, M, i, j, k},
p = n/2;(*p is odd*)
M = oddOrderMagicSquare[p];
M = ArrayFlatten@{{M, M + 2 p^2}, {M + 3 p^2, M + p^2}};
If[n == 2, Return[M]];
i = Transpose@{Range@p};
k = (n - 2)/4;
j = Range[k]~Join~Range[n - k + 2, n];
M[[Flatten@{i, i + p}, j]] = M[[Flatten@{i + p, i}, j]];
i = k + 1;
j = {1, i};
M[[Flatten@{i, i + p}, j]] = M[[Flatten@{i + p, i}, j]];
M
];

oddOrderMagicSquare[n_?OddQ] :=
Module[{p},
p = Range[n];
Outer[Plus, p, p - (n + 3)/2]~Mod~n*n +
Outer[Plus, p, 2 p - 2]~Mod~n + 1
];

magic; // AbsoluteTiming
magic; // AbsoluteTiming
magic; // AbsoluteTiming
`
• good you tried to speed it up. I translated Matlab code line by line, did not try to make any changes since I wanted first to get it working correctly. But speeding it up is interesting exercise on its own. Feb 3, 2015 at 11:11
• Thank you very much. It's really hard for me as a beginner to decide whose answer is the right one. Since @Nasser is the original answer, I marked his right answer. I'm sorry, I can only vote up for yours.
– Nick
Feb 6, 2015 at 12:45
• @Nick Don't worry about that. Feb 7, 2015 at 1:34
• @chyaong, Could you give some explanation(or your algorithm ) to help me to understand your high efficiency implementation when $n$ is doubly-even? thanks a lot:)
– xyz
May 6, 2015 at 8:34