Cleve Moler is the author of the first MATLAB, one of the founders of MathWorks, and is currently Chief Mathematician at the company. He is the author of two books about MATLAB that are available online. He writes here about MATLAB, scientific computing and interesting mathematics.
The algorithms used by MATLAB for generating magic squares of order n fall into three cases:
odd, n is odd.
doubly-even, n is divisible by 4.
singly-even, n is divisible by 2, but not by 4.
The best known algorithm for generating magic squares of odd order is de La Loubere's method. Simon de La Loubere was the French ambassador to Siam in the late 17th century. I sometimes refer to his method as the "nor'easter algorithm", after the winter storms that move northeasterly up the coast of New England. You can see why if you follow the integers sequentially through magic(9).
The integers from $1$ to $n^2$ are inserted along diagonals, starting in the middle of first row and heading in a northeasterly direction. When you go off an edge of the array, which you do at the very first step, continue from the opposite edge. When you bump into a cell that is already occupied, drop down one row and continue.
We used this algorithm in MATLAB for many years. Here is the code.
A = zeros(n,n);
i = 1;
j = (n+1)/2;
for k = 1:n^2
is = i;
js = j;
A(i,j) = k;
i = n - rem(n+1-i,n);
j = rem(j,n) + 1;
if A(i,j) ~= 0
i = rem(is,n) + 1;
j = js;
A big difficulty with this algorithm and resulting program is that it inserts the elements one at a time---it cannot be vectorized.
A New Algorithm
A few years ago we discovered an algorithm for generating the same magic squares of odd order as de La Loubere's method, but with just four MATLAB matrix operations.
[I,J] = ndgrid(1:n);
A = mod(I+J+(n-3)/2,n);
B = mod(I+2*J-2,n);
M = n*A + B + 1;
Let's see how this works with n = 5. The statement
[I,J] = ndgrid(1:n)
produces a pair of matrices whose elements are just the row and column indices, i and j.
Both A and B are fledgling magic squares. They have equal row, column and diagonal sums. But their elements are not the integers from $1$ to $n^2$. Each has duplicated elements between $0$ and $n-1$. The final statement
M = n*A+B+1
produces a matrix whose elements are integers between $1$ and $n^2$ and which has equal row, column and diagonal sums. What is not obvious, but is true, is that there are no duplicates. So M must contain all of the integers between $1$ and $n^2$ and consequently is a magic square.
The doubly-even algorithm is also short and sweet, and tricky.
M = reshape(1:n^2,n,n)';
[I,J] = ndgrid(1:n);
K = fix(mod(I,4)/2) == fix(mod(J,4)/2);
M(K) = n^2+1 - M(K);
Let's look at our friend magic(4). The matrix M is initially just the integers from 1 to 16 stored sequentially in 4 rows.
1 2 3 4
5 6 7 8
9 10 11 12
13 14 15 16
The logical array K is true for half of the indices and false for the other half in a pattern like this.
1 0 0 1
0 1 1 0
0 1 1 0
1 0 0 1
The elements where K is false, that is 0, are left alone.
. 2 3 .
5 . . 8
9 . . 12
. 14 15 .
The elements where K is true, that is 1, are reversed.
16 . . 13
. 11 10 .
. 7 6 .
4 . . 1
The final result merges these two matrices to produce the magic square.
Singly Even Order
The algorithm for singly even order is the most complicated and so we will give just a glimpse of how it works. If n is singly even, then $n/2$ is odd and magic(n) can be constructed from four copies of magic(n/2). For example, magic(10) is obtained from A = magic(5) by forming a block matrix.
[ A A+50
A+75 A+25 ]
The column sums are all equal because sum(A) + sum(A+75) equals sum(A+50) + sum(A+25). But the rows sums are not quite right. The algorithm must finish by doing a few swaps of pieces of rows to clean up the row sums. For the details, issue the command.
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.