A Computer Program to Find Magic Squares

I share an office with Charlie Gilderdale, author of "Magic squares for special occasions". His article describes an algorithm (method) for creating a 4 by 4 magic square with a particular date across the top. One day in February, Charlie was trying to use this algorithm to produce a magic square with top row 10 2 20 01, to use at a workshop on that date. He had no problem producing one, but he was being fussy: he wanted his square to contain no negative numbers, and have no number used more than once. With some dates, this is quite easy to do, by experimenting a little at the points where the algorithm gives a choice. However for this particular date, I spent half the morning running systematically through hundreds of possibilities before I found one of the four with no repeats.

One of the good things about a neatly described algorithm is that a computer can be taught to follow it. This algorithm involves a choice at four points: a computer can quickly run through the algorithm for each possibility, check it, and print out only those squares which meet Charlie's requirements. Such a program is quite straightforward to write: with many algorithms of this type, it is quicker to write a program and run it than to do things by hand.

This program is written in Python. Python is a reasonably straightforward language which is not too difficult for a beginner to learn, but is also quite sophisticated and powerful. It is free, and runs on PCs (under MS-DOS, Windows, OS/2), Macs, and many brands of UNIX. You can download it or read much more about it.

So, let's look at how this program works. You will probably find it helpful to have Charlie's article in front of you. As in Charlie's article, the numbers in the square are represented by the letters a to p. I will present the program as it is, and try and explain each line as we go along, and then at the end I will try and summarise what's going on.

def up_to(last): return range(last+1)

def has_no_repeats(square):
  for i in up_to(15):
    for j in up_to(i-1):
      if square[i]==square[j]: return 0
  return 1

These two sections are setting up functions which we will use later in the program. I will explain them later.

a,b,c,d = 10,2,20,01

This just tells the program the values we want across the top of the square.

for m in up_to(b+c):
  p = b+c-m

This relates to step (ii) in Charlie's article. We know that we can choose the values for m and p, as long as they add up to (b+c). So m could be anything from 0 to (b+c), and then p will be (b+c)-m. The first of these two lines starts a "loop". This is a bit of code that will be repeated over and over again, starting with m being 0, then with m being 1, and so on up to m being b+c. The second line just calculates the value of p. We show which code is within the loop by indenting it.

for g in up_to(a+p):
    j = a+p-g

This deals with step (iii), and starts another loop, for different values of g. This loop is inside the first loop: we call them nested loops. So the computer will work first with m=0,g=0, then m=0,g=1, then m=0,g=2 and so on up to m=0,g=(a+p). then when it gets to the end of the g loop, it will go back to the start of the m loop and work with m=1,g=0, then m=1,g=1, then m=1,g=2, etc. This sounds like an awful lot of work, but that's what a computer is for!

for f in up_to(d+m):
      k = d+m-f

Here we start a third nested loop, corresponding to step (iv) in Charlie's article.

n = g+k-b
      if n< 0: continue

In step (v), no choice is involved, as three of the numbers in the equation b+n=g+k have already been decided. All we need to do is check that n doesn't turn out negative. If it is negative, we need to give up on the set of values it's currently working with, as they don't work. "continue" tells the computer to give up working with the current value of f, and go back to the start of the f loop and try the next value of f.

o = f+j-c
      if o< 0: continue

If the value of n was okay, we carry on to step (vi), which tells us how to work out o. Again, if o is negative, we give up and go back to the next value of f.

for h in up_to(a+m):
        l = a+m-h
        e = j+k-h
        i = f+g-l
        if e< 0 or i< 0: continue

The last choice, in stage (vii), means we begin a last loop. After that, we can work out the remaining values, and check that they're not negative. If they are, we go back and try the next value of h.
[You might be wondering why we don't need to check whether l is negative. Remember that we knew that h and l added up to (a+m). Since we're only trying values of h up to (a+m), the lowest l can be is 0.]

At this point, we have given instructions to work out the whole square. If the computer has made it past all those "continue"s, then we have a set of values which are all positive. All we need to do now is to check whether any numbers occur twice, and print out the list values if it's okay.

if has_no_repeats([a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p]):
          print a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p

has_no_repeats is the function which was defined at the start of the program, and if its output is 1, the program will follow the print instruction to print out a list of 16 numbers. If its output is 0, it will ignore the print line, and go back to the start of the h loop.

Here is the program all in one go:

def up_to(last): return range(last+1)

def has_no_repeats(square):
  for i in up_to(15):
    for j in up_to(i-1):
      if square[i]==square[j]: return 0
  return 1

a,b,c,d = 10,2,20,01
for m in up_to(b+c):
  p = b+c-m
  for g in up_to(a+p):
    j = a+p-g
    for f in up_to(d+m):
      k = d+m-f
      n = g+k-b
      if n< 0: continue
      o = f+j-c
      if o< 0: continue
      for h in up_to(a+m):
        l = a+m-h
        e = j+k-h
        i = f+g-l
        if e< 0 or i< 0: continue
        if has_no_repeats([a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p]):
          print a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p

So, the computer sets gaily off, taking its first value for m, its first value for g, its first value for f and its first value for h. As it goes along, it works out the values for all the other letters as well. If they're all positive, it succeeds in getting to the last two lines, and if they're all different, it'll print them out as well. It will then go back and try the next value of h. Once it has tried all the values of h, it will go back to the start of the f loop, and try the next value of f (and all the values of h that could work with that value of f). It will still be on its first values for m and g!

You can perhaps visualise this a bit better if you picture a tree which forks four times along each branch:

The picture shows just a bit of the tree, and there will probably be more branches at each fork, but hopefully you can get the idea:

The computer explores all the possibilities, starting with taking the topmost branch at each fork. Each time it reaches the end of a branch, or a problem, it goes back to the last fork, and takes the next branch. Eventually it will get down to the twig in the bottom right, having printed out all the branches which actually worked!

I did promise to explain the first six lines at some point:

def up_to(last): return range(last+1)

def has_no_repeats(square):
  for i in up_to(15):
    for j in up_to(i-1):
      if square[i]==square[j]: return 0
  return 1

The first function is a bit of a fiddle. Python has a built-in function called "range" for making lists of numbers; for instance, range(3) is the list [0,1,2]. But it's a bit unhelpful that range(3) doesn't include 3, so our function "up_to" is the same, but goes one higher than "range". So up_to(3) is [0,1,2,3].

The second function checks for repeated numbers, using two nested loops. "square[i]" means the i ^th number in the square. (We are using i and j here because that's traditional in this sort of situation. They are not the same i and j as we used earlier in the program, and we could just as easily have used x and y, or p and q.)
If you think through the loops, first the computer checks whether the 0^th number matches any of those before it (there aren't any, but never mind), then whether the 1^st matches the 0^th, then whether the 2^nd matches the 0^th or 1^st, then whether the 3^rd matches the 0^th, 1^st or 2^nd, and so on. If at any point there is a match, the function's output is 0. If we get all the way through with no matches, the output is 1.

You can download the program (you'll need Python first) and play with it. It would be interesting to know which dates have no magic square to suit Charlie. His article gives one restriction: are there others? If you are feeling adventurous, you could try extending the program to run through a series of dates: do use Ask NRICH if you need help with this.

This is not the only way to write a program to follow this algorithm. A future article will show how a recursive function can be used instead of all the loops.