rec.puzzles Archive

Are there large squares, containing only consecutive integers, all of whose rows, columns and diagonals have the same sum?

arithmetic/magic.squares Added material

Problem

Are there large squares, containing only consecutive integers, all of whose
rows, columns and diagonals have the same sum?  How about cubes?
Show solution
These are called magic squares.  A magic square of order n (integers
from 1 to n*n) has only one possible sum: (n*n+1)*n/2.

Odd and even order squares must be constructed by different approaches.
For odd orders, the most common algorithm is a recursive scheme
devised by de la Loubere about 300 years ago.  For even orders, one
procedure is the Devedec algorithm, which treats even orders not
divisible by 4 slightly differently from those which are divisible by
4 (doubly even).

For squares with odd-length sides, the following algorithm builds a magic
square:

Put 1 in the middle box in the upper row. From then on, if it's
possible to put the next number one box diagonally up and to the right
(wrapping around if the edge of the grid is reached), do so, otherwise,
put it directly below the last one.

               17 24  1  8 15
               23  5  7 14 16 
                4  6 13 20 22
               10 12 19 21  3 
               11 18 25  2  9 

...or even

         47 58 69 80  1 12 23 34 45
         57 68 79  9 11 22 33 44 46
         67 78  8 10 21 32 43 54 56
         77  7 18 20 31 42 53 55 66
          6 17 19 30 41 52 63 65 76
         16 27 29 40 51 62 64 75  5
         26 28 39 50 61 72 74  4 15
         36 38 49 60 71 73  3 14 25
         37 48 59 70 81  2 13 24 35

See archive entry knight.tour for magic squares that are knight's tours.

To get a 4x4 square, write the numbers in order across each row, filling
the square...

1  2  3  4
5  6  7  8
9  10 11 12
13 14 15 16

then use the following pattern as a mask:

.  X  X  .
X  .  .  X
X  .  .  X
.  X  X  .

Everywhere there is an X, complement the number (subtract it from
n*n+1).  For the 4x4 you get:

1  15 14 4
12 6  7  9
8  10 11 5
13 3  2  16

For n even (n>4):

Make an initial magic square by writing an n/2 magic square four times
(the same one each time).  Now, although this square adds up right we
have the numbers 1 to n*n/4 written four times each.  To fix this,
simply add to it n*n/4 times one of the following magic squares:

if n/2 is odd (example: n/2=7),

3 3 3 0 0 0 0 2 2 2 2 2 1 1   (there are int(n/4) 3s, int(n/4-1) 1s on each
3 3 3 0 0 0 0 2 2 2 2 2 1 1    row)
3 3 3 0 0 0 0 2 2 2 2 2 1 1
0 3 3 3 0 0 0 2 2 2 2 2 1 1   (this is row int(n/4)+1.  It starts with just
3 3 3 0 0 0 0 2 2 2 2 2 1 1    the one 0)
3 3 3 0 0 0 0 2 2 2 2 2 1 1
3 3 3 0 0 0 0 2 2 2 2 2 1 1
0 0 0 3 3 3 3 1 1 1 1 1 2 2   (the lower half is the same as the upper half
0 0 0 3 3 3 3 1 1 1 1 1 2 2    with 3<->0 and 1<->2 swapped.  This guarantees
0 0 0 3 3 3 3 1 1 1 1 1 2 2    that each number 1-n*n will appear in the
3 0 0 0 3 3 3 1 1 1 1 1 2 2    completed square)
0 0 0 3 3 3 3 1 1 1 1 1 2 2
0 0 0 3 3 3 3 1 1 1 1 1 2 2
0 0 0 3 3 3 3 1 1 1 1 1 2 2

if n/2 is even (example: n/2=4),

0 0 3 3 2 2 1 1  (there are n/4 of each number on each row)
0 0 3 3 2 2 1 1
0 0 3 3 2 2 1 1
0 0 3 3 2 2 1 1
3 3 0 0 1 1 2 2
3 3 0 0 1 1 2 2
3 3 0 0 1 1 2 2
3 3 0 0 1 1 2 2

References:
    "Magic Squares and Cubes"
    W.S. Andrews
    The Open Court Publishing Co.
    Chicago, 1908

    "Mathematical Recreations"
    M. Kraitchik
    Dover
    New York, 1953

Additional material

Here is a formula that defines a magic square of order n, when n is
odd:

M(i,j) = 1 + (((n-1)/2-i+j) mod n) + n*(((n-3)/2+i+j) mod n)

-- Jan Karrman <jan@tdb.uu.se>, Uppsala University

magic square of the sun :

