# Speed up the auxiliary function magicSquare when $n$ is doubly-even

Today, I saw a question about magic square (please see here). When $n$ is doubly-even, namely for Mod[n,4]==0

magicSquare[n_] /; Mod[n, 4] == 0 :=
Module[{mat1, mat2, square, bools},
mat1 = Table[Range[n], {i, n}];
mat2 = Transpose@mat1;
square = Partition[Range[1, n^2], n];
bools =
Equal,
{Floor[Mod[mat2, 4]/2], Floor[Mod[mat1, 4]/2]}, 2];
Apply[
If[#2, n^2 + 1 - #1, #1] &,
]


Algorithm for doubly-even

For example, let $n=4$

Test

magicSquare[3000]; // AbsoluteTiming (* n is doubly-even*)

{40.0029297, Null}


Obviously, it is rather time-consuming. In addition, I write a function magicQ[A] that checks if A is a magic square

magicQ[mat_?MatrixQ] /; Equal @@ Dimensions@mat :=
Module[
{n = Length@mat, row, col, diagonal},
row = Total[mat, {2}];
col = Total[mat, {1}];
diagonal = {Tr@mat, Tr@Transpose@mat};
Union[row, col, diagonal] == {(n^3 + n)/2}
]

magicSquare[3000] // magicQ // AbsoluteTiming

{6.8847656, True}


Question

• How to speed up this magicSquare function?
• Search for more efficient solution for magicQ.
• Where's the definition of magic1? By the way Tr@Transpose@mat should be something like Tr@Reverse@mat. – xzczd May 6 '15 at 11:09

The main reason for the slowness of the magicSquare is that your failed to insist on vecterization. (you're already aware of its importance, right? ) Making use of the internal functions owning Listable attribute and treating lists as a whole as much as possible, it's not hard to come to the following:

magicSquare2[n_] /; Mod[n, 4] == 0 :=
Module[{mat, square, bools}, mat = ConstantArray[Range[1., n], n];
square = Partition[Range[n^2], n];
bools = UnitStep@-Abs[# - #\[Transpose] &@Floor[Mod[mat, 4]/2]];
n^2 + 1 - #1 &@square bools + (1 - bools) square]

ans1 = magicSquare[3000]; // AbsoluteTiming
ans2 = magicSquare2[3000]; // AbsoluteTiming

{23.320632, Null}
{0.796866, Null}


Two other tiny changes I've made in magicSquare2 are

1. Use ConstantArray instead of Table.
2. Change Range[n] to Range[1., n] so mat will become a MachinePrecision matrix and the calculation of Mod[mat, 4]/2 will be sped up.

If you have a C compiler installed:

magicSquare3[n_] /; Mod[n, 4] == 0 := cfmagicSquare3[n]
cfmagicSquare3 =
With[{g = CompileGetElement},
Compile[{{n, _Integer}},
Module[{mat1, mat2, square},
mat1 = Floor[Mod[Table[Range[1., n], {n}], 4]/2];
mat2 = mat1\[Transpose];
square = Partition[Range[n^2], n];
Table[
If[g[mat1, i, j] == g[mat2, i, j], n^2 + 1 - g[square, i, j],
g[square, i, j]], {i, n}, {j, n}]], CompilationTarget -> C,
RuntimeOptions -> "Speed"]];

ans3 = magicSquare3[3000]; // AbsoluteTiming

{0.281241, Null}


Tricks used in the above piece of code are almost all inherited from this post. The only thing I can add here is, mat2 = mat1\[Transpose]; is necessary for speed, I'm not sure about the reason, maybe it's because that if one eliminate this line and use g[mat1, i, j] == g[mat1, j, i] in the last Table, g[mat1, j, i] will cause more frequent access for the higher dimension of mat1, which is slower than that for lower dimension. (You remember the truth that for n = 10^3; test = RandomReal[1, {n, n}];, test[[All, 1]] is slower than test[[1]]? )

As to the magicQ, despite of the mistake mentioned in the comment above, I think it's already quite efficient and simple. Of course one can create a faster one, still, by compiling to C:

With[{g = CompileGetElement}, Block[{check1, check2},
check1[mat_, n_, tot_]:= Do[If[Sum[g[mat, i, j], {j, n}] != tot, Throw@False], {i, n}];
check2[mat_, n_, tot_]:= If[Sum[g[mat, j, j], {j, n}] != tot, Throw@False];
cfmagicQ2 =
Hold@Compile[{{mat, _Integer, 2}},
Module[{n = Length@mat, mat2 = mat\[Transpose], mat3 = Reverse[mat], tot},
tot = (n^3 + n)/2;
Catch[check1[mat, n, tot]; check1[mat2, n, tot]; check2[mat, n, tot];
check2[mat3, n, tot]; True]], CompilationTarget -> C,
RuntimeOptions -> "Speed"] //. Flatten[DownValues /@ {check1, check2}] //
ReleaseHold];
magicQ2[mat_?MatrixQ] /; Equal @@ Dimensions@mat := cfmagicQ2[mat]];

test = ans3; test[[100, 100]] = 0;
ans3 // magicQ // AbsoluteTiming
test // magicQ // AbsoluteTiming
ans3 // magicQ2 // AbsoluteTiming
test // magicQ2 // AbsoluteTiming

{0.156215, True}
{0.156223, False}
{0.062461, True}
{0.046865, False}


Tricks for speeding up the compiled code is the same as those used in magicSquare3. Throw and Catch will make the calculation stop once the matrix is detected as "non-magic" thus judgement for non-magic square is sped up. Notice that I've used a code auto-generation technique for simplicity, which is inherited from this post.

• Thanks a lot. And I need some time to understand your code:) – xyz May 7 '15 at 0:32
• @ShutaoTang I've added some brief explanation, hope it helps. – xzczd May 7 '15 at 3:29
• UnitStep@-Abs[# - #\[Transpose] &@Floor[Mod[mat, 4]/2]] is much efficient than MapThread[ Equal, {Floor[Mod[mat2, 4]/2], Floor[Mod[mat1, 4]/2]}, 2. For instance, mat1 = RandomReal[{1, 10}, {3000, 3000}];MapThread[ Equal, {Floor[Mod[mat1\[Transpose], 4]/2], Floor[Mod[mat1, 4]/2]}, 2]; // AbsoluteTiming UnitStep@-Abs[# - #\[Transpose] &@ Floor[Mod[mat1, 4]/2]]; // AbsoluteTiming gives the following result {5.6494141, Null},{0.7744141, Null} – xyz May 7 '15 at 6:24
• BTW, How can I find the documentation about CompileGetElement`? – xyz May 7 '15 at 7:05
• @ShutaoTang It's undocumented, but mentioned in a few posts here, for example this. You can have a search for more…Well, I admit I found no detailed explanation for so long. Maybe you can consider asking a separate question for this issue? – xzczd May 7 '15 at 7:19