The magic square function
@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
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]
];