06 32 03 34 35 01
07 11 27 28 08 30
19 14 16 15 23 24
18 20 22 21 17 13
25 29 10 09 26 12
36 05 33 04 02 31

-- Stefano Dupre <saxtv@mbox.vol.it>

> These are called magic squares.  A magic square of order n (integers
> from 1 to n*n) has only one possible sum: (n*n+1)*n/2.
>
> Odd and even order squares must be constructed by different approaches.
> For odd orders, the most common algorithm is a recursive scheme
> devised by de la Loubere about 300 years ago.

The algorithm you described is correct, but it isn't recursive.

-- David Grabiner <grabiner@math.harvard.edu>

There is an enormous quantity of uh, literature?  well, written matter --
about magic squares. These books have good chapters about them:

M. Kraitchik, _Mathematical Recreations_
H. S. M. Coxeter & W. W. R. Ball, _Mathematical Recreations and Essays_

There are several books devoted to magic squares and generalizations:

W. S. Andrews, _Magic Squares and Cubes_
W. H. Benson & O. Jacoby, _New Recreations With Magic Squares_
same authors: _Magic Cubes: New Recreations_

These are all *in print* in well-produced trade paperbacks from Dover Press.

The original poster asked if there were relations between magic squares and
other mathematical structures.  For a positive answer, see my paper,
"Magic Squares and Linear Algebra", in the _American Mathematical Monthly_,
June-July 1991.

-- Chris Henrich <cjhenric@news.delphi.com>

Here are references to several articles that might be of interest to you:

D.B.  Widdis and R.B. Richter, _It's magic! Multiplication theorems for
magic squares_, College Mathematics Journal, v.20, 301-306. (1989)

In this article a "multiplication" for magic squares is given, which is (in TeX
notation): For A and B magic squares of order m and n, respectively,
the composition * is defined so that the ijth entry of A*B is

      (A*B)_{ij} = B_{i\mod n, j\mod m}+n^2(A_{[i/n],[j/n]}-1).

Here [ ] is the greatest integer function.

A mathematician here at UT at San Antonio has shown that with this product the
set of all magic squares becomes a semigroup. He also discusses some special
classes of magic squares (such as pandiagonals) and shows that these classes
form sub-semigroups. Furthermore, he gives sufficient conditions under which a
magic square of order mn can be factored as the product of squares of orders
m  and n. The paper is:

G. P. Wene, _The semigroup of magic squares_. Currently in submission to J. Rec.
Math.

Two other articles, both from J. Rec. Math, that are related to Wene's paper are:

 J. R. Hendricks, _Inlaid odd-order magic squares_, v. 24, 6-11. (1992)

 F.E. Hruska, _Magic squares, matrices, planes and angles_, v. 23, 183-189. (1991).

-- Francis A. (Sandy) Norman <fnorman@lonestar.utsa.edu>  (March 1994)

[ Editor: for an upside-down magic square of two-digit numbers, see
arithmetic/upside.down.magic.square. ]

[ Editor: a magic cube of order 6, in answer to Col. G. L. Sicherman. ]

The following example comes from [1].  Each of the 108 orthogonal rows
sums to 651 as well as the 4 space diagonals.  It consists of all the
integers from 1 to 216.  Other magic cubes of order 6 are also given in
[1] as well as methods of construction.

Reference:

[1] W.S. Andrews, Magic Squares and Cubes, Dover Publications, Inc.
    New York: 1960. Page 197.

Example of a magic cube of order 6:

                  4 139 161  26 174 147
                 85 166 107 188  93  12
                 98 152 138   3 103 157
Level 1         179  17  84 165 184  22
                183  21  13 175  89 170
                102 156 148  94   8 143

                193  58  80 215  39  66
                112  31 134  53 120 201
                125  71  57 192 130  76
Level 2          44 206 111  30  49 211
                 48 210 202  40 116  35
                129  75  67 121 197  62

                 18 153 136 163  23 158
                 99 180   1  82 104 185
                181  19  95 176 171   9
Level 3         100 154 149  14  90 144
                167   5 108 189 172  10
                 86 140 162  27  91 145

                207  72  55  28 212  77
                126  45 190 109 131  50
                 46 208 122  41  36 198
Level 4         127  73  68 203 117  63
                 32 194 135  54  37 199
                113  59  81 216 118  64

                155  20 150  15 169 142
                101 182  96 177  88   7
                  6  87 106 187  92 173
Level 5         141 168 160  25  11 146
                151  16 137  83 105 159
                 97 178   2 164 186  24

                 74 209  69 204  34  61
                128  47 123  42 115 196
                195 114 133  52 119  38
Level 6          60  33  79 214 200  65
                 70 205  56 110 132  78
                124  43 191  29  51 213

