The magic square function

Question

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?

Answer 1

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]
   ];

Answer 2

@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)
    ] // Experimental`CompileEvaluate;

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