# 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 '15 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/… – Mr.Wizard Feb 3 '15 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 '15 at 16:21
• @Hans thanks. I changed the signature to check for n!=2 – Nasser Feb 3 '15 at 17:52
• Thank you very much. The output is also beautiful. – Nick Feb 6 '15 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[3000]; // AbsoluteTiming
magic[3001]; // AbsoluteTiming
magic[3002]; // 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. – Nasser Feb 3 '15 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 '15 at 12:45
• @Nick Don't worry about that. – chyanog Feb 7 '15 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 '15 at 8:34