-- Stanley Rabinowitz <stan%tools.DEC@DECWRL.ARPA>

[ Editor: Roberto Sierra's posting "magicsquare - NxN Magic Square
Solver" (comp.sources.misc, volume 40, issue 4, October 1993), with his
program copied in full.  The program was posted as a shell archive; the
archive wrapper and the "X" that began each line of the program are
left out. ]

rqureshi@girtab.usc.edu (Rauf) wrote:
>  I was wondering if any of the esteemed reader of this group
>  could answer my question or show me a way to find it.
>
>  How you can write the integers from 1 to 36, in a big square with
>  36 small squares in it, such that the sum of the numbers in diagonal,
>  column or rows is same. with no repitition of numbers.
>
>  I know how to do it for squares of order 5X5, 7X7,... but this 6X6 is
>  a bit tricky for me.

[*Long* reply with C code to follow -- be forwarned!!]

I love puzzles, so this seemed an ideal problem to tackle last
night while eating cold spaghetti with the headphones blaring.
Unfortunately, after much puzzling and a ream of scratch paper
consumed, I couldn't come up with a direct approach to solving
the problem for abitrary NxN cases, nor could I answer the specific
6x6 question posed.

What I *was* able to do fairly rapidly (being a programmer and
not a mathematician by trade) was to hack out an exhaustive search
program in C which finds the solutions to the NxN problem if you throw
enough CPU cycles at the problem.  I'm fairly certain that if you
ran this program overnight, or perhaps for a day or two you'd get
the answers you want.  [Warning: the algorithm is probably of
exponential order, so don't try anything too large early on
unless you have unlimited access to a Cray!].

The code is actually not too bad.  For 3x3 and 4x4 cases, it
generates the solutions quite rapidly (or at least faster than
the *worst* algorithm I could conceive of).  It starts taking
more time for 5x5 and 6x6 and so on.  Still, there's a recursive
pruning algorithm that does a pretty good job in reducing the
number of cases tested.

What would definitely help would be to understand the properties
of these squares in more detail.  I have some hunches about the
basics of these beasts, but not enough to trust to code.

For example, I suspect that all odd-sized squares have fairly
simple solutions in which the 'middle' number in the sequence
must be located at the center of the square.  For 3x3 squares,
the 5 goes in the middle.  For 5x5 square, I'm pretty sure the
13 has to go in the middle -- it makes 'sense' that way.  For
even-sized squares (4x4, 6x6, etc.), there are a lot more
solutions, and they are much weirder, since there is a greater
degree of freedom to where values are placed.

Also, I suspect that all of the solutions for a given square size
have the same sum, but have no proof of this phenomenon.  For example,
the 3x3 solution (and rotations thereof) adds up to 15.  All of
the 4x4 solutions add up to 34.  I'm not sure what the 5x5
solutions should add up to, but I have a hunch the magic number
is 65 (the middle number '13' times the order '5') -- yet to
be verified, so don't quote me.

[ Editor: as the main entry says, the sum is always (n*n+1)*n/2: the n
rows together hold 1+2+...+n*n = n*n*(n*n+1)/2.  For n = 5 that is
65, for n = 6 it is 111. ]

Again, these are only vague hunches... no firm conclusions.

I actually had some code that was already pretty close to what
I needed, so it only took a few hours to whack it into shape.
You should be able to get this to run on any machine with an ANSI-
compliant C compiler.  There's instructions for how to fire the
thing up embedded in the code and a detailed description of the
algorithm.

One nice thing about the code is that you can tell the program
to look for cases only exhibiting a certain sum value, which
will limit the number of cases tested by orders of magnitude.
For example, it takes a couple of minutes on a Sun to come up
with the first 4x4 solution unless you tell it that the 4x4
sum is 34 (all solutions add to 34), in which case you'll
get the same answer back in about a second!!  If you know
what the magic number is for larger systems, 5x5, 6x6 and
7x7 and so on, then that'll help a great deal.  The problem
is, I don't know what the magic numbers are, or even if all
of the solutions share the same sum.

If anyone can come up with a better algorithm, let me know,
I'd be interested in hearing about it.  I know I could do
a better job if I spent more time on it and understood the
properties of these 'magic' squares better.  The nice thing
about having an exhaustive search program lying around is
that you don't have to know too much about the underlying
nature of the problem.  Just leave your system running the
next time you go home for spaghetti.  :-)

-- Roberto Sierra <bert@netcom.com>, Tempered MicroDesigns,
   San Francisco

