5
$\begingroup$

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

enter image description here enter image description here

For example, let $n=4$

enter image description here enter image description here

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.
$\endgroup$
  • $\begingroup$ Where's the definition of magic1? By the way Tr@Transpose@mat should be something like Tr@Reverse@mat. $\endgroup$ – xzczd May 6 '15 at 11:09
7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thanks a lot. And I need some time to understand your code:) $\endgroup$ – xyz May 7 '15 at 0:32
  • $\begingroup$ @ShutaoTang I've added some brief explanation, hope it helps. $\endgroup$ – xzczd May 7 '15 at 3:29
  • $\begingroup$ 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} $\endgroup$ – xyz May 7 '15 at 6:24
  • $\begingroup$ BTW, How can I find the documentation about CompileGetElement`? $\endgroup$ – xyz May 7 '15 at 7:05
  • $\begingroup$ @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? $\endgroup$ – xzczd May 7 '15 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.