Magic squares are cool. Algorithms that create magic squares are even cooler.

You probably remember magic squares from your childhood: they are *n* x *n* matrices that contain the numbers 1,2,...,n^{2} and for which the row sum, column sum, and the sum of both diagonals are the same value.
There are many algorithms to generate
magic squares. I remember learning as a child how to construct a magic square for any odd number, *n*, by using the "Siamese method." I also remember being fascinated by Ben Franklin's construction of "semi-magic squares" in which he used the sum of "bent diagonals" instead of straight diagonals.

Last week I had email discussions with several readers who had the idea to generate magic squares in SAS. The readers noted that my recent blog post on "How to return multiple values from a SAS/IML function," which features a module that computes row sums and column sums of a matrix, could be extended to check whether a matrix is a magic square.

When I first came to SAS, I was disappointed because there is no function in the SAS/IML language that generates magic squares. Prior to arriving at SAS, I had used MATLAB for matrix computations, and the MATLAB documentation is full of calls to the `magic` function, which generates magic squares of any size. An internet search for "magic squares in SAS" did not yield any useful results. However, when I searched for "MATLAB magic squares," I was thrilled to discover that Cleve Moler—cofounder of The MathWorks and numerical analyst extraordinaire—has written an e-book (C. Moler, 2011, *Experiments with MATLAB*)
that includes a chapter on the construction of magic squares!

The following program implements Moler's algorithm in the SAS/IML language. For the description of the algorithm, see Chapter 10 (p. 7–10) of Moler's e-book.

proc iml; start Magic(n); /* Moler's algorithm for constructing a magic square for any n. See Moler, C. (2011) "Magic Square," Chapter 10 in Experiments with MATLAB, p. 7-10 E-Book: http://www.mathworks.com/moler/exm/chapters/magic.pdf */ /* define helper functions */ start ModPos(a, q); /* always return positive value for mod(a,q) */ return( a - q*floor(a/q) ); finish; /* Magic square when n is odd */ start MagicOdd(n); I = repeat( T(1:n), 1, n ); J = repeat( 1:n, n ); /* or T(I) */ A = ModPos(I+J-(n+3)/2, n); /* modify formula to match MATLAB output */ B = mod(I+2*J-2, n); return( n*A + B + 1 ); finish; /* Magic square when n is a multiple of 4 */ start Magic4(n); M = shape(1:n##2, n, n); I = repeat( T(1:n), 1, n ); J = repeat( 1:n, n ); /* or T(I) */ idx = loc(floor(mod(I,4)/2) = floor(mod(J,4)/2)); M[idx] = n##2 + 1 - M[idx]; return( M ); finish; /* Magic square when n>4 is even but not a multiple of 4 */ start Magic2(n); s = n/2; /* s is odd */ A = MagicOdd(s); M = (A || (A + 2*s##2)) // /* 2x2 blocks of magic squares */ ((A+3*s##2) || (A + s##2)); /* col sums are correct; permute within cols to adjust row sums */ i = 1:s; r = (n-2)/4; h = do(n-r+2,n,1); /* might be empty */ j = (1:r) || h; M[i || i+s, j] = M[i+s || i, j]; /* swap rows for cols in j */ i = r+1; j = 1 || i; M[i || i+s, j] = M[i+s || i, j]; /* swap rows for cols in j */ return( M ); finish; /* main algorithm */ if n<=1 then return(1); if n=2 then return( {. ., . .} ); /* no magic square exists */ if mod(n,2)=1 then return( MagicOdd(n) ); /* n is odd */ if mod(n,4)=0 then return( Magic4(n) ); /* n is divisible by 4 */ return( Magic2(n) ); /* n is divisible by 2, but not by 4 */ finish;

The implementation uses two interesting features. First, I define a helper function that always returns a positive value for the expression mod(a,q), because the MOD function in SAS can sometimes return a negative value. Secondly, the computation for the Magic2 function constructs a block matrix of size 2*n* from a magic square of size *n*. Ian Wakeling told me that this is an application of a general technique to construct a large magic square from smaller components.

To test the implementation, you can construct and print a few magic squares, as follows:

do n = 3 to 6; M = Magic(n); L = cat("Magic Square, n=", char(n,1)); /* form label */ print M[label=L]; end;

At long last, I have my wish: a function in PROC IML that generates magic squares!

Do you have a favorite memory regarding magic squares? Do you have a favorite algorithm that constructs a magic square? Post a comment!