/*
**  MagicSquare.c   --        Find solutions to the Magic Squares problem.
**                        Roberto Sierra        7/25/93  v1.0
**
**  Description:
**        This program finds all the possible ways that the integers
**        1 through N^2 can be arranged in an NxN square so that the
**        sum of any row, column, or diagonal is identical.  By default,
**        the program prints the first solution it finds.  You can
**        use command line options to have all solutions printed or
**        merely counted.  The program allows the square to be any
**        size from 1x1 (trivial case) to 10x10.  WARNING: the time
**        required by the algorithm blows up sharply as the dimensions
**        increase.
**
**        The code makes no attempt to eliminate symmetrical solutions,
**        so the number of distinct solutions will be much lower than
**        the number of solutions reported by the program.
**
**
**  Usage:
**        MagicSquare [-ac] n [sum]
**
**        n        The dimensions of the square (the length of a side).
**                Must be an integer value in the range 1<=n<=10.
**        sum        Optional sum constraint.  By default, the program
**                will find arrangements matching any sum that appears along
**                all rows, columns and diagonals.  By specifying a specific
**                sum, only those solutions which match that specific sum
**                will be searched for and reported.  This should speed up
**                the program quite significantly.  The sum value, if
**                supplied, must be an integer in the proper range --
**                [FYI -- n*(n+1)/2 <= sum <= n*(2n^2-n+1)/2]
**        -a        Find (and print) ALL solutions.
**        -c        Count all solutions, but do not print them.  Overrides -a.
**
**        The output is sent to stdout.  All error messages are
**        sent to stderr.  If a problem arises, -1 is returned.
**
**
**  Examples:
**
**        MagicSquare 3                ## Show a 3x3 solution
**        3x3 board...
**        sum=15
**          2  7  6
**          9  5  1
**          4  3  8
**
**        MagicSquare -c 3        ## Count all 3x3 solutions
**        3x3 board...
**        ...there are 8 solutions.
**
**        MagicSquare 4 34        ## Show a 4x4 solution whose sum is 34
**        4x4 square...
**        sum = 34
**          1 14  8 11
**         15  4 10  5
**         12  7 13  2
**          6  9  3 16
**
**        MagicSquare -a 4 10        ## Show all 4x4 solutions whose sum is 10
**        4x4 square...
**        ...there are 0 solutions whose sum is 10
**
**
**  Build Instructions:
**        You'll need an ANSI C compiler (or the willingness to edit
**        the program a bit).  If you've got Gnu C, then you can
**        compile and load the program as follows:
**
**                gcc MagicSquare.c -ansi -o MagicSquare
**
**        [If you're using MPW on the Mac, define '-d MPW' on the
**        compile line so that background processing will occur.]
**
**
**  Algorithm:
**        A brute-force approach would be to attempt to find all possible
**        ways that N^2 integers can be arranged within an NxN square,
**        compute the row/column/diagonal sums for each position, and
**        report only those arrangments which have equal sums all around.
**        The problem is that N^2 integers can be arranged in (N^2)!
**        ways, which blows up rapidly as N increases.  For N=2,
**        only 24 positions need be tested.  For N=3, some 362,880
**        positions exist (still relatively reasonable).  However,
**        for N=4, roughly 2x10^13 possible positions exist and the
**        numbers get much worse beyond that.  In addition, the cost of
**        computing row/column/diagonal sums for each position will
**        gobble up lots of CPU time.
**
**        Although the approach I've taken is still fairly simplistic
**        and still blows up quite sharply as the size of the square
**        grows, it's not quite on the order of (N^2)!  [I wouldn't
**        hazard to guess what the actual order is...].  Fortunately,
**        I already had a piece of recursive code which I was able to
**        hack into shape in the space of an evening to do the job, so
**        I'm quite satisfied with the result for that reason alone.
**        [FYI -- The other code computed solutions to the Eight Queens
**        chess problem in an extremely rapid, straightforward manner.]
**
**        The algorithm makes no assumptions about the nature of these
**        'magic' squares, and therefore is able to find ALL possible
**        solutions.  A better algorithm would exploit the specific
**        nature of the magic these squares exhibit.
**
**        My algorithm is as follows:
**        (1)  If the user has supplied a constraining sum, perform
**             step 3, otherwise perform step 2 instead.
**        (2)  Along the 'forward' diagonal, recursively generate
**             all possible ways that N integers can be arranged
**             along the diagonal.  For each position, record the
**             sum of the diagonal and use it as the 'target' for
**             the remainder of the algorithm.  Go on to step 4.
**        (3)  If a constraining sum has been supplied by the user,
**             then along the 'forward' diagonal, recursively generate
**             all possible ways that N integers that add up to the
**             target sum can be arranged.  Go on to step 4.
**        (4)  Now that the 'target' sum value has been computed
**             and a set of diagonal values are in place, pivot on
**             the value in the top-left corner and recursively
**             generate all possible ways that the remaining (N^2)-N
**             integers can be arranged such that the sum of the row
**             matches the target.  Stop the recursion at the first
**             point that we know we're at a dead end (for instance,
**             if the row sum exceeds the target sum or there are
**             no numbers left which are small enough to do the job).
**        (5)  Once the top row has been generated, pivot on the
**             same top-left value and recursively generate a column
**             from the remaining set of integers such that the
**             column sum matches the target.  Once again, stop
**             the recursion at any point where it is clear that
**             the sum won't match the target or the necessary
**             integers are unavailable.
**        (6)  Repeat steps 4 and 5 on successive diagonal pivots
**             until all rows and columns have been generated.
**        (7)  Lastly, when all rows and columns have been generated
**             with matching sums, check that the sum of the reverse
**             diagonal also matches.  If so, print the board layout
**             as a solution and continue in this manner until all
**             possible solutions have been found (-a or -c option).
**
**        Note:
**        The loops are structured to try smaller numbers first
**        and the larger numbers later.  I experimented with flipping
**        the allocation order around and discovered that it typically
**        takes *much* longer to find the first solution if the larger
**        numbers are allocated first.  I guess this is because there
**        is more low-level recursion going on with low numbers last.
**
**        The initial code was somewhat ugly in the way that
**        the sums were computed and the way that unique integers
**        were tracked.  Hopefully I'll be able to spend another
**        evening putting in some better code to make these 'inner
**        loops' significantly faster.
**
**
**  Contact:
**        For queries regarding this program, or to obtain source
**        for the Eight Queens chess program mentioned, contact
**        the author at any of the following addresses:
**
**                Roberto Sierra
**                bert@netcom.com   (preferred address)
**                73557.2101@compuserve.com
**
**                Tempered MicroDesigns
**                P.O. Box 170638
**                San Francisco, CA  94117
**
**
**  Fine Print:
**        This program is in the public domain and can be used for
**        any purpose whatsoever, including commercial application.
**        [I'd like to hear what you do with it, though.]
**        Absolutely no warranty or liability is implied or extended
**        by the author.
**
**
**  Modification History:
**        PRS  7/25/93  v1.0 -- Original ANSI C version.  Inefficient.
**
*/


