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

Question

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 =
   MapThread[
    Equal,
    {Floor[Mod[mat2, 4]/2], Floor[Mod[mat1, 4]/2]}, 2];
  Apply[
   If[#2, n^2 + 1 - #1, #1] &,
  MapThread[List, {square, bools}, 2], {2}]
 ]

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.

Answer 1

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 = Compile`GetElement}, 
   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 = Compile`GetElement}, 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.