The Number of “Magic” Squares, Cubes, and Hypercubes 1 Matthias Beck, Moshe Cohen, Jessica Cuomo, and Paul Gribelyuk 1 INTRODUCTION. 5 0 0 The peculiar interest of magic squares and all lusus numerorum in general lies in the fact that 2 they possess the charm of mystery. They appear to betray some hidden intelligence which by a n preconceived plan produces the impression of intentional design, a phenomenon which finds its a close analogue in nature. J Paul Carus [3, p. vii] 2 Magic squares have turned up throughout history, some in a mathematical context, others in ] O philosophical or religious contexts. According to legend, the first magic square was discovered in C Chinaby an unknownmathematician sometime before the firstcentury A.D. It was a magic square . of order three thought to have appeared on the back of a turtle emerging from a river. Other magic h t squares surfaced at various places aroundthe world in thecenturies following their discovery. Some a m of the more interesting examples were recorded in Europe during the 1500s. Cornelius Agrippa wrote De Occulta Philosophia in 1510. In it he describes the spiritual powers of magic squares [ and produces some squares of orders from three up to nine. His work, although influential in the 3 mathematical community, enjoyed only brief success, for the counter-reformation and the witch v 3 hunts of the Inquisition began soon thereafter: Agrippa himself was accused of being allied with 1 the devil. Although this story seems outlandish now, we cannot ignore the strange mystical ties 0 magic squares seem to have with the world and nature surrounding us, above and beyond their 1 0 mathematical significance. 2 Despite the fact that magic squares have been studied for a long time, they are still the subject 0 / of research projects. These include both mathematical-historical research, such as the discovery of h unpublished magic squares of Benjamin Franklin [12], and pure mathematical research, much of t a which is connected with the algebraic and combinatorial geometry of polyhedra (see, for example, m [1], [4], and [11]). Aside from mathematical research, magic squares naturally continue to be an : v excellent source of topics for “popular” mathematics books (see, for example, [3] or [13]). In this i X paper we explore counting functions that are associated with magic squares. r Wedefineasemi-magic square tobeasquarematrixwhoseentries arenonnegative integers and a whose rows and columns (called lines in this setting) sum to the same number. A magic square is a semi-magic square whose main diagonals also add up to the line sum. A symmetric magic square is a magic square that is a symmetric matrix. A pandiagonal magic square is a semi-magic square whose diagonals parallel to the main diagonal from the upper left to the lower right, wrapped around (i.e., continued to a duplicate square placed to the left or right of the given one), add up to the line sum. Figure 1 illustrates our various definitions. We caution the reader about clashing definitions in the literature. For example, some people would reserve the term “magic square” for whatwewillcallatraditional magic square, amagicsquareof ordernwhoseentries aretheintegers 1,2,...,n2. Our goal is to count these various types of squares. In the traditional case, this is in some sense not very interesting: for each order there is a fixed number of traditional magic squares. For 1Appeared in American Mathematical Monthly 110, no. 8 (2003), 707–717. The number of “magic” squares, cubes and hypercubes 2 1 0 1 2 3 0 0 1 2 0 0 1 2 1 0 1 2 0 1 2 1 2 1 0 0 2 1 2 0 1 2 1 0 1 Figure 1: A semi-magic, a magic, and a symmetric pandiagonal magic square. example, there are 880 traditional 4×4 magic squares. The situation becomes more interesting if we drop the condition of traditionality and study the number of magic squares as a function of the line sum. We denote the total number of semi-magic, magic, symmetric magic, and pandiagonal magic squares of order n and line sum t by H (t),M (t),S (t), and P (t), respectively. n n n n We illustrate these notions for the case n = 2, which is not very complicated: here a semi-magic square is determined once we know one entry (denoted by ∗ in Figure 2); a magic square has to have identical entries in each coefficient. Hence ∗ t−∗ t/2 t/2 t−∗ ∗ t/2 t/2 Figure 2: Semi-magic and magic squares for n = 2. 1 if t is even, H (t)= t+1 , M (t) = S (t) = P (t) = 2 2 2 2 0 if t is odd. (cid:26) These easy results already hint at something: the counting function H is of a different character n than the functions M , S , and P . n n n The oldest nontrivial results on this subject were first published in 1915: Macmahon [9] proved that t+3 t+2 2t2+ 2t+1 if 3|t, H (t) = 3 + , M (t) = 9 3 3 4 2 3 0 otherwise. (cid:18) (cid:19) (cid:18) (cid:19) (cid:26) The first structural result for general n was proved in 1973 by Ehrhart [7] and Stanley [15]. Theorem 1 (Ehrhart, Stanley) . The function H (t) is a polynomial in t of degree (n − 1)2 n that satisfies the identities H (−n−t)= (−1)n−1H (t) , H (−1) = H (−2) = ··· = H (−n+1) = 0 . n n n n n The fact that H is a polynomial of degree (n−1)2 was conjectured earlier by Anand, Dumir, and n Gupta [2]. An elementary proof of Theorem 1 can be found in [14]. Stanley also proved analogous results for the counting functions for symmetric semi-magic squares. In this paper, we establish analogues of these theorems for the other counting functions mentioned earlier. A quasi-polynomial Q of degree d is an expression of the form Q(t) =c (t) td+···+c (t) t+c (t) , d 1 0 The number of “magic” squares, cubes and hypercubes 3 where c ,c ,...,c are periodic functions of t and c 6≡ 0. The least common multiple of the 0 1 d d periods of the c is called the period of Q. With this definition we can state the first of our two j main theorems. Theorem 2 ThefunctionsM (t), S (t), andP (t)are quasi-polynomials intofdegrees n2−2n−1, n n n n2/2−n/2−2, and n2−3n+2, respectively, that satisfy the identities M (−n−t) =(−1)n−1 M (t) , M (−1) = M (−2) = ··· = M (−n+1) = 0 , n n n n n S (−n−t)= (−1)n(n−1)/2 S (t) , S (−1) = S (−2) = ··· = S (−n+1) = 0 , (1) n n n n n P (−n−t) = P (t) , P (−1) = P (−2) = ··· = P (−n+1) = 0 . n n n n n We can extend the counting function H for semi-magic squares in the following way. Define n a semi-magic hypercube (also called a quasi-magic hypercube) to be a d-dimensional n × ··· ×n array of nd nonnegative integers that sum to the same number t parallel to any axis; that is, if we denote the entries of the array by m (1 ≤ j ≤ n), then we require n m = t for all j1...jd k jk=1 j1...jd k = 1,...,d. Again we count all such cubes in terms of d, n, and t; we denote the corresponding P enumerating function by Hd(t). So H2 = H , whose properties are stated in Theorem 1. Except n n n for the case d = n = 3, which seems to have appeared first in [5], we could not find any other references. We prove the following result for general d. Theorem 3 The function Hd(t) is a quasi-polynomial in t of degree (n − 1)d that satisfies the n identities Hd(−n−t) = (−1)n−1 Hd(t) , Hd(−1) = Hd(−2) = ··· = Hd(−n+1) = 0 . n n n n n We prove Theorem 2 in Sections 2 and 3. To be able to compute the counting functions M , n S , and P for specific n, we need the periods of these quasi-polynomials. We describe in Section 4 n n methodsforfindingtheseperiodsandhencefortheactual computation ofM ,S , andP . Finally, n n n we prove Theorem 3 in Section 5. All of our methods are based on the idea that one can interpret thevariouscounting functionsas enumeratinginteger points (“lattice points”)incertain polytopes. 2 SOME GEOMETRIC COMBINATORICS. A convex polytope P in Rd is the convex hull of finitely many points in Rd. Alternatively (and this correspondenceis nontrivial [16]), one can defineP as the boundedintersection of affinehalfspaces. If there is a hyperplane H = x∈ Rd : a·x= b such that P ∩H consists of a single point, then this point is a vertex of P. A polytope is rational if all of its vertices have rational coordinates. (cid:8) (cid:9) SupposethatP isaconvex rationalpolytopeinRd. Forapositiveinteger t,letL (t)denotethe P number of lattice points in the dilated polytope tP = {tx :x∈ P} (see Figure 3). Similarly, define L∗(t) to be the number of lattice points in the (relative) interior of tP. Ehrhart, who initiated the P study of the lattice point count in dilated polytopes, proved the following [6]. Theorem 4 (Ehrhart) . If P is a convex rational polytope, then the functions L (t) and L∗(t) P P are quasi-polynomials in t whose degree is the dimension of P and whose period divides the least common multiple of the denominators of the vertices of P. In particular, if P has integer vertices, then L and L∗ are polynomials. Ehrhart conjectured and P P partially proved thefollowing reciprocity law, which was proved byMacdonald [8](for thecasethat P has integer vertices) and McMullen [10] (for the case that P has arbitrary rational vertices). The number of “magic” squares, cubes and hypercubes 4 Figure 3: Dilation of a polytope. Theorem 5 (Ehrhart-Macdonald reciprocity law) . L (−t)= (−1)dimPL∗(t) P P Theorems 4 and 5 mark the beginning of our journey towards a proof of Theorem 2. The connection to our counting functions M , S , and P is the following: all the various magic square n n n conditions are linear inequalities (the entries are nonnegative) and linear equalities (the entries in each line sum up to t) in the n2 variables forming the entries of the square. In other words, what we are computing is the number of nonnegative integer solutions to the linear system A x= b , where x is a column vector in Rn2 whose components are the entries of our square, A denotes the matrix determining the magic-sum conditions, and b is the column vector each of whose entries is t. We will use the convention that the entries of x are ordered by rows of the original square, and from left to right in each row, that is, the entries x of the square are being arranged as jk x= (x ,x ,...,x ,x ,x ,...,x ,...,x ,x ,...,x ) . 11 12 1n 21 22 2n n1 n2 nn The number of “magic” squares, cubes and hypercubes 5 The matrix A is constructed accordingly. For example, if we study magic 3×3 squares, 1 1 1 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0   0 0 0 0 0 0 1 1 1  1 0 0 1 0 0 1 0 0  A =   .  0 1 0 0 1 0 0 1 0     0 0 1 0 0 1 0 0 1     1 0 0 0 1 0 0 0 1     0 0 1 0 1 0 1 0 0      The first row of Ax = b represents the first row sum in the square: x + x + x = t; the 11 12 13 second and third row (x +x +x = t and x +x +x = t) are the sums of the second and 21 22 23 31 32 33 third row in the square. The fourth row of Ax = b represents the first column sum in the square: x +x +x = t; and the next two rows take care of the second and third column in our square. 11 21 31 Finally the last two rows of Ax= b represent the two diagonal constraints x +x +x = t and 11 22 33 x +x +x =t. 13 22 31 Furthermore, by writing b = t 1, where 1 signifies a column vector each of whose entries is 1, we can see that our counting function enumerates nonnegative integer solutions to the linear system A x= t 1, that is, it counts the lattice points in the t-dilate of the polytope P = x = (x ,...,x ) ∈Rn2 : x ≥ 0 for k = 1,2,...,n2, A x = 1 , 1 n2 k n o provided that we choose the matrix A according to the magic-sum conditions. Note that P is an intersection of half-spaces and hyperplanes, and is therefore convex. No matter which counting function the matrix A corresponds to, the entries of A are 0s and 1s. We obtain the vertices of P by converting some of the inequalities x ≥ 0 to equalities. It is easy to conclude from this k that the vertices of P are rational. Hence by Theorem 4, M , S , and P are quasi-polynomials n n n whose degrees are the dimensions of the corresponding polytopes. Geometrically, the “magic-sum variable” t is the dilation factor of these polytopes. The Ehrhart-Macdonald reciprocity law (Theorem 5) connects the lattice-point count in P to that of the interior of P. In our case, this interior (again, using any matrix A suitable for one of the counting functions) is described by Pint = x= (x ,...,x ) ∈ Rn2 : x > 0 for k = 1,2,...,n2, A x = 1 . 1 n2 k n o The lattice-point count is the same as before, with the difference that we now allow only positive integers as solutions to the linear system. This motivates us to define by M∗(t),S∗(t), and P∗(t) n n n thecounting functions formagic squares,symmetric magic squares, andpandiagonal magic squares as before, but with the restriction that the entries be positive integers. Theorem 5 then asserts, for example, that M (−t) =(−1)deg(Mn)M∗(t) . (2) n n On the other hand, by its very definition M∗(t) = 0 for t = 1,2,...,n − 1. Hence we obtain n M (−t) = 0 for such t. Also, since each row of the matrix A definedby some magic-sum conditions n has exactly n 1s and all other entries 0, it is not hard to conclude that M∗(t) = M (t − n). n n Combining this with the reciprocity law (2), we obtain M (t) = (−1)deg(Mn)M (−t−n) . n n The number of “magic” squares, cubes and hypercubes 6 There are analogous statements for S and P . Once we verify the degree formulas of Theorem 2, n n (1) follows from the observation that (−1)n2−2n−1 = (−1)n−1 , (−1)12n2−21n−2 = (−1)n(n−1)/2 , (−1)n2−3n+2 = (−1)n(n−3) = 1 . 3 PROOF OF THE DEGREE FORMULAS. Let’s start with magic squares: by Theorem 4, the degree of M is n2 (the number of places we n have to fill) minus the number of linearly independent constraints. In other words, we have to find the rank of the n2×(2n+2) matrix 1 | 1 | | 1 ... | ... | ··· | ...   1 | 1 | | 1 − − − − − − − − − − − − −    1 ··· 1 | | |    A = | 1 ··· 1 | |  .  | | ... |     | | | 1 ··· 1    − − − − − − − − − − − − −    1 | 1 | ··· | 1     1 | 1 | ··· | 1      Here we only show the entries that are 1; every other entry is 0. Similarly as in our example, the first n rows of A represent the n column constraints of our square (x +···+x = t, ..., x + 11 n1 1n ···+x = t), the next n rows represent the row contraints of the square, and the last two rows nn take care of the diagonal constraints. The sum of the first n rows of A is the same as the sum of the next n rows, so one of the first 2n rows is redundant and we can eliminate it; we choose the (n+1)th. Furthermore, we can add the difference of the first and (n+2)th row to the (2n+1)th and subtract the nth row from the (2n+2)th. These operations yield the n2×(2n+1) matrix 1 | 1 | | 1 ... | ... | ··· | ...   1 | 1 | | 1 − − − − − − − − − − − − −    | 1 ··· 1 | |   | | ... |  .    | | | 1 ··· 1    − − − − − − − − − − − − −    | 0 2 1 ··· 1 | ∗ | ∗     | 1 −1 | ∗ | ∗      These represent 2n+1 linearly independent restrictions on our magic square. The polytope corre- spondingto M therefore has dimension n2−2n−1: it lies in an affine subspace of that dimension, n The number of “magic” squares, cubes and hypercubes 7 and the point (1/n,...,1/n) is an interior point of the polytope. Hence, by Theorem 4, the degree of M is n2−2n−1. n The case of symmetric magic squares is simpler: here we have n k = n(n+1)/2 places to k=1 fill, and the row-sum condition is equivalent to the column-sum condition. Hence there are n+2 P constraints, whichareeasily identified aslinearly independent. Thedimensionofthecorresponding polytope, and the degree of S , is therefore n2/2−n/2−2. n Finally, we discuss pandiagonal magic squares. Again we have n2 places to fill; this time the constraints are represented by the n2×3n matrix 1 | 1 | | 1  ... | ... | ··· | ...  1 | 1 | | 1  − − − − − − − − − − − − − − − −     1 | 1 | | 1   1 | ... | ··· | 1   ... | 1 | | ...  .    1 | 1 | | 1   − − − − − − − − − − − − − − − −     1 ··· 1 | | |     | 1 ··· 1 | |   ..   | | . |     | | | 1 ··· 1      Herethefirstnrowsrepresentthecolumnconstraintsinoursquare,thenextnrowsthepandiagonal constraints, and the last n rows the row constraints. As in the first case, the sum of the first n rows equals the sum of the next n rows, as well as the sum of the last n rows. We can therefore eliminate the (2n)th and (2n+1)th rows to get the n2×(3n−2) matrix 1 | 1 | | 1  ... | ... | ··· | ...  1 | 1 | | 1  − − − − − − − − − − − − − − − −     1 | 1 | | 1   ... | ... | ··· | 1  .  1 | 1 | | ... 1   − − − − − − − − − − − − − − − −     | 1 ··· 1 | |   ..   | | . |     | | | 1 ··· 1      In this new matrix, we can now replace the (n+1)th row by the difference of the (n+1)th and The number of “magic” squares, cubes and hypercubes 8 first rows, the (n+2)th row by the difference of the (n+2)th and second rows, and so on: 1 | 1 | | 1 ... | ... | ··· | ...   1 | 1 | | 1 − − − − − − − − − − − − −    | −1 1 | |   | ... | ∗ | ∗  .    | −1 1 | |   − − − − − − − − − − − − −     | 1 ··· 1 | |   ..   | | . |     | | | 1 ··· 1      Finally, we can replace the (n + 2)th row by the sum of the (n + 1)th and (n + 2)th rows, the (n+3)th row by the sum of the (n+2)th and (n+3)th rows, etc. Subtracting the (2n)th row from the (n + 1)th row gives a matrix of full rank, that is, rank 3n − 2. The dimension of the corresponding polytope, which is the degree of P , is thus n2−3n+2. n This finishes the proof of Theorem 2. 4 COMPUTATIONS. To interpolate a quasi-polynomial of degree d and period p, we need to compute p(d+1) values. The periods of M , S , and P are not as simple to derive as their degrees. What we can do, n n n however, is to compute for fixed n the vertices of the respective polytopes, whose denominators give the periods of the quasi-polynomials (Theorem 4). This is easy for very small n but gets complicated very quickly. For example, we can practically find by hand that the vertices of the polytope corresponding to S are (2/3,0,1/3,1/3,2/3,0) and (0,2/3,1/3,1/3,0,2/3), whereas the 3 polytope corresponding to S has seventy-four vertices. To make computational matters worse, the 5 least common multiple of the denominators of these vertices is 60; one would have to do a lot of interpolation to obtain S . 5 The reciprocity laws in Theorem 2 essentially halve the number of computations for each quasi- polynomial. Nevertheless, the task of computing our quasi-polynomials becomes impractical with- out further tricks or large computing power. The following table contains data about the vertices of the polytopes corresponding to M , S , and P for n = 3 and 4. It was produced using Maple. n n n Polytope corresponding to Number of vertices l.c.m. of denominators M 4 3 3 S 2 3 3 P 3 1 3 M 20 2 4 S 12 4 4 P 28 2 4 With this information, it is now easy (that is, easy with the aid of a computer) to interpolate The number of “magic” squares, cubes and hypercubes 9 each quasi-polynomial. Here are the results: 2t2+ 2t+1 if 3|t, 9 3 M (t) = 3 ( 0 otherwise; 2t+1 if 3|t, 3 S (t) = 3 ( 0 otherwise; P (t) = 1t2+ 2 +1; 3 2 3 1 t7+ 7 t6+ 89 t5+ 11t4+ 49t3+ 38t2+ 71t+1 if t is even, 480 240 480 16 30 15 30 M (t) = 4 ( 1 t7+ 7 t6+ 89 t5+ 11t4+ 779t3+ 593t2+ 1051t+ 13 if t is odd; 480 240 480 16 480 240 480 16 5 t4+ 5 t3+t2+ 3t+1 if t ≡ 0 mod 4, 128 16 2 S (t) = 5 t4+ 5 t3+t2+ 3t+ 7 if t ≡ 2 mod 4, 4  128 16 2 8   0 if t is odd;   7 t6+ 7 t5+ 23t4+t3+ 341t2+ 31t+1 if t is even, 1440 120 72 180 15 P (t) = 4 ( 0 if t is odd. 5 SEMI-MAGIC HYPERCUBES. We now prove Theorem 3. The fact that Hd is a quasi-polynomial, the reciprocity law for this n function and Hd∗, and the location of special zeros for Hd follow in exactly the same way as the n n respective statements in Theorem 2. As in that case, theremainingtask is to findthedegree ofHd, n that is, the dimension of the corresponding polytope. Again we have only to find the dimension of the affine subspace of Rnd in which this polytope lives. We do so by counting linearly independent constraint equations. To this end, let us write the coordinates of a point in the polytope corresponding to Hd as n c(a ,...,a ) with 1 ≤ a ≤ n. These are real numbers that satisfy the constraints 1 d j n c(a ,...,a ) ≥ 0 , c(a ,...,a ,j,a ,...,a ) = 1 (1 ≤ k ≤ d) . 1 d 1 k−1 k+1 d j=1 X Oncethe(n−1)d coordinatesc(a ,...,a )with1 ≤ a ≤ n−1arechosen,theremainingcoordinates 1 d j areclearlydeterminedbytheforegoingconditions. Therefore,thereareatleastnd−(n−1)d linearly independent constraint equations. On the other hand, every constraint involves at least one of the nd − (n − 1)d coordinates c(a ,...,a ) for which at least one a equals n. Consider, say, the coordinate c(a ,...,a ,n,...,n) 1 d j 1 k The number of “magic” squares, cubes and hypercubes 10 in which 0 ≤ k < d and a ,...,a < n. We compute: 1 k n−1 c(a ,...,a ,n,...,n) = 1− c(a ,...,a ,n,...,n,j ) 1 k 1 k d jXd=1 n−1 n−1 = 1− 1− c(a ,...,a ,n,...,n,j ,j ) 1 k d−1 d   jXd=1 jdX−1=1   = 1−(n−1)+ c(a ,...,a ,n,...,n,j ,j ) 1 k d−1 d 1≤jd−X1,jd≤n−1 = ... = 1−(n−1)+(n−1)2−···+(−1)d−k−1(n−1)d−k−1 +(−1)d−k c(a ,...,a ,j ,...,j ) 1 k k+1 d 1≤jk+1X,...,jd≤n−1 1−(1−n)d−k = +(−1)d−k c(a ,...,a ,j ,...,j ) . 1 k k+1 d n 1≤jk+1X,...,jd≤n−1 We could have changed the order of summation, that is, used the d−k constraints that are in force here in a different order. However, we would always end up with the same expression. This implies that all constraints involving c(a ,...,a ,n,...,n) with 0 ≤ k < d and a ,...,a < n, but not 1 k 1 k involving any coordinate with a = 1 for some j with j ≤ k, are equivalent to the one constraint j 1−(1−n)d−k c(a ,...,a ,n,...,n) = +(−1)d−k c(a ,...,a ,j ,...,j ) . 1 k 1 k k+1 d n 1≤jk+1X,...,jd≤n−1 Therefore there are at most nd−(n−1)d linearly independent constraint equations. Accordingly, the dimension of the polytope and the degree of Hd are (n−1)d. n 6 CLOSING REMARKS. One big open problem is to determine the periods of our counting functions. The evidence gained from our data seems to suggest that the periods increase in some fashion with n. We believe the following is true: Conjecture 1 . The functions M , S , and P are not polynomials when n ≥ 5. n n n Resultsin [1]lendfurthercredencetothevalidity of thisconjecture. Asforsemi-magic hypercubes, B´ona showed in [5]that H3 really is a quasi-polynomial, not a polynomial. We challenge the reader 3 to prove: Conjecture 2 . The function Hd is not a polynomial when n,d≥ 3. n ACKNOWLEDGEMENTS.WearegratefultoJesu´sDeLoera,JosephGallian,RichardStanley, Thomas Zaslavsky, and the referees for helpful comments and corrections on earlier versions of this paper.