#include <stdio.h>                                /* Need standard I/O functions */
#include <stdlib.h>                                /* Need exit() routine interface */
#include <string.h>                                /* Need strcmp() interface */
#ifdef        MPW                                        /* Macintosh MPW ONLY */
#include <CursorCtl.h>                                /* Need cursor control interfaces */
#endif

#define MAXSIZE     10                                /* Maximum size of square */
#define MAXVALUE    (MAXSIZE * MAXSIZE)         /* Maximum integer value */

                                                /* ONE of these must be #define'd */
#define        ORIG_SCAN                                /* Select original 'uniqueness' method */
#undef        LUT_SCAN                                /* Faster scanner -- limited to 5x5 max */




/******************************/
/****        GLOBAL VARIABLES   ****/
/******************************/

int                n = 0;                                /* Size of square */
int                nsquare = 0;                        /* Maximum integer value (n*n) */
int                printing = 1;                        /* TRUE if printing positions */
int                findall = 0;                        /* TRUE if finding all solutions */

unsigned long        solutions = 0;                        /* Number of solutions found */
int                board[MAXSIZE][MAXSIZE];        /* Shows current contents of board */
int                targ = 0;                        /* Target col/row/diagonal sum value */
char                *progname = 0;                        /* The name of this program */
char                usage[] = {"usage: %s [-ac] n [sum]\n"};
                                                /* Usage string */


/*******************/
/****        MACROS        ****/
/*******************/

/*
**  ORIG_SCAN scanning method
**  This is a series of macros which scan through integers
**  such that no integer is ever used more than once.  Since
**  this stuff sits at the 'bottom' of the recursion, the loops
**  tend to slow everything down.  [This was the original scanning
**  method, and was the easiest to implement.]
*/

#ifdef        ORIG_SCAN                                /* Start of ORIG_SCAN implementation */

int        inuse[MAXVALUE+1];                        /* Shows which numbers are in use */

#define RESERVE(i)            (inuse[i] = i)        /* Mark value as in use */
#define RELEASE(i)            (inuse[i] = 0)        /* Mark value as not in use */
/* Looping primitives */
                                                /* Loop through all values */
#define FOR_ALL(i)            for(i=1; i<=nsquare; ++i)
                                                /* Find all available values */
#define FOR_AVAILABLE(i)    FOR_ALL(i) if (!inuse[i])
                                                /* Find all available thru maximum */
