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,...,n2 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 2n 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!
Update: As of SAS/IML 12.3, the MAGIC function is built into SAS.
Do you have a favorite memory regarding magic squares? Do you have a favorite algorithm that constructs a magic square? Post a comment!
12 Comments
Very nice indeed Rick, thank you for sharing your code with us. It does bring back old memories.
I would like to mention another method for constructing magic squares which is based on obtaining Latin squares from PROC FACTEX (SAS/QC) with additional code in IML. For those who would like to try, I have posted my program to the IML Discussion Forum here.
https://communities.sas.com/thread/37262
The IML part of the code is based on constructing larger squares from smaller ones, the first and coolest, is a simple multiplication of two squares of size m and n, to give a square of size m x n. The direct product operator '@' available in IML makes this trivially easy:
start MagicProduct(X,Y);
m=nrow(X);
n=nrow(Y);
return( j(n,n,1)@X + (Y-1)@j(m,m,m#m) );
finish;
There is also a doubling method which I believe is essentially equivalent to your module called Magic2(n) above (for the method see this web page http://www.grogono.com/magic/10x10.php). I have not done extensive tests, but you could try using the following code instead:
start Magic2(n);
m=n/2;
v=0 || j(1, m, 1) || j(1, m-1, 0);
v=repeat( v//(1-v), m, 1); /* basic pattern */
v[(n-1):n, ]=v[(n-1):n, n:1]; /* reverse last two rows */
return( (MagicOdd(m)-1)@j(2,2,4) + 2#v + t(v)[,n:1] +1);
finish;
In my program this module is called Magic2Q, and the last line looks slightly different as my program makes squares where the smallest element is a zero rather than one.
I have found a rearrangement of the formula for the MagicOdd construction that will guarantee that the first argument of the two mod functions can never go negative, so in this instance the helper function is not necessary.
start MagicOdd(n);
I = repeat( T(1:n), 1, n );
J = repeat( ((n-1)/2) : (3#(n-1)/2), n );
A = mod(I + J, n);
B = mod(I + 2#J + 1, n);
return( n*A + B + 1);
finish;
Thanks, Ian! I programmed it the way I did because that's the way that it was constructed in Moler's book. In MATLAB mod(a,q) is always positive when q>0, so I wrote the helper funtion to mimic the MATLAB mod() function. But you're right: in practice there is no need for the helper function.
Very cool Rick....I bet you can solve my "evil" level Sudoku's with your magic too!
This post is on the sas.com main page! It should really increase the ranking of you page in google:) I think every language that uses operations with matrices should have this example. I thought that magic square was really cool when I started learning matlab.
The wiki on Associative Magic Squares provides two references for recent programs that create magic squares. These programs utilize Local Search and Tabu search. The program on Sourceforge searches for the maximum water retention on magic squares. The Complete utility on the download section of Harry White's website allows one to specify certain values in a magic square and it will complete the skeleton and make it magic. Check out the 14 x 14 magic square with water writing in the square at the end of the wiki noted above.
Thanks for all the coments above !
Craig
Hi Rick,
Here's a rather unexpected way of making new magic squares from old.
A ring of magic powers!
x = magic(4);
print x [f=2.0];
do cycle = 1 to 8;
x = mod( x**5, 17);
print x [f=2.0];
end;
The original square, and the one from cycle 1, are rotated anticlockwise in steps of 90 degrees.
Can I make an animated heatmap to visualise something like this?
Amazing. You'll have to explain why multiplying the matrix by itself 5 times (why 5?) mod 17 results in another magic square. Very cool.
Speaking of rotating matrices, check out this blog post: http://blogs.sas.com/content/iml/2013/10/18/rotating-matrices/
In SAS/IML 13.1 you can create heat maps to visualize matrices with the HEATMAPCONT and HEATMAPDISC subroutines. I need to think about how to animate the results.
I am not totally sure, other than to say that is to do with the pattern of eigenvalues. I read "Odd Magic Powers", A. C. Thompson, The American Mathematical Monthly, 101(1994), pp. 339-342, and then started experimenting in IML. Only odd powers will work and it seems that squares of order 4n produced by your function are specially symmetrical (but not pan-diagonal).
.
I have an even more fantastic example. This produces a set of 6 magic squares!
.
x = magic(8);
print x [f=3.0];
do cycle=1 to 6;
x = rank( mod( x**5, 346) );
print x [f=3.0];
end;
.
If you try larger squares or larger powers, be careful not to exceed the largest integer that can be exactly represented.
Dr Rick,
In your program that implements Moler's algorithm in the SAS/IML language, the computation for the Magic2 function constructs a block matrix of size 2n from a magic square of size n, as you said it. Don’t you think that the amount of time needed when n is big (n > 100) can be reduced if a well-defined function f(i, j) is stored for any variable n? Don’t you think that this renders the generation of magic squares more straight?
Contact me if necessary, I mean if you are in search of such a function.
Best wishes
Bay Ima
Pingback: Overview of new features in SAS/IML 12.3 - The DO Loop
Latin Square of the Moon can generate a Magic Square. If we replace the numbers in the first column of “moon” with their respective root numbers we commence filling a root Square. Read more about it on my blog.