#define FOR_AVAILABLE_RESTRICTED(i,high)    \
                            for(i=1; i<=(high); ++i) if (!inuse[i])

void init_scan(void)                                /* Initializes scanning variables */
{
    register int i;                                /* Loop variable */
    FOR_ALL(i) RELEASE(i);                        /* Release all values */
}                                                /* End of init_scan() */

#endif        /*ORIG_SCAN*/                                /* End of ORIG_SCAN implementation */



/*
**  LUT_SCAN scanning method
**  [not yet implemented]
*/

#ifdef        LUT_SCAN                                /* Start of LUT_SCAN implementation */
#error        LUT_SCAN not yet implemented.  Use ORIG_SCAN.
#endif        /*LUT_SCAN*/                                /* End of LUT_SCAN implementation */





/***********************/
/****        ROUTINES    ****/
/***********************/

/* Internal prototypes */
void        main(int argc,char **argv);                    /*        Main program */
void        free_diagonal(int i,int sum);                    /*        Diagonal recursion */
void        constrained_diagonal(int i,int sum);            /*        Diagonal to match constraint
*/
void        row(int i);                                    /*        Main row routine */
void        column(int i);                                    /*        Main column routine */
void        row_recurse(int i,int j,int sum,int *p);    /*        Row recursion routine */
void        col_recurse(int i,int j,int sum,int *p);    /*        Column recursion routine
*/
void        final_check(void);                            /*        Bottom level solution check */
void        pboard(void);                                    /*        Print a solution */





/*---------------------------- main() -------------------------
**
**  The job of the main program is mostly to decode the command
**  line arguments and initialize all the arrays and variables
**  needed.  Then, the appropriate free or constrained diagonal
**  recursion routine is called.
*/

void main(int argc,char **argv)
{
    register int    i;                                /* Loop variables */
    register char   *p;                         /* Pointer to argument */

#ifdef        MPW                                        /* Macintosh MPW ONLY */
    InitCursorCtl(0);                                /* Enable cursor control */
#endif

    progname = argv[0];                         /* The name of the program */

    /****   DECODE COMMAND LINE ARGUMENTS   ****/

    for (i=1; i<argc; ++i) {                        /* Scan through arguments */
        p = argv[i];                                /* Pointer to base of argument */
        if (*p == '-') {                        /* Command line option? */
            while (*++p) {                        /* Loop through characters */
                switch (*p) {                        /* What is the character */
                case 'a':                        /* '-a' option */
                    findall = 1;                /* Set flag to find all solutions */
                    break;
                case 'c':                        /* '-c' option */
                    printing = 0;                /* Counting, not printing */
                    findall = 1;                /* Also forces findall option */
                    break;
                default:                        /* Illegal option */
                    fprintf(stderr,"%s: Illegal option '%s'\n",progname,argv[i]);
                    fprintf(stderr,usage,progname);
                    exit(-1);
                }                                /* End of switch */
            }                                        /* End of loop */
        } else {                                /* End of option test */
            if (n == 0) {                        /* Is N unitialized? */
                if (sscanf(p,"%d",&n) != 1) {        /* Read N value */
                    fprintf(stderr,"%s: '%s' is not an integer\n", progname, p);
                    exit(-1);
                }
                if (n <= 0 || n > MAXSIZE) {        /* Range check */
                    fprintf(stderr,"%s: 1st number, %d, not in proper range\n",
                        progname, n);
                    fprintf(stderr,"%s: should be in range 1..%d\n",progname,MAXSIZE);
                }
            } else if (targ == 0) {                /* Is target sum uninitialized? */
                    int lowsum = n*(n+1)/2;                /* Sum of integers 1..N */
                int highsum = n*(2*n*n-n+1)/2;        /* Sum of integers N^2-N+1..N^2 */

                if (sscanf(p,"%d",&targ) != 1) {
                    fprintf(stderr,"%s: '%s' is not an integer\n", progname, p);
                    exit(-1);
                }
                if (targ < lowsum || targ > highsum) {
                    fprintf(stderr,
                        "%s: 2nd number, %d, not in proper range\n",
                        progname, targ);
                    fprintf(stderr,
                        "%s: should be in the range %d..%d for a %dx%d square\n",
                        progname, lowsum, highsum, n, n);
                    exit(-1);
                }
            } else {                                /* Too many arguments */
                fprintf(stderr,"%s: unexpected argument '%s'\n", progname, p);
                fprintf(stderr,usage,progname);
                exit(-1);
            }                                        /* End of value if-then-else */
        }                                        /* End of argument test */
    }                                                /* End of argument scan loop */
    if (n == 0) {
        fprintf(stderr,"%s: parameter N is missing\n",progname);
        fprintf(stderr,usage,progname);
        exit(-1);
    }


    /*        INITIALIZATION        */
    nsquare = n*n;                                /* This is the highest integer */
    init_scan();                                /* Initialize scanning variables */

    printf("%dx%d square...\n", n, n);
    fflush(stdout);                                /* Force prompt to output */
    if (targ == 0) free_diagonal(0,0);                /* Begin unconstrained recursion */
    else constrained_diagonal(0,0);                /* Otherwise, use constrained sum */


    /*        REPORT NUMBER OF SOLUTIONS FOUND  */
    if (printing && solutions) putchar('\n');
    if (solutions == 1) {
        printf("...there is %ld solution", solutions);
    } else {
        printf("...there are %ld solutions", solutions);
    }
    if (targ) printf(" whose sum is %d", targ);
    printf("\n");

    exit(0);                                        /* No errors */
}                                                /* End of main() */




/*---------------------- free_diagonal() -----------------------
**
**  This fills in the forward diagonal values and computes
**  the sum of the diagonal as it proceeds.  It is recursive
**  so that all possible combinations are attempted.  When
**  the diagonal is full, then row/column filling is performed.
*/

void free_diagonal(register int pivot, int sum)
{
    register int    x;                                /* Value under test */
    register int    *p = &board[pivot][pivot];        /* Address of pivot point */

    if (pivot++ == n) {                                /* Stop recursion */
#ifdef        MPW                                        /* Macintosh MPW ONLY */
        SpinCursor(1);                                /* Enable background processing */
#endif
        targ = sum;                                /* Set target sum for rows/columns */
        row(0);                                        /* Begin row/column fill */
        targ = 0;                                /* Reset to 'floating' sum */
    } else {                                        /* Diagonal not full yet */
        FOR_AVAILABLE(x) {                        /* For all available values... */
            *p = RESERVE(x);                        /* Reserve value for pivot */
            free_diagonal(pivot,sum+x);                /* Recurse to next pivot */
            *p = RELEASE(x);                        /* Release pivot value */
        }
    }
}                                                /* End of free_diagonal() */




/*-------------------- constrained_diagonal() -------------------
**
**  This fills in the forward diagonal values in a constrained
**  fashion so that their sum is equal to that of the target
**  value supplied by the user.  It is recursive in nature so
**  that all possible combinations are attempted.  When the
**  diagonal is full, then row/column filling is started.
*/

void constrained_diagonal(register int pivot, int sum)
{
    register int    x;                                /* Value under test */
    register int    *p = &board[pivot][pivot];        /* Address of pivot point */
    register int    range = targ - sum;                /* Range of possible values */

    if (pivot++ == n) {                                /* Stop recursion */
#ifdef        MPW                                        /* Macintosh MPW ONLY */
        SpinCursor(1);                                /* Enable background processing */
#endif
        if (sum == targ) row(0);                /* Begin row/column fill */
    } else {                                        /* Diagonal not full yet */
        if (range>nsquare) range = nsquare;        /* Clip range to m        FOR_AVAILABLE_RESTRICTED(x,range) {        /* For all available <= range */
            *p = RESERVE(x);                        /* Reserve value for pivot */
            constrained_diagonal(pivot,sum+x);        /* Recurse to next pivot */
            *p = RELEASE(x);                        /* Release pivot value */
        }
    }
}                                                /* End of constrained_diagonal() */




/*------------------------- row() ------------------------
**
**  row() is the top of the row recursion process.
**  The current sum of the specified row is computed.
**  Then the row recursion routine is called to fill
**  in the remaining slots with all possible values.
*/

void row(register int i)
{
    register int    j,sum=0;                        /* Column index, row sum */
    register int    *p = &board[i][0];                /* Pointer to first item */

    for (j=0; j<=i; ++j) {                        /* For all columns thru pivot */
        sum += *p++;                                /* Add values in row */
    }
    row_recurse(i,i+1,sum,p);                        /* Begin row recursion */
}                                                /* End of row() */



/*----------------------- column() -----------------------
**
**  COLUMN is the top of the column recursion process.
**  The current sum of the specified column is computed.
**  Then the column recursion routine is called to fill
**  in the remaining slots with all possible values.
*/

void column(register int j)
{
    register int    i,sum=0;                        /* Row index, column sum */
    register int    *p = &board[0][j];                /* Pointer to first item */

    for (i=0; i<=j; ++i) {                        /* For all rows thru pivot */
        sum += *p;                                /* Add values in column */
        p += MAXSIZE;                                /* Advance to next row */
    }
    col_recurse(j+1,j,sum,p);                        /* Begin column recursion */
}                                                /* End of column() */





/*------------------- row_recurse() ----------------------
**
**  ROW_RECURSE is the row recursion routine.
**  This routine fills in the remaining items in a row
**  with all possible values so that the row sum does
**  not exceed the target sum.        If the target sum is
**  acheived exactly, then the row is accepted and column
**  recursion is initiated.
*/

void row_recurse(int i, int j, int sum, register int *p)
{
    register int    x;                                /* Value under test */
    register int    range = targ - sum;                /* Range of possible values */

    if (j++ == n) {                                 /* Stop recursion */
#ifdef        MPW                                        /* Macintosh MPW ONLY */
        SpinCursor(1);                                /* Enable background processing */
#endif

        if (range == 0) {                        /* Accept row if sum is correct */
            if (i+1 == n) final_check();        /* Last row done, final checks */
            else column(i);                        /* Otherwise, do column fill */
        }
    } else {                                        /* Haven't reached bottom yet */
        if (range>nsquare) range = nsquare;        /* Clip range to maximum */
        FOR_AVAILABLE_RESTRICTED(x,range) {        /* For all available <= range */
            *p = RESERVE(x);                        /* Reserve value at position */
            row_recurse(i,j,sum+x,p+1);                /* Recurse to next column */
            *p = RELEASE(x);                        /* Release value from position */
        }
    }
}                                                /* End of row_recurse() */




/*--------------------- col_recurse() -----------------------
**
**  COL_RECURSE is the column recursion routine.
**  This routine fills in the remaining items in a column
**  with all possible values so that the column sum does
**  not exceed the target sum.        If the target sum is
**  acheived exactly, then the recursion ceases and a
**  final check of the reverse diagonal is performed before
**  the solution can be accepted.
*/

void col_recurse(int i,int j, int sum, register int *p)
{
    register int    x;                                /* Value under test */
    register int    range = targ - sum;                /* Range of possible values */

    if (i++ == n) {                                 /* Stop recursion */
        if (range == 0) {                        /* Accept column if sum is correct */
            if (j+1 == n) final_check();        /* Last column done, final checks */
            else row(j+1);                        /* Otherwise do next row fill */
        }
    } else {                                        /* Haven't reached bottom yet */
        if (range>nsquare) range = nsquare;        /* Clip range to maximum */
        FOR_AVAILABLE_RESTRICTED(x,range) {        /* For all available <= range */
            *p = RESERVE(x);                        /* Reserve value at position */
            col_recurse(i,j,sum+x,p+MAXSIZE);        /* Recurse to next row */
            *p = RELEASE(x);                        /* Release value from position */
        }
    }
}                                                /* End of col_recurse() */




/*--------------------- final_check() ---------------------
**
**  FINAL_CHECK tests the reverse diagonal sum.
**  This occurs at the bottom of all pivot, row and column recursion,
**  so if the sum is OK, then a new solution has been found.
*/

void final_check(void)
{
    register int    i, sum=0;                        /* Loop variable, diagonal sum */

    for(i=0; i<n; ++i) {                        /* Loop through all rows (cols) */
        sum += board[i][n-i-1];                        /* Add values along diagonal */
    }
    if (sum == targ) {                                /* Diagonal sum matches desired? */
        ++solutions;                                /* Congrats.  We have a solution */
        if (printing) pboard();                 /* Print the solution if desired */
        if (!findall) exit(0);                        /* Stop if only looking for first */
    }
}                                                /* End of final_check() */




/*----------------------- pboard() ----------------------
**
**  This routines prints the board for a particular solution.
**
*/

void pboard(void)
{
    register int    i,j;                        /* Rank/File indices */

    if (findall) {                                /* Print solution number */
        printf("\nSolution #%lu -- ",solutions);
    }
    printf("sum = %d\n",targ);                        /* Print characteristic sum */
    for(i=0; i<n; ++i) {                        /* Loop through all rows */
        for(j=0; j<n; ++j) {                        /* Loop through all columns */
            printf(" %2d",board[i][j]);         /* Print values at each position */
        }                                        /* End of column loop */
        putchar('\n');                                /* Break line */
    }                                                /* End of row loop */
    fflush(stdout);                                /* Force solution to output */
}                                                /* End of pboard() */





/****        End of MagicSquare.c        ****/

[ Editor: in constrained_diagonal() one line is garbled in the posting:
the comment "/* Clip range to m..." runs into the next line,
"FOR_AVAILABLE_RESTRICTED(x,range) {", apparently two lines run
together in transmission (compare the same two lines in row_recurse()).
It is kept as posted. ]