Six Little Squares and How Their Numbers Grow

Matthias Beck^{1} Department of Mathematics San Francisco State University 1600 Holloway Avenue San Francisco, CA 94132 USA

                  beck@math.sfsu.edu

Thomas Zaslavsky^{2} Department of Mathematical Sciences Binghamton University Binghamton, NY 13902-6000 USA

              zaslav@math.binghamton.edu

Abstract

We count the 3 × 3 magic, semimagic, and magilatin squares, as functions either of the magic sum or of an upper bound on the entries in the square. Our results on magic and semimagic squares differ from previous ones, in that we require the entries in the square to be distinct from each other and we derive our results not by ad hoc reasoning, but from the general geometric and algebraic method of our paper “An enumerative geometry for magic and magilatin labellings”. Here we illustrate that method with a detailed analysis of 3 × 3 squares.

1Research partially supported by National Science Foundation grant DMS-0810105.

2Research partially supported by National Science Foundation grant DMS-0070729 and by the SGPNR.

1 Introduction

“Today, the study of magic squares is not regarded as a subject of mathematics, but many earlier mathematicians in China and Japan studied it.” These words from Shigeru’s history of old Japanese mathematics [14, p. 435] are no longer completely true. While the construction of magic squares remains for the most part recreational, their counting has become part of the mainstream of enumerative combinatorics, as an example of quasipolynomial counting formulas and as an application of Ehrhart’s theory of lattice points in polytopes. There are several classical [12, 11, 16] and recent [1, 2] mathematical works on counting something like magic squares, but without the requirement that the entries be distinct, and often omitting the diagonals. In previous articles [5, 6] we established the groundwork for an enumerative theory of magic squares with distinct entries. Here we apply those geometrical and algebraic methods to solve the problem of counting three kinds of magical 3 × 3 squares.

Each square x = (x_{jk})_{3×3} has positive integral entries that satisfy certain line-sum equations and distinctness conditions. In a weakly semimagic square every row and column sum is the same (their common value is called the magic sum); in a weakly magic square each of the two diagonals also adds up to the magic sum. Such squares have been studied before (see, e.g., Beck et al. [2] and Stanley [17]); the difference here is that we count strongly magic or semimagic squares, where all entries of the square are distinct. (Since strongly magic squares are closest to what are classically known as “magic squares”—see the introduction to our general magic article [6]—we call strong squares simply “magic” or “semimagic” without qualification.) The third type we count is a magilatin square; this is a weakly semimagic square with the restriction that the entries be distinct within a row or column. The numbers of standard magic squares (with entries 1, 2, …, n^{2}) and latin squares (in which each row or column has entries 1, 2, …, n) are special evaluations of our counting functions.

We count the squares in two ways: by magic sum (an affine count), and by an upper bound on the numbers in the square (a cubic count). Letting N(t) denote the number of squares in terms of a parameter t which is either the magic sum or a strict upper bound on the entries, we know by our previous work [6] that N is a quasipolynomial, that is, there are a positive integer p and polynomials N_{0}, N_{1}, …, N_{p−1} so that

N(t) = N_{t (mod p)}(t) .

The minimal such p is the period of N; the polynomials N_{0}, N_{1}, …, N_{p−1} are the constituents of N, and N_{0} is the principal constituent. Here we find an explicit list of constituents and also the explicit rational generating function N(x) = ∑_{t>0} N(t)x^{t} (from which the quasipolynomial is easily extracted).

Each magic and semimagic square also has an order type, which is the arrangement of the cells in order of increasing value of their entries. The order type is a linear ordering of the cells because all entries are distinct in these squares. There are 9! = 362880 possible linear orderings but only a handful are order types of squares. Our approach finds the actual number of order types for each kind of square; it is the absolute value of the constant term of the principal constituent, that is, |N_{0}(0)|. (See Theorems 3.4 and 3.14 and Examples 3.11, 3.12, and 3.21 in our paper on magic labellings [6].) Obviously, this number will be the same for cubic and affine counts of the same kind of square. (There is also an order type for magilatin squares, which is a linear ordering only within each row and column. As it is not a simple permutation of the cells, we shall not discuss it any further.)

One of our purposes is to illustrate the technique of our general treatment [6]. Another is to provide data for the further study of magic squares and their relatives; to this end we list the exact numbers of each type for small values of the parameter and also the numbers of symmetry types, reduced squares, and reduced symmetry types of each type (and we refer to the On-Line Encyclopedia of Integer Sequences (OEIS) [15] for the first 10,000 values of each counting sequence). A square is reduced by subtracting the smallest entry from all entries; thus, the smallest entry in a reduced square is 0. A square is normalized by being put into a form that is unique in each symmetry class. Clearly, the number of normalized squares, i.e., of equivalence classes under symmetry, is fundamental; and the number of reduced, normalized squares is more fundamental yet.

There are other ways to find exact formulas. Xin [18] tackles 3×3 magic squares, counted by by magic sum, using MacMahon’s partition calculus. He gets a generating function that agrees with ours (thereby confirming both). Stanley’s idea of Möbius inversion over the partition lattice [17, Exercise 4.10] is similar to ours in spirit, but it is less flexible and requires more computation. Beck and van Herick [3] have counted 4 × 4 magic squares using the same basic geometrical setup as ours but with a more direct counting method.

Our paper is organized as follows. Section 2 gives an outline of our theoretical and computational setup, as well as some comments on checks and feasibility. In Section 3 we give a detailed analysis of our computations for counting magic 3 × 3 squares. Sections 4 and 5 contain the setup and the results of similar computations for 3 × 3 semimagic and magilatin squares. We conclude in Section 6 with some questions and conjectures.

We hope that these results, and still more the method, will interest both magic squares enthusiasts and mathematicians.

2 The technique

2.1 General method

The means by which we solve the specific examples of 3×3 magic, semimagic, and magilatin squares is inside-out Ehrhart theory [5]. That means counting the number of 1/t-fractional points in the interior of a convex polytope P that do not lie in any of a certain set ℋ of hyperplanes. The number of such points is a quasipolynomial function E_{P^{∘},ℋ}^{∘} (t), the open Ehrhart quasipolynomial of the open inside-out polytope (P^{∘}, ℋ). The exact polytope and hyperplanes depend on which of the six problems it is, but we can describe the general picture. First, there are the equations of magic; they determine a subspace s of all 3 × 3 real matrices which we like to call the magic subspace—though mostly we work in a smaller overall space ℝ^{d} that results from various reductions. Then there is the polytope P of constraints, which is the intersection with s of either a hypercube [0, 1]^{3^{2}} or a standard simplex {x ∈ ℝ^{3^{2}} : x ≥ 0, ∑ x_{jk} = 1}: the former when we impose an upper bound on the magic square entries and the latter when we predetermine the magic sum. The parameter t is the strict upper bound in the former case (which we call cubical due to the shape of P), the magic sum in the latter (which we call affine as P lies in a proper affine subspace). Finally there are the strong magical exclusions, the hyperplanes that must be avoided in order to ensure the entries are distinct—or in the magilatin examples, as distinct as they ought to be. These all have the form x_{ij} = x_{kl}. The combination of P and the excluded hyperplanes forms the vertices of (P, ℋ), which are all the points of intersection of facets of P and hyperplanes in ℋ that lie in or on the boundary of P. Thus, we count as a vertex every vertex of P itself, each point that is the intersection of some facets and some hyperplanes in ℋ, and any point that is the intersection of some hyperplanes and belongs to P, but not intersection points that are outside P. (Points of each kind do occur in our examples.) The denominator of (P, ℋ) is the least common denominator of all the coordinates of all the vertices of (P, ℋ). The period of E_{P^{∘},ℋ}^{∘} (t) divides the denominator; this gives us a known bound on it.

This geometry might best be explained with an example. Let us consider magic 3 × 3 squares,

[ \begin{matrix}x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{matrix} ] ∈ ℤ_{>0}^{3^{2}}.

The magic subspace is

s = \{[ \begin{matrix}x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{matrix} ] ∈ ℝ^{3^{2}} : \begin{matrix}x_{11} + x_{12} + x_{13} = x_{21} + x_{22} + x_{23} \\ = x_{31} + x_{32} + x_{33} = x_{11} + x_{21} + x_{31} \\ = x_{12} + x_{22} + x_{32} = x_{13} + x_{23} + x_{33} \\ = x_{11} + x_{22} + x_{33} = x_{13} + x_{22} + x_{31} \end{matrix} \} .

The hyperplane arrangement ℋ that captures the distinctness of the entries is

ℋ = {(x_{11} = x_{12}) ∩ s, (x_{11} = x_{21}) ∩ s, …, (x_{32} = x_{33}) ∩ s} .

Finally, there are two polytopes associated to magic 3 × 3 squares, depending on whether we count them by an upper bound on the entries:

P_{c} = s ∩ [0, 1]^{3^{2}},

or by magic sum:

P_{a} = s ∩ \{[ \begin{matrix}x_{11} & x_{12} & x_{13} \\ x_{21} & x_{22} & x_{23} \\ x_{31} & x_{32} & x_{33} \end{matrix} ] ∈ ℝ_{≥0}^{3^{2}} : x_{11} + x_{12} + x_{13} = 1\} .

Our cubical counting function computes the number of magic squares all of whose entries satisfy 0 < x_{ij} < t, in terms of an integral parameter t. These squares are the lattice points in

(P_{c}^{∘} \ ⋃ ℋ) ∩ (\frac{1}{t}ℤ)^{3^{2}}.

Our second, affine, counting function computes the number of magic squares with positive entries and magic sum t. These squares are the lattice points in

(P_{a}^{∘} \ ⋃ ℋ) ∩ (\frac{1}{t}ℤ)^{3^{2}}.

In general, the number of squares we want to count, N(t), is the Ehrhart quasipolynomial E_{P^{∘},ℋ}^{∘} (t) of an open inside-out polytope (P^{∘}, ℋ). We obtain the necessary Ehrhart quasipolynomials by means of the computer program LattE [9]. It computes the closed Ehrhart generating function

E_{P} (x) := 1 + ∑_{t=1}^{∞} E_{P} (t)x^{t}  of the values  E_{P} (t) := #(P ∩ (\frac{1}{t}ℤ)^{d} ).

Counting only interior points gives the open Ehrhart quasipolynomial E_{P^{∘}}(t) and its generating function

E_{P^{∘}}(x) := ∑_{t=1}^{∞} E_{P^{∘}}(t)x^{t}.

Since we want the open inside-out Ehrhart generating function E_{P^{∘},ℋ}^{∘} (x) = ∑_{t=1}^{∞} E_{P^{∘},ℋ}^{∘} (t)x^{t}, we need several transformations. One is Ehrhart reciprocity [10], which is the following identity of rational generating functions:

E_{P^{∘}}(t) = (−1)^{1+dim P} E_{P} (x^{−1}) .  (1)

The inside-out version [5, Equation (4.6)] is

E_{P^{∘},ℋ}^{∘} (x) = (−1)^{1+dim P} E_{P,ℋ} (x^{−1}) .  (2)

We need to express the inside-out generating functions in terms of ordinary Ehrhart generating functions. To do that we take the intersection poset

ℒ(P^{∘}, ℋ) := {P^{∘} ∩ ⋂𝒮 : 𝒮 ⊆ ℋ} \ {∅},

which is ordered by reverse inclusion. Note that ℒ(P^{∘}, ℋ) and ℒ(P, ℋ), defined similarly but with P instead of P^{∘}, are isomorphic posets because ℋ is transverse to P; specifically, ℒ(P, ℋ) = {ū : u ∈ ℒ(P^{∘}, ℋ)}, where ū is the (topological) closure of u. Now we have the Möbius inversion formulas [5, Equations (4.7) and (4.8)]

E_{P^{∘},ℋ}^{∘} (x) = ∑_{u∈ℒ(P^{∘},ℋ )} μ(ˆ0, u) E_{u}(x)  (3)

and

E_{P,ℋ} (x) = ∑_{u∈ℒ(P,ℋ )} |μ(ˆ0, u)| E_{u}(x)  (4)

(since ℋ is transverse to P; see our general paper [6]). Here μ is the Möbius function of ℒ(P^{∘}, ℋ) [17].

Thus we begin by getting all the cross-sectional generating functions E_{u}(x) from LattE. Then we either sum them by (4) and apply inside-out reciprocity (2), or apply ordinary reciprocity (1) first and then sum by (3). (We did whichever of these seemed more convenient.) In the semimagic and magilatin counts we need a third step because the generating functions we computed pertain to a reduced problem; those of the original problem are obtained through multiplication by another generating function.

Once we have the generating function we extract the quasipolynomial, essentially by the binomial series. If an Ehrhart quasipolynomial q of a rational convex polytope has period p and degree d, then its generating function q can be written as a rational function of the form

q(x) := ∑_{t≥1} q(t) x^{t} = \frac{a_{p(d+1)} x^{p(d+1)} + a_{p(d+1)−1} x^{p(d+1)−1} + ⋯ + a_{1} x}{(1 − x^{p})^{d+1}}  (5)

for some nonnegative integers a_{1}, a_{2}, …, a_{p(d+1)}. Grouping the terms in the numerator of (5) according to the residue class of the degree modulo p and expanding the denominator, we get

q(x) = \frac{∑_{r=1}^{p} ∑_{j=0}^{d} a_{pj+r} x^{pj+r}}{(1 − x^{p})^{d+1}} = ∑_{r=1}^{p} ∑_{k≥0} [∑_{j=0}^{d} a_{pj+r}\binom{d + k − j}{d}] x^{pk+r} .

Hence the rth constituent of the quasipolynomial q is

q_{r}(t) = ∑_{j=0}^{d} a_{pj+r}\binom{d + \frac{t−r}{p} − j}{d} for r = 1, …, p.

2.2 How we apply the method

The initial step is always to reduce the size of the problem by applying symmetry. Each problem has a normal form under symmetry, which is a strong square. The number of all magic or semimagic squares is a constant multiple of the number of symmetry types, because every such square has the same symmetry group. For magilatin squares, there are several symmetry types with symmetry groups of different sizes, so each type must be counted separately.

Semimagic and magilatin squares also have an interesting reduced form, in which the values are shifted by a constant so that the smallest cell contains 0; and a reduced normal form; the latter two are not strong but are aids to computation. Reduced squares are counted either by magic sum (the “affine” counting rule) or by the largest cell value (the “cubic” count). All reduced normal semimagic squares correspond to the same number of unreduced squares, while the different symmetry types of magilatin square give reduced normal squares whose corresponding number of unreduced squares depends on the symmetry type.

The total number of squares, N(t), and the number of reduced squares, R(t), are connected by a convolution identity N(t) = ∑_{s} f(t − s)R(s), where f is a periodic constant (by which we mean a quasipolynomial of degree 0; we say constant term for the degree-0 term of a quasipolynomial, even though the “constant term” may vary periodically) or a linear polynomial. Writing for the generating functions N(x) := ∑_{t>0} N(t)x^{t} and similarly R(x), and f(x) := ∑_{t≥0} f(t)x^{t}, we have N(x) = f(x)R(x). It follows from the form of the denominator in Equation (5) that the period of N divides the product of the periods of f and R. The reduced number R(t) is, in the semimagic case, a constant multiple of the number n(t) of reduced, normal semimagic squares; in the magilatin case it is a sum of different multiples, depending on a symmetry group, of the number of reduced, normal magilatin squares of each different type T. Each n(t) is the open Ehrhart quasipolynomial E_{Q^{∘},ℐ}^{∘}(t) of an inside-out polytope (Q, ℐ) which is smaller than the original polytope P.

We compute n(x) := ∑_{t>0} n(t)x^{t} from the Ehrhart generating functions E_{u}(x) of the nonvoid sections u of Q^{∘} by flats of ℐ through the following procedure:

1. We calculate the flats and sections by hand.

2. We feed each u into the computer program LattE [9], which returns the closed generating function Eū(x), whose constant term equals 1 because u is nonvoid and convex.

3. With semimagic squares, by Equations (1)–(4) we have the Möbius-inversion formulas

\begin{aligned}n(x) = E_{Q^{∘},ℐ}^{∘}(x) & = ∑_{u∈ℒ}μ_{ℒ} (ˆ0, u)E_{u}(x) \\ & = (−1)^{1+dim Q} ∑_{u∈ℒ}|μ_{ℒ} (ˆ0, u)| Eū(x^{−1}). \end{aligned}  (6)

where ℒ := ℒ(Q^{∘}, ℐ), the intersection poset of (Q^{∘}, ℐ).

The procedure for magilatin squares is similar but taking account of the several types.

2.3 Checks

We check our results in a variety of ways.

The degree is the dimension of the polytope, or the number of independent variables in the magic-sum equations.

The leading coefficient is the volume of the polytope. (The volume is normalized so that a fundamental domain in the affine space spanned by the polytope has unit volume.) This check is also not difficult. The volume is easy to find by hand in the magic examples. In affine semimagic the polytope is the Birkhoff polytope B_{3}, whose volume is well known (Section 4.2). The cubical semimagic volume is not well known but it was easy to find (Section 4.1.1). The magilatin polytopes are the same as the semimagic ones.

The firmest verification is to compare the results of the generating function approach with those of direct enumeration. If we count the squares individually for t ≤ t_{1} where t_{1} ≥ pd, only the correct quasipolynomial can agree with the counts (given that we know the degree d and period p from the geometry). Though t_{1} = pd is too large to reach in some of the examples, still we gain considerable confidence if even a smaller value of t_{1} yields numbers that agree with those derived from the quasipolynomial or generating function. We performed this check in each case.

2.4 Feasibility

Based on our solutions of the six 3×3 examples we believe our counting method is practical. The calculations are simple and readily verified. Linear algebra tells us the degree, geometry tells us the period; we obtain the generating function using the Ehrhart package LattE [9] and then apply reciprocity (Equation (1) or (2)) and Möbius inversion (Equation (3) or (4)), and extract the constituents, all with Maple. The programming is not too difficult.

In the magic square problems we found the denominator by calculating the vertices of the inside-out polytope. Then we took two different routes. In one we applied LattE and Equation (3). In the other we calculated N(t) for small values of t by generating all magic squares, taking enough values of t that we could fit the quasipolynomial constituents to the data. This was easy to program accurately and quick to compute, and it gave the same answer. The programs can be found at our “Six Little Squares” Web site [7].

In principle the semimagic and magilatin problems can be solved in the same two ways. The geometrical method with Möbius inversion gave complete answers in a few minutes of computer time after a simple hand analysis of the geometry (see Section 4). Direct enumeration on the computer proved unwieldy (at best), especially in the affine case, where the period is largest. A straightforward computer count of semimagic squares by magic sum (performed in Maple—admittedly not the language of choice—on a personal computer) seemed destined to take a million years. Switching to a count of reduced normal squares, the calculation threatened to take only a thousand years. These programs are at our “Six Little Squares” Web site [7], as is a complicated “supernormalized formula” that greatly speeds up affine semimagic counting (see Section 4.2.7).

2.5 Notation

We use a lot of notation. To keep track of it we try to be reasonably systematic.

M, m refer to magic squares (Section 3).

S, s and subscript s refer to semimagic squares (Section 4).

L, l and subscript ml refer to magilatin squares (Section 5).

R, r refer to reduced squares (the minimum entry is 0), while M, S, L, et al. refer to ordinary squares (all positive entries).

c refers to “cubic” counting, by an upper bound on the entries.

a refers to “affine” counting, by a specified magic sum.

X (capital) refers to all squares of that type.

x (minuscule) refers to symmetry types of squares, or equivalently normalized squares.

3 Magic squares of order 3

The standard form of a magic square of order 3 is well known; it is

\begin{matrix}α + γ & −α − β + γ & β + γ \\ −α + β + γ & γ & α − β + γ \\ −β + γ & α + β + γ & −α + γ \end{matrix}  (7)

where the magic sum is s = 3γ. Taking account of the 8-fold symmetry, under which we may assume the largest corner value is α + γ and the next largest is β + γ, and the distinctness of the values, we have α > β > 0 and α ≠ 2β. One must also have γ > α + β to ensure positivity.

In this pair of examples, the dimension of the problem is small enough that there is no advantage in working with the reduced normal form (where γ = 0).

3.1 Magic squares: Cubical count (by upper bound)

Here we count by a strict upper bound t on the permitted values; since the largest entry is α + β + γ, the bound is α + β + γ < t. The number of squares with upper bound t is M_{c}(t). We think of each magic square as a t^{−1}-lattice point in P_{c}^{∘} \ ⋃ ℋ_{c}, the (relative) interior of the inside-out polytope

P_{c} := {(x, y, z) : 0 ≤ y ≤ x, x + y ≤ z, x + y + z ≤ 1},
ℋ_{c} := {h} where h : x = 2y,

but multiplied by t to make the entries integers. Here we use normalized coordinates x = α/t, y = β/t, and z = γ/t. The semilattice of flats is ℒ(P_{c}^{∘}, ℋ_{c}) = {P_{c}^{∘}, h∩P_{c}^{∘}} with P_{c}^{∘} < h∩P_{c}^{∘}. The vertices are

O = (0, 0, 0), C = (0, 0, 1), D = (\frac{1}{2}, 0, \frac{1}{2}), E = (\frac{1}{3}, \frac{1}{6}, \frac{1}{2}), F = (\frac{1}{4}, \frac{1}{4}, \frac{1}{2}),

of which O, C, D, F are the vertices of P_{c} and O, C, E are those of h ∩ P_{c}. (Both these polytopes are simplices.) From Equation (3),

M_{c}(x) = 8E_{P_{c}^{∘},ℋ_{c}}^{∘} (x) = 8[E_{P_{c}^{∘}}(x) − E_{h∩P_{c}^{∘}}(x)]

which we evaluate by LattE and Ehrhart reciprocity, Equation (1):

= 8 [\frac{x^{8}}{(1 − x)^{2}(1 − x^{2})(1 − x^{4})} − \frac{x^{8}}{(1 − x)^{2}(1 − x^{6})}]

= \frac{8x^{10}(2x^{2} + 1)}{(1 − x)^{2}(1 − x^{4})(1 − x^{6})}

= \frac{8x^{10}(2x^{2} + 1)(x^{4} − x^{2} + 1)(x^{11} + x^{10} + ⋯ + x + 1)^{2}(x^{10} + x^{8} + ⋯ + x^{2} + 1)}{(1 − x^{12})^{4}} .

From this generating function we extract the quasipolynomial

M_{c}(t) = \begin{cases}\frac{t^{3}−16t^{2}+76t−96}{6} = \frac{(t−2)(t−6)(t−8)}{6} , & if t ≡ 0, 2, 6, 8 (mod 12); \\ \frac{t^{3}−16t^{2}+73t−58}{6} = \frac{(t−1)(t^{2}−15t+58)}{6} , & if t ≡ 1 (mod 12); \\ \frac{t^{3}−16t^{2}+73t−102}{6} = \frac{(t−3)(t^{2}−13t+34)}{6} , & if t ≡ 3, 11 (mod 12); \\ \frac{t^{3}−16t^{2}+76t−112}{6} = \frac{(t−4)(t^{2}−12t+28)}{6} , & if t ≡ 4, 10 (mod 12); \\ \frac{t^{3}−16t^{2}+73t−90}{6} = \frac{(t−2)(t−5)(t−9)}{6} , & if t ≡ 5, 9 (mod 12); \\ \frac{t^{3}−16t^{2}+73t−70}{6} = \frac{(t−7)(t^{2}−9t+10)}{6} , & if t ≡ 7 (mod 12); \end{cases}  (8)

and the first few nonzero values for t > 0:

\begin{matrix}t & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23 & 24 \\ M_{c}(t) & 8 & 16 & 40 & 64 & 96 & 128 & 184 & 240 & 320 & 400 & 504 & 608 & 744 & 880 & 1056 \\ m_{c}(t) & 1 & 2 & 5 & 8 & 12 & 16 & 23 & 30 & 40 & 50 & 63 & 76 & 93 & 110 & 132 \end{matrix}

The last row is the number of symmetry classes, or normal squares, i.e., M_{c}(t)/8. The rows are sequences \underline{A108576} and \underline{A108577} in the OEIS [15].

The symmetry of the constituents about residue 1 is curious.

The principal constituent is

\frac{t^{3} − 16t^{2} + 76t − 96}{6} = \frac{(t − 2)(t − 6)(t − 8)}{6} .

Its unsigned constant term, 16, is the number of linear orderings of the cells that are induced by magic squares. Thus, up to the symmetries of a magic square, there are just two order types, even allowing arbitrarily large cell values. (The order types are illustrated in Example 3.11 of our general magic article, [6].)

We confirmed these results by direct enumeration, counting the strongly magic squares for t ≤ 60 [7].

Compare this quasipolynomial to the weak quasipolynomial:

\frac{t^{3}−3t^{2}+5t−3}{6} = \frac{(t−1)(t^{2}−2t+3)}{6} , if t is odd;
\frac{t^{3}−3t^{2}+8t−6}{6} = \frac{(t−1)(t^{2}−2t+6)}{6} , if t is even;

with generating function

\frac{(x^{2} + 2x − 1)(2x^{3} − x^{2} + 1)}{(1 − x)^{2}(1 − x^{2})^{2}} .

3.1.1 Reduced magic squares

A more fundamental count than the number of magic squares with an upper bound is the number of reduced squares. Let R_{mc}(t) be the number of 3 × 3 reduced magic squares with maximum cell value t, and r_{mc}(t) the number of reduced symmetry types, or equivalently of normalized reduced squares with maximum t. Then we have the formulas

M_{c}(t) = ∑_{k=0}^{t−1}(t − 1 − k)R_{mc}(k) and m_{c}(t) = ∑_{k=0}^{t−1}(t − 1 − k)r_{mc}(k),

since every reduced square with maximum k gives t−1−k unreduced squares with maximum < t (and positive entries) by adding l to each entry where 1 ≤ l ≤ t − 1 − k. In terms of generating functions,

M_{c}(x) = \frac{x^{2}}{(1 − x)^{2}} R_{mc}(x) and m_{c}(x) = \frac{x^{2}}{(1 − x)^{2}} r_{mc}(x).  (9)

We deduce the generating functions

r_{mc}(x) = \frac{x^{8}(2x^{2} + 1)}{(1 − x^{4})(1 − x^{6})}

and R_{mc}(x) = 8r_{mc}(x), and from the latter the quasipolynomial

R_{mc}(t) = \{ \begin{matrix}\begin{matrix}2t − 16, & if t ≡ 0 \\ 2t − 4, & if t ≡ 2, 10 \\ 2t − 8, & if t ≡ 4, 8 \\ 2t − 12, & if t ≡ 6 \end{matrix} \} & mod 12; \\ 0,  if t is odd \end{matrix}  (10)

(1/8-th of these for r_{mc}(t)) as well as the first few nonzero values:

\begin{matrix}t & 8 & 10 & 12 & 14 & 16 & 18 & 20 & 22 & 24 & 26 & 28 & 30 & 32 & 34 & 36 & 38 & 40 & 42 \\ R_{mc}(t) & 8 & 16 & 8 & 24 & 24 & 24 & 32 & 40 & 32 & 48 & 48 & 48 & 56 & 64 & 56 & 72 & 72 & 72 \\ r_{mc}(t) & 1 & 2 & 1 & 3 & 3 & 3 & 4 & 5 & 4 & 6 & 6 & 6 & 7 & 8 & 7 & 9 & 9 & 9 \end{matrix}

The sequences are \underline{A174256} and \underline{A174257} in the OEIS [15].

The principal constituent is 2t − 16, whose constant term in absolute value, 16, is the number of linear orderings of the cells that are induced by magic squares—necessarily, the same number as with M_{c}(t).

We confirmed the constituents by testing them against the coefficients of the generating function for several periods.

Our way of reasoning, from all squares to reduced squares, is backward; logically, one should count reduced squares and then deduce the ordinary magic square numbers from them via Equation (9). Counting magic squares is not hard enough to require that approach, but in treating semimagic and magilatin squares we follow the logical progression since then reduced squares are much easier to handle.

3.2 Magic squares: Affine count (by magic sum)

The number of magic squares with magic sum t = 3γ is M_{a}(t). We take the normalized coordinates x = α/t and y = β/t. The inside-out polytope is

P_{a} : 0 ≤ y ≤ x, x + y ≤ \frac{1}{3},
ℋ_{a} : {h} where h : x = 2y.

The semilattice of flats is ℒ(P_{a}^{∘}, ℋ_{a}) = {P_{a}^{∘}, h ∩ P_{a}^{∘}} with P_{a}^{∘} < h ∩ P_{a}^{∘}. The vertices are

O = (0, 0), D = (\frac{1}{3}, 0), E = (\frac{2}{9}, \frac{1}{9}), F = (\frac{1}{6}, \frac{1}{6}),

of which O, D, F are the vertices of P_{a} and O, E are the vertices of h ∩ P_{a}. From Equations (1)–(3),

M_{a}(x) = 8E_{P_{a}^{∘},ℋ_{a}}^{∘} (x) = 8[E_{P_{a}^{∘}}(x) − E_{h∩P_{a}^{∘}}(x)]

= 8 [\frac{x^{12}}{(1 − x^{3})^{2}(1 − x^{6})} − \frac{x^{12}}{(1 − x^{3})(1 − x^{9})}]

= \frac{8x^{15}(2x^{3} + 1)}{(1 − x^{3})(1 − x^{6})(1 − x^{9})}

= \frac{8x^{15}(2x^{3} + 1)(x^{9} + 1)(x^{12} + x^{6} + 1)(x^{15} + x^{12} + ⋯ + x^{3} + 1)}{(1 − x^{18})^{3}} .

From this generating function we extract the quasipolynomial

M_{a}(t) = \begin{cases}\frac{2t^{2}−32t+144}{9} = \frac{2}{9}(t^{2} − 16t + 72), & if t ≡ 0 (mod 18); \\ \frac{2t^{2}−32t+78}{9} = \frac{2}{9}(t − 3)(t − 13), & if t ≡ 3 (mod 18); \\ \frac{2t^{2}−32t+120}{9} = \frac{2}{9}(t − 6)(t − 10), & if t ≡ 6 (mod 18); \\ \frac{2t^{2}−32t+126}{9} = \frac{2}{9}(t − 7)(t − 9), & if t ≡ 9 (mod 18); \\ \frac{2t^{2}−32t+96}{9} = \frac{2}{9}(t − 4)(t − 12), & if t ≡ 12 (mod 18); \\ \frac{2t^{2}−32t+102}{9} = \frac{2}{9}(t^{2} − 16t + 51), & if t ≡ 15 (mod 18); \\ 0, & if t ≢ 0 (mod 3); \end{cases}  (11)

and the first few nonzero values for t > 0:

\begin{matrix}t & 15 & 18 & 21 & 24 & 27 & 30 & 33 & 36 & 39 & 42 & 45 & 48 & 51 & 54 \\ M_{a}(t) & 8 & 24 & 32 & 56 & 80 & 104 & 136 & 176 & 208 & 256 & 304 & 352 & 408 & 472 \\ m_{a}(t) & 1 & 3 & 4 & 7 & 10 & 13 & 17 & 22 & 26 & 32 & 38 & 44 & 51 & 59 \end{matrix}

The last row is the number of symmetry classes, or normalized squares, which is M_{a}(t)/8. The two sequences are \underline{A108578} and \underline{A108579} in the OEIS [15].

The principal constituent is

\frac{2t^{2} − 32t + 144}{9} = \frac{2}{9}(t^{2} − 16t + 72),

whose constant term, 16, is the number of linear orderings of the cells that are induced by magic squares—the same number as with M_{c}(t).

We verified our results by direct enumeration, counting the strong magic squares for t ≤ 72 [7].

Compare the magic-square quasipolynomial to the weak quasipolynomial:

\begin{cases}\frac{2t^{2}−6t+9}{9} , & if t ≡ 0 (mod 3); \\ 0, & if t ≢ 0 (mod 3); \end{cases}

due to MacMahon [12, Vol. II, par. 409, p. 163], with generating function

\frac{5x^{6} − 2x^{3} + 1}{(1 − x^{3})^{3}} .

3.2.1 Reduced magic squares

Let R_{ma}(t) be the number of 3 × 3 reduced magic squares with magic sum t, and r_{ma}(t) the number of reduced symmetry types, or equivalently of normalized reduced squares with magic sum t. Then

M_{a}(t) = ∑_{\begin{matrix}0<s<t \\ s≡t (mod 3) \end{matrix}} R_{ma}(s) and m_{a}(t) = ∑_{\begin{matrix}0<s<t \\ s≡t (mod 3) \end{matrix}} r_{ma}(s),

since every reduced square with sum s = t − 3k, where 0 < 3k ≤ t − 3, gives one unreduced square with sum t (and positive entries) by adding 3k to each entry. In terms of generating functions,

m_{a}(x) = \frac{x^{3}}{1 − x^{3}} r_{ma}(x) ;

thus,

r_{ma}(x) = \frac{x^{12}(2x^{3} + 1)}{(1 − x^{6})(1 − x^{9})} ;

and R_{ma}(x) = 8r_{ma}(x). The quasipolynomial is

R_{ma}(t) = 8r_{ma}(t) = \{ \begin{matrix}\begin{matrix}\frac{4}{3}t − 16 = \frac{4}{3}(t − 12), & if t ≡ 0 \\ \frac{4}{3}t − 4 = \frac{4}{3}(t − 3), & if t ≡ 3, 15 \\ \frac{4}{3}t − 8 = \frac{4}{3}(t − 6), & if t ≡ 6, 12 \\ \frac{4}{3}t − 12 = \frac{4}{3}(t − 9), & if t ≡ 9 \end{matrix} \} & mod 18; \\ 0,  if t ≢ 0 & mod 18. \end{matrix}  (12)

The initial nonzero values:

\begin{matrix}t & 12 & 15 & 18 & 21 & 24 & 27 & 30 & 33 & 36 & 39 & 42 & 45 & 48 & 51 & 54 & 57 & 60 & 63 \\ R_{ma}(t) & 8 & 16 & 8 & 24 & 24 & 24 & 32 & 40 & 32 & 48 & 48 & 48 & 56 & 64 & 56 & 72 & 72 & 72 \\ r_{ma}(t) & 1 & 2 & 1 & 3 & 3 & 3 & 4 & 5 & 4 & 6 & 6 & 6 & 7 & 8 & 7 & 9 & 9 & 9 \end{matrix}

These sequences are \underline{A174256} and \underline{A174257} in the OEIS [15].

One of the remarkable properties of magic squares of order 3 is that R_{mc}(2k) = R_{ma}(3k). The reason is that the middle term of a reduced 3 × 3 magic square equals s/3, if s is the magic sum, and the largest entry is 2s/3. Thus, the reduced squares of cubic and affine type, allowing for the difference in parameters, are the same, and although the counts of magic squares by magic sum and by upper bound differ, the only reason is that the reduced squares are adjusted differently to get all squares.

The principal constituent \frac{4}{3}t − 16 has constant term whose absolute value is the same as with all other reduced magic quasipolynomials.

We confirmed the constituents by comparing their values to the coefficients of the generating function for several periods.

4 Semimagic squares of order 3

Now we apply our approach to counting semimagic squares. Here is the general form of a reduced, normalized 3 × 3 semimagic square, in which the magic sum is s = 2α + 2β + γ:

\begin{matrix}0 & β & 2α + β + γ \\ α + β & α + β + γ − δ & δ \\ α + β + γ & α + δ & β − δ \end{matrix}  (13)

Proposition 1. A reduced and normalized 3 × 3 semimagic square has the form (13) with the restrictions

\begin{aligned}0 & < α, β, γ; \\ 0 & < δ < β; \end{aligned}  (14)

and

δ ≠ \begin{cases}\frac{β−α}{2}, \frac{β}{2} , \frac{β+γ}{2} , \frac{β+α+γ}{2} ; \\ β − α, α + γ; \\ γ. \end{cases}  (15)

The largest entry in the square is w := x_{13} = 2α + β + γ.

Each reduced normal square with largest entry w corresponds to exactly 72(t − w − 1) different magic squares with entries in the range (0, t), for 0 < w < t. Each reduced normal square with magic sum s corresponds to exactly 72 different magic squares with magic sum equal to t, if t ≡ s (mod 3), and none otherwise, for 0 < s < t.

Proof. By permuting rows and columns we can arrange that x_{11} = min x_{ij} and that the top row and left column are increasing. By flipping the square over the main diagonal we can further force x_{21} > x_{12}. By subtracting the least entry from every entry we ensure that x_{11} = 0. Thus we account for the 72(t − w − 1) semimagic squares that correspond to each reduced normal square.

The form of the top and left sides in (13) is explained by the fact that x_{11} < x_{12} < x_{21} < x_{31} < x_{13}. The conditions x_{ij} > x_{11} for i, j = 2, 3, together with the row-sum and column-sum equations, imply that x_{13} is the largest entry and that x_{23} < x_{12}.

The only possible equalities amongst the entries are ruled out by the following inequations:

x_{22} ≠ x_{12}, x_{21}, x_{23};
x_{32} ≠ x_{12}, x_{22};
x_{33} ≠ x_{23}, x_{32}.

These correspond to the restrictions (15):

x_{22} ≠ x_{12} ⟺ δ ≠ α + γ;
x_{22} ≠ x_{21} ⟺ δ ≠ γ;
x_{22} ≠ x_{23} ⟺ δ ≠ \frac{β + α + γ}{2} ;
x_{32} ≠ x_{12} ⟺ δ ≠ β − α;
x_{32} ≠ x_{22} ⟺ δ ≠ β + γ − δ ⟺ δ ≠ \frac{β + γ}{2} ;
x_{33} ≠ x_{23} ⟺ δ ≠ β − δ ⟺ δ ≠ \frac{β}{2} ;
x_{33} ≠ x_{32} ⟺ δ ≠ α + δ ≠ β − δ ⟺ δ ≠ \frac{β − α}{2}.  ■

4.1 Semimagic squares: Cubical count (by upper bound)

We are counting squares by a strict upper bound on the allowed value of an entry; this bound is the parameter t. Let S_{c}(t), for t > 0, be the number of semimagic squares of order 3 in which every entry belongs to the range (0, t).

4.1.1 Counting the weak squares

The polytope P_{c} is the 5-dimensional intersection of [0, 1]^{9} with the semimagic subspace in which all row and column sums are equal. This polytope is integral because it is the intersection with an integral polytope of a subspace whose constraint matrix is totally unimodular; so the quasipolynomial is a polynomial. By contrast, the inside-out polytope (P, ℋ) for enumerating strong semimagic squares has denominator 60. We verified by computer counts for t ≤ 18 that the weak polynomial is

\frac{3t^{5} − 15t^{4} + 35t^{3} − 45t^{2} + 32t − 10}{10} = \frac{(t − 1)(t^{2} − 2t + 2)(3t^{2} − 6t + 5)}{10}

with generating function

\frac{(7x^{2} − 2x + 1)(2x^{3} + x^{2} + 4x − 1)}{(1 − x)^{6}} .

We conclude that P_{c} has volume 3/10.

4.1.2 Reduction of the number of strong squares

We compute S_{c}(t) via R_{c}(w), the number of reduced squares with largest entry equal to w. The formula is

S_{c}(t) = ∑_{w=0}^{t−1}(t − 1 − w)R_{c}(w).  (16)

The value R_{c}(w) = 72r_{c}(w) where r_{c}(w) is the number of normal reduced squares (in which we know the largest entry to be x_{13}). Thus r_{c}(s) counts the number of \frac{1}{s}-integral points in the interior of the 3-dimensional polytope Q_{c} defined by

0 ≤ x, y;  0 ≤ z ≤ y;  x + y ≤ \frac{1}{2}  (17)

with the seven excluded (hyper)planes

z = \frac{y − x}{2}, \frac{y}{2}, \frac{1 − y − 2x}{2}, \frac{1 − x − y}{2} , y − x, 1 − x − 2y, 1 − 2x − 2y ,  (18)

the three coordinates being x = α/w, y = β/w, and z = δ/w.

The hyperplane arrangement for reduced normal squares is that of (18). We call it ℐ_{c}. Thus r_{c}(s) = E_{Q_{c}^{∘},ℐ_{c}}^{∘}(s).

4.1.3 Geometrical analysis of the reduced normal polytope

We apply Möbius inversion, Equation (3), over the intersection poset ℒ(Q_{c}^{∘}, ℐ_{c}). We need to know not only ℒ(Q_{c}^{∘}, ℐ_{c}) but also all the vertices of (Q_{c}, ℐ_{c}), since they are required for computing the Ehrhart generating function and estimating the period of the reduced normal quasipolynomial.

We number the planes:

π_{1} : x − y + 2z = 0,
π_{2} :  y − 2z = 0,
π_{3} :  2x + 2z = 1,
π_{4} :  x + 2z = 1,
π_{5} : x − y + z = 0,
π_{6} : x + y + z = 1,
π_{7} : 2x + y + z = 1.

The intersection of two planes, π_{j} ∩ π_{k}, is a line we call l_{jk}; π_{3} ∩ π_{5} ∩ π_{6} is a line we also call l_{356}. The intersection of three planes is, in general, a point but not usually a vertex of (Q_{c}, ℐ_{c}).

Our notation for the line segment with endpoints X, Y is XY, while \overline{XY} denotes the entire line spanned by the points. The triangular convex hull of three noncollinear points X, Y, Z is XY Z. We do not need quadrilaterals, as the intersection of each plane with Q is a triangle.

We need to find the intersections of the planes with Q_{c}^{∘}, separately and in combination. Here is a list of significant points; we shall see it is the list of vertices of (Q_{c}, ℐ_{c}). The first column has the vertices of Q_{c}, the second the vertices of (Q_{c}, ℐ_{c}) that lie in open edges, the third the vertices that lie in open facets, and the last is the sole interior vertex.

O = (0, 0, 0), D_{c} = (0, \frac{1}{2}, \frac{1}{2}) ∈ OC,  F_{c} = (0, \frac{2}{3}, \frac{1}{3}) ∈ OBC, H_{c} = (\frac{1}{5}, \frac{2}{5}, \frac{1}{5}).
A = (\frac{1}{2}, 0, 0), E_{c} = (\frac{1}{3}, \frac{1}{3}, 0) ∈ AB, G_{c} = (\frac{1}{4}, \frac{1}{2}, \frac{1}{4}) ∈ ABC,
B_{c} = (0, 1, 0), E′_{c} = (\frac{1}{3}, \frac{1}{3}, \frac{1}{3}) ∈ AC, G′_{c} = (\frac{1}{5}, \frac{3}{5}, \frac{1}{5}) ∈ ABC,
C_{c} = (0, 1, 1), E′′_{c} = (0, 1, \frac{1}{2}) ∈ BC, G′′_{c} = (\frac{1}{5}, \frac{3}{5}, \frac{2}{5}) ∈ ABC,

(19)

The denominator of (Q_{c}, ℐ_{c}) is the least common denominator of all the points; it evidently equals 4 · 3 · 5 = 60.

The intersections of the planes with the edges of Q_{c} are in Table 1. The subscript c is omitted.

Table 2 shows the lines generated by pairwise intersection of planes.

Plane  Intersection with edge  Intersection
OA OB OC AB AC BC  with Q
π_{1}  O O O E ∉ Q E′′  OEE′′
π_{2}  OA O O A A E′′  OAE′′
π_{3}  A ∅ D A A E′′  ADE′′
π_{4}  ∉ Q ∅ D ∉ Q E′ E′′  DE′′E′
π_{5}  O O OC E  C C  OCE
π_{6}  ∉ Q B D B E′ B  BDE′
π_{7}  A B D AB A B  ABD

Table 1: Intersections of planes of ℐ with edges of Q. A vertex of Q contained in π_{j} will show up three times in the row of π_{j}. In order to clarify the geometry, we distinguish between a plane’s meeting an edge line outside Q and not meeting it at all (i.e., their being parallel).

π_{2}  π_{3}  π_{4}  π_{5}  π_{6}  π_{7}
π_{1}  \begin{aligned}x & = 0 \\ y & = 2z \end{aligned}  \begin{aligned}x + z & = \frac{1}{2} \\ y − z & = \frac{1}{2} \end{aligned}  \begin{aligned}y & = 1 \\ x + 2z & = 1 \end{aligned}  \begin{aligned}z & = 0 \\ x & = y \end{aligned}  \begin{aligned}z & = 2y − 1 \\ x & = 2 − 3y \end{aligned}  \begin{aligned}x + z & = \frac{1}{3} \\ y − z & = \frac{1}{3} \end{aligned}
π_{2}  \begin{aligned}y & = 2z \\ x + z & = \frac{1}{2} \end{aligned}  \begin{aligned}y & = 2z \\ x + y & = 1 \end{aligned}  \begin{aligned}y & = 2z \\ x & = z \end{aligned}  \begin{aligned}y & = 2z \\ x + 3z & = 1 \end{aligned}  \begin{aligned}y & = 2z \\ 2x + 3z & = 1 \end{aligned}
π_{3}  \begin{aligned}x & = 0 \\ z & = \frac{1}{2} \end{aligned}  l_{356}:  \begin{aligned}x + z & = \frac{1}{2} \\ y & = \frac{1}{2} \end{aligned}  \begin{aligned}x + z & = \frac{1}{2} \\ y & = z \end{aligned}
π_{4}  \begin{aligned}x & = 1 − 2z \\ y & = 1 − z \end{aligned}  \begin{aligned}x & = 1 − 2z \\ y & = z \end{aligned}  \begin{aligned}x & = 1 − 2z \\ y & = 3z − 1 \end{aligned}
π_{5}  l_{356}  \begin{aligned}x & = 1 − 2y \\ z & = 3y − 1 \end{aligned}
π_{6}  \begin{aligned}x & = 0 \\ y + z & = 1 \end{aligned}

Table 2: The equations of the pairwise intersections of planes of ℐ_{c}.

In Table 3 we describe the intersection of each line with Q_{c} and with its interior. The subscript c is omitted.

π_{2}  π_{3}  π_{4}  π_{5}  π_{6}  π_{7}
π_{1}  \begin{matrix}OE′′ \\ (OBC) \end{matrix}  \begin{matrix}E′′ \\ (BC) \end{matrix}  \begin{matrix}E′′ \\ (BC) \end{matrix}  \begin{matrix}OE \\ (OAB) \end{matrix}  F G′  EF
π_{2}  \begin{matrix}AE′′ \\ (ABC) \end{matrix}  \begin{matrix}E′′ \\ (BC) \end{matrix}  OG  F G  AF
π_{3}  \begin{matrix}DE′′ \\ (OBC) \end{matrix}  l_{356}: DG  \begin{matrix}AD \\ (OAC) \end{matrix}
π_{4}  DG′′  \begin{matrix}DE′ \\ (OAC) \end{matrix}  \begin{matrix}D \\ (OC) \end{matrix}
π_{5}  l_{356}  ED
π_{6}  \begin{matrix}BD \\ (OBC) \end{matrix}

Table 3: The intersections of lines with Q and Q^{∘}. The second (parenthesized) row in each box shows the smallest face of Q to which the intersection belongs, if that is not Q itself; these intersections are not part of the intersection poset of (Q^{∘}, ℐ).

Last, we need the intersection points of three planes of ℐ_{c}; or, of a plane and a line. Some are not in Q_{c} at all; them we can ignore. Some are on the boundary of Q_{c}; they are necessary in finding the denominator, but all of them are points already listed in (19). It turns out that

π_{2} ∩ π_{5} ∩ π_{7} = H_{c}

is the only vertex in Q_{c}^{∘}, so it is the only one we need for the intersection poset.

Here, then, is the intersection poset (Figure 1). The subscript c is omitted. In the figure, for simplicity, we write π_{j}, etc., when the actual element is the simplex π_{j} ∩ Q^{∘}, etc.; we also state the vertices of the simplex. The Möbius function μ(ˆ0, u) equals (−1)^{codim u} with the exception of μ(ˆ0, l_{356}) = 2.

4.1.4 Generating functions and the quasipolynomial

That was the first half of the work. The second half begins with finding r_{c}(w) = E_{(Q_{c}^{∘},ℐ_{c})}^{∘}(w) from the Ehrhart generating functions E_{u}(t) of the intersections by means of (6). The next step, then, is to calculate all necessary generating functions. This is done by LattE. r_{c}(x) is the result of applying reciprocity to the sum of all these rational functions (with the appropriate Möbius-function multiplier −1, excepting E_{DG^{∘}}(x) whose multiplier is −2). The PSfrag replacements

ℒ(Q^{∘}, ℐ)

ℝ^{3} (OABC) π_{1} (OEE′′) π_{2} (OAE′′) π_{3} (ADE′′) π_{4} (DE′′E′) π_{5} (OCE) π_{6} (BDE′) π_{7} (ABD) l_{16} (F G′) l_{17} (EF) l_{356} (DG) l_{57} (ED) l_{25} (OG) l_{26} (F G) l_{27} (AF) l_{45} (DG′′)

H Figure 1: The intersection poset ℒ(Q^{∘}, ℐ) for semimagic squares. The diagram shows both the flats and (in parentheses) their intersections with Q. result is:

−r_{c}(1/x) = E_{OABC}(x) + E_{OEE′′}(x) + E_{OAE′′}(x)

+ E_{ADE′′}(x) + E_{DE′E′′}(x) + E_{OCE}(x) + E_{BDE′}(x) + E_{ABD}(x) + E_{F G′}(x) + E_{EF} (x) + E_{OG}(x) + E_{F G}(x) + E_{AF} (x) + 2E_{DG}(x) + E_{DG′′}(x) + E_{DE}(x) + E_{H}(x)

= \frac{1}{(1 − x)^{3}(1 − x^{2})} + \frac{1}{(1 − x)(1 − x^{2})(1 − x^{3})} + \frac{1}{(1 − x)(1 − x^{2})^{2}}

+ \frac{1}{(1 − x^{2})^{3}} + \frac{1}{(1 − x^{2})^{2}(1 − x^{3})} + \frac{1}{(1 − x)^{2}(1 − x^{3})}

+ \frac{1}{(1 − x)(1 − x^{2})(1 − x^{3})} + \frac{1}{(1 − x)(1 − x^{2})^{2}} + \frac{1}{(1 − x^{3})(1 − x^{5})}

+ \frac{1}{(1 − x^{3})^{2}} + \frac{1}{(1 − x)(1 − x^{4})} + \frac{1}{(1 − x^{3})(1 − x^{4})}

+ \frac{1}{(1 − x^{2})(1 − x^{3})} + 2\frac{1}{(1 − x^{2})(1 − x^{4})} + \frac{1}{(1 − x^{2})(1 − x^{5})}

+ \frac{1}{(1 − x^{2})(1 − x^{3})} + \frac{1}{1 − x^{5}} .

(20)

Then by (16) the generating function for cubically counted semimagic squares is

\begin{aligned}S_{c}(x) & = 72\frac{x^{2}}{(1 − x)^{2}} r_{c}(x) \\ & = \frac{72x^{10} [18 x^{9} + 46 x^{8} + 69 x^{7} + 74 x^{6} + 65 x^{5} + 46 x^{4} + 26 x^{3} + 11 x^{2} + 4 x + 1]}{(1 − x^{2})^{2}(1 − x^{3})^{2}(1 − x^{4})(1 − x^{5})} . \end{aligned}

(21)

From the geometrical denominator 60 or the (standard-form) algebraic denominator (1− x^{60})^{5} we know the period of S_{c}(t) divides 60. We compute the constituents of S_{c}(t) by the method of Section 2.1; the result is that

S_{c}(t) = \begin{cases}\frac{3}{10}t^{5} − \frac{75}{8}t^{4} + \frac{331}{3}t^{3} − \frac{5989}{10}t^{2} + c_{1}(t)t − c_{0}(t), & if t is even; \\ \frac{3}{10}t^{5} − \frac{75}{8}t^{4} + \frac{331}{3}t^{3} − \frac{11933}{20}t^{2} + c_{1}(t)t − c_{0}(t), & if t is odd; \end{cases}  (22)

where c_{1} varies with period 6, given by

c_{1}(t) = \begin{cases}1464, & if t ≡ 0, 2 \\ 1456, & if t ≡ 4 \\ \frac{2831}{2} , & if t ≡ 1 \\ \frac{2847}{2} , & if t ≡ 3, 5 \end{cases} (mod 6);

and c_{0}, given by Table 4, varies with period 60. (It is curious that the even constant terms have half the period of the odd terms.) Thus the period of S_{c} turns out to be 60, the largest it could be.

\begin{matrix}t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) \\ 0 & 1296 & 12 & 1296 & 24 & \frac{6624}{5} & 36 & \frac{6192}{5} & 48 & \frac{6624}{5} \\ 1 & \frac{110413}{120} & 13 & \frac{120781}{120} & 25 & \frac{23465}{24} & 37 & \frac{23465}{24} & 49 & \frac{120781}{120} \\ 2 & \frac{3824}{3} & 14 & \frac{19552}{15} & 26 & \frac{18256}{15} & 38 & \frac{19552}{15} & 50 & \frac{3824}{3} \\ 3 & \frac{47727}{40} & 15 & \frac{9315}{8} & 27 & \frac{9315}{8} & 39 & \frac{47727}{40} & 51 & \frac{44271}{40} \\ 4 & \frac{18152}{15} & 16 & \frac{16856}{15} & 28 & \frac{18152}{15} & 40 & \frac{3544}{3} & 52 & \frac{3544}{3} \\ 5 & \frac{25705}{24} & 17 & \frac{25705}{24} & 29 & \frac{131981}{120} & 41 & \frac{121613}{120} & 53 & \frac{131981}{120} \\ 6 & \frac{6192}{5} & 18 & \frac{6624}{5} & 30 & 1296 & 42 & 1296 & 54 & \frac{6624}{5} \\ 7 & \frac{25193}{24} & 19 & \frac{129421}{120} & 31 & \frac{119053}{120} & 43 & \frac{129421}{120} & 55 & \frac{25193}{24} \\ 8 & \frac{19552}{15} & 20 & \frac{3824}{3} & 32 & \frac{3824}{3} & 44 & \frac{19552}{15} & 56 & \frac{18256}{15} \\ 9 & \frac{44847}{40} & 21 & \frac{41391}{40} & 33 & \frac{44847}{40} & 45 & \frac{8739}{8} & 57 & \frac{8739}{8} \\ 10 & \frac{3544}{3} & 22 & \frac{3544}{3} & 34 & \frac{18152}{15} & 46 & \frac{16856}{15} & 58 & \frac{18152}{15} \\ 11 & \frac{130253}{120} & 23 & \frac{140621}{120} & 35 & \frac{27433}{24} & 47 & \frac{27433}{24} & 59 & \frac{140621}{120} \end{matrix}

Table 4: Constant terms (without the negative sign) of the constituents of the semimagic cubical quasipolynomial S_{c}(t).

The principal constituent (for t ≡ 0) is

\frac{3}{10}t^{5} − \frac{75}{8}t^{4} + \frac{331}{3}t^{3} − \frac{5989}{10}t^{2} + 1464t − 1296.

Its unsigned constant term, 1296, is the number of order types of semimagic squares. Allowing for the 72 symmetries of a semimagic square, there are just 18 symmetry classes of order types.

For the first few nonzero values of S_{c}(t) see the following table. (This is sequence \underline{A173546} in the OEIS [15].) The third row is the number of normalized squares, or symmetry classes (sequence \underline{A173723}), which equals S_{c}(t)/72. The other lines give the numbers of reduced squares (sequence \underline{A173727}) and of reduced normal squares (i.e., symmetry types of reduced squares; sequence \underline{A173724}), which may be of interest.

\begin{matrix}t & 8 & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 \\ S_{c}(t) & 0 & 0 & 72 & 288 & 936 & 2592 & 5760 & 11520 & 20952 & 35712 & 57168 & 88272 \\ s_{c}(t) & 0 & 0 & 1 & 4 & 13 & 36 & 80 & 160 & 291 & 496 & 794 & 1226 \\ R_{c}(t) & 72 & 144 & 432 & 1008 & 1512 & 2592 & 3672 & 5328 & 6696 & 9648 & 11736 & 15552 \\ r_{c}(t) & 1 & 2 & 6 & 14 & 21 & 36 & 51 & 74 & 93 & 134 & 163 & 216 \end{matrix}

Compare the strong to the weak quasipolynomial. The leading coefficients agree and the strong coefficient of t^{4} is constant. These facts, of which we made no use in deducing the quasipolynomial, provide additional verification of the correctness of the counts and constituents.

4.1.5 Another method: Direct counting

We checked the constituents by directly counting (in Maple) all semimagic squares for t ≤ 100. The numbers agreed with those derived from the generating function and quasipolynomial above.

4.2 Semimagic squares: Affine count (by magic sum)

Now we count squares by magic sum: we compute S_{a}(t), the number of squares with magic sum t.

4.2.1 The Birkhoff polytope

The polytope P for semimagic squares of order 3, counted by magic sum, is 4-dimensional and integral. (It is the polytope of doubly stochastic matrices of order 3, i.e., a Birkhoff polytope [8, 4].)

4.2.2 Affine weak semimagic

The polytope for weak semimagic squares of order 3 is the same P.

The weak quasipolynomial, or rather, polynomial, first computed by MacMahon [12, Vol. II, par. 407, p. 161], is

\frac{t^{4} − 6t^{3} + 15t^{2} − 18t + 8}{8} = \frac{(t − 1)(t − 2)(t^{2} − 3t + 4)}{8}

with generating function

\frac{6x^{4} − 9x^{3} + 10x^{2} − 5x + 1}{(1 − x)^{5}} .

4.2.3 Reduction

The count is via R_{a}(s), the number of reduced squares with magic sum s. The formula is

S_{a}(t) = ∑_{\begin{matrix}0<s≤t−3 \\ s≡t (mod 3) \end{matrix}} R_{a}(s)  if t > 0.  (23)

We have R_{a}(s) = 72r_{a}(w), where r_{a}(s) is the number of reduced, normalized squares with magic sum s, equivalently the number of \frac{1}{s}-integral points in the interior of the 3-dimensional polytope Q_{a} defined by

0 ≤ x, y;  0 ≤ z ≤ y;  x + y ≤ \frac{1}{2}  (24)

with the seven excluded (hyper)planes

z = \frac{y − x}{2}, \frac{y}{2}, \frac{1 − y − 2x}{2}, \frac{1 − x − y}{2} , y − x, 1 − x − 2y, 1 − 2x − 2y ,  (25)

the three coordinates being x = α/s, y = β/s, and z = δ/s.

The hyperplane arrangement for reduced, normalized squares is that of (25). We call it ℐ_{a}. Thus r_{a}(s) = E_{Q_{a}^{∘},ℐ_{a}}^{∘}(s).

4.2.4 The reduced, normalized weak polytopal quasipolynomial

This function simply counts \frac{1}{s}-lattice points in Q_{a}^{∘}. The counting formula is ∑_{α} ∑_{β} ∑_{δ} 1, summed over all triples that satisfy (14). It simplifies to

∑_{α} \binom{⌊\frac{s−1}{2}⌋ − α}{2},

which gives the Ehrhart quasipolynomial

E_{Q_{a}^{∘}}(s) = \binom{⌊\frac{s−1}{2}⌋}{3} = \begin{cases}\binom{\frac{s−1}{2}}{3} = \frac{1}{48}(s − 1)(s − 3)(s − 5), & for odd s; \\ \binom{\frac{s−2}{2}}{3} = \frac{1}{48}(s − 2)(s − 4)(s − 6), & for even s. \end{cases}  (26)

The leading coefficient is vol Q_{a}.

We deduce from (26) that

E_{P^{∘}}(x) = \frac{x^{7}(1 + x)}{(1 − x^{2})^{4}}

and by reciprocity that

E_{P} (x) = \frac{1 + x}{(1 − x^{2})^{4}} .

4.2.5 Geometrical analysis of the reduced, normalized polytope

We apply Möbius inversion, Equation (6), over the intersection poset ℒ(Q_{a}^{∘}, ℐ). We number the planes:

π_{1} : x − y + 2z = 0,
π_{2} :  y − 2z = 0,
π_{3} : 2x + y + 2z = 1,
π_{4} : x + y + 2z = 1,
π_{5} : x − y + z = 0,
π_{6} : x + 2y + z = 1,
π_{7} : 2x + 2y + z = 1.

The intersection of two planes, π_{j} ∩ π_{k}, is a line we call l_{jk}; π_{3} ∩ π_{5} ∩ π_{6} is a line we also call l_{356}. The intersection of three planes is, in general, a point but not usually a vertex of (Q_{a}, ℐ_{a}). Our geometrical notation is as in the cubical analysis.

We need to find the intersections of the planes with Q_{a}^{∘}, separately and in combination. Here is a list of significant points; we shall see it is the list of vertices of (Q_{a}, ℐ_{a}). The first column has the vertices of Q_{a}, the second the vertices of (Q_{a}, ℐ_{a}) that lie in open edges, the third the vertices that lie in open facets, and the last is the sole interior vertex.

O = (0, 0, 0), D_{a} = (0, \frac{1}{3}, \frac{1}{3}) ∈ OC, F_{a} = (0, \frac{2}{5}, \frac{1}{5}) ∈ OBC, H_{a} = (\frac{1}{7}, \frac{2}{7}, \frac{1}{7}),
A = (\frac{1}{2}, 0, 0), E_{a} = (\frac{1}{4}, \frac{1}{4}, 0) ∈ AB, G_{a} = (\frac{1}{6}, \frac{2}{6}, \frac{1}{6}) ∈ ABC,
B_{a} = (0, \frac{1}{2}, 0), E′_{a} = (\frac{1}{4}, \frac{1}{4}, \frac{1}{4}) ∈ AC, G′_{a} = (\frac{1}{8}, \frac{3}{8}, \frac{1}{8}) ∈ ABC,
C_{a} = (0, \frac{1}{2}, \frac{1}{2}), E′′_{a} = (0, \frac{1}{2}, \frac{1}{4}) ∈ BC, G′′_{a} = (\frac{1}{8}, \frac{3}{8}, \frac{2}{8}) ∈ ABC.

(27)

The least common denominator of O, A, B_{a}, C_{a} explains the period 2 of E_{Q_{a}^{∘}}. The denominator of (Q_{a}, ℐ_{a}) is the least common denominator of all the points; it evidently equals 8 · 3 · 5 · 7 = 840.

The intersections of the planes with the edges of Q_{a} are in Table 1. Table 5 shows the lines generated by pairwise intersection of planes. Table 3 describes the intersection of each line with Q_{a} and with Q_{a}^{∘}.

π_{2}  π_{3}  π_{4}  π_{5}  π_{6}  π_{7}
π_{1}  \begin{aligned}x & = 0 \\ y & = 2z \end{aligned}  \begin{aligned}x & = \frac{1−4z}{3} \\ y & = \frac{1+2z}{3} \end{aligned}  \begin{aligned}x + 2z & = \frac{1}{2} \\ y & = \frac{1}{2} \end{aligned}  \begin{aligned}x & = y \\ z & = 0 \end{aligned}  \begin{aligned}x & = 2 − 5y \\ z & = 3y − 1 \end{aligned}  \begin{aligned}x & = \frac{1−5z}{4} \\ y & = \frac{1+3z}{4} \end{aligned}

π_{2}  \begin{aligned}x + y & = \frac{1}{2} \\ y & = 2z \end{aligned}  \begin{aligned}x & = 1 − 2y \\ y & = 2z \end{aligned}  \begin{aligned}x & = z \\ y & = 2z \end{aligned}  \begin{aligned}x & = 1 − 5z \\ y & = 2z \end{aligned}  \begin{aligned}x & = \frac{1−5z}{2} \\ y & = 2z \end{aligned}

π_{3}  \begin{aligned}x & = 0 \\ y & = 1 − 2z \end{aligned}  l_{356}:  \begin{aligned}x + z & = \frac{1}{3} \\ y & = \frac{1}{3} \end{aligned}  \begin{aligned}2x + 3y & = 1 \\ z & = y \end{aligned}

π_{4}  \begin{aligned}x & = 3y − 1 \\ z & = 1 − 2y \end{aligned}  \begin{aligned}x & = 1 − 3y \\ z & = y \end{aligned}  \begin{aligned}x + y & = \frac{1}{3} \\ z & = \frac{1}{3} \end{aligned}

π_{5}  l_{356}  \begin{aligned}x & = 1 − 3y \\ z & = 4y − 1 \end{aligned}

π_{6}  \begin{aligned}x & = 0 \\ z & = 1 − 2y \end{aligned}

Table 5: The equations of the pairwise intersections of planes of ℐ_{a}.

Last, we need the intersection points of three planes of ℐ_{a}; or, of a plane and a line. Some are not in Q_{a} at all; them we can ignore. Some are on the boundary of Q_{a}; they are necessary in finding the denominator, but all of them are points already listed in (27). It turns out that

π_{2} ∩ π_{5} ∩ π_{7} = H_{a}

is the only vertex in Q_{a}^{∘}, so it is the only one we need for the intersection poset.

The combinatorial structure and the intersection poset (Figure 1) for the affine count are identical to those for the cubical count. The reason is that the affine polytope P_{a} is the 4-dimensional section of P_{c} by the flat in which the magic sum equals 1, and this flat is orthogonal to the line of intersection of the whole arrangement ℋ_{a}.

4.2.6 Generating functions and the quasipolynomial

The second half of the affine solution is to find r_{a}(s) = E_{(Q_{a}^{∘},ℐ_{a})}^{∘}(s) by applying Equations (1)– (4) after finding the Ehrhart generating functions E_{u}(s) for u ∈ ℒ(Q_{a}^{∘}, ℐ_{a}). The next step, then, is to calculate those generating functions. This is done by LattE. Then (−1)^{3}r_{a}(x^{−1}) is the sum of all these rational functions; that is,

−r_{a}(1/x) = E_{OABC}(x) + E_{OEE′′}(x) + E_{OAE′′}(x)

+ E_{ADE′′}(x) + E_{DE′′E′}(x) + E_{OCE}(x) + E_{BDE′}(x) + E_{ABD}(x) + E_{F G′}(x) + E_{EF} (x) + E_{OG}(x) + E_{F G}(x) + E_{AF} (x) + 2E_{DG}(x) + E_{DG′′}(x) + E_{DE}(x) + E_{H}(x)

= \frac{1}{(1 − x)(1 − x^{2})^{3}} + \frac{1}{(1 − x)(1 − x^{4})^{2}} + \frac{1}{(1 − x)(1 − x^{2})(1 − x^{4})}

+ \frac{1}{(1 − x^{2})(1 − x^{3})(1 − x^{4})} + \frac{1}{(1 − x^{3})(1 − x^{4})^{2}} + \frac{1}{(1 − x)(1 − x^{2})(1 − x^{4})}

+ \frac{1}{(1 − x^{2})(1 − x^{3})(1 − x^{4})} + \frac{1}{(1 − x^{2})^{2}(1 − x^{3})} + \frac{1}{(1 − x^{5})(1 − x^{8})}

+ \frac{1}{(1 − x^{4})(1 − x^{5})} + \frac{1}{(1 − x)(1 − x^{6})} + \frac{1}{(1 − x^{5})(1 − x^{6})}

+ \frac{1}{(1 − x^{2})(1 − x^{5})} + 2\frac{1}{(1 − x^{3})(1 − x^{6})} + \frac{1}{(1 − x^{3})(1 − x^{8})}

+ \frac{1}{(1 − x^{3})(1 − x^{4})} + \frac{1}{1 − x^{7}} .

(28)

The generating function for the affine count of semimagic squares, by (23), is

\begin{matrix}S_{a}(x) = 72 \frac{x^{3}}{1 − x^{3}} r_{a}(x) = \\ \frac{72x^{15} \left\{{18x^{21} + 5x^{20} + 15x^{19} + 11x^{17} − 8x^{16} + x^{15} − 23x^{14} − 13x^{13} − 22x^{12} − 9x^{11} \atop − 16x^{10} + x^{9} − 3x^{8} + 7x^{7} + 7x^{6} + 9x^{5} + 7x^{4} + 6x^{3} + 4x^{2} + 2x + 1}\right\}}{(1 − x^{3})^{2}(1 − x^{4})(1 − x^{5})(1 − x^{6})(1 − x^{7})(1 − x^{8})} . \end{matrix}

(29)

From the geometrical or generating-function denominator we know that the period of S_{a}(t) divides 840 = lcm(3, 4, 6, 7, 8). This is long, but it can be simplified. The factor 7 in the period is due to a single term in (28). If we treat it separately we have r_{a} as a sum of the H-term x^{7}/(1 − x^{7}) and a “truncated” generating function for r_{a}(x) + \frac{x^{7}}{1−x^{7}} , and a corresponding truncated expression

S_{a}(x) − 72\frac{x^{10}}{(1 − x^{3})(1 − x^{7})} =

\frac{−72x^{10} \left\{{17x^{19} + 5x^{18} + 15x^{17} + x^{16} + 12x^{15} − 7x^{14} + 2x^{13} − 7x^{12} \atop − 8x^{11} − 9x^{10} − 9x^{9} − 6x^{8} − 6x^{7} − x^{6} + x^{4} + x^{3} − 1}\right\}}{(1 − x^{3})^{2}(1 − x^{4})(1 − x^{5})(1 − x^{6})(1 − x^{8})} .

We extract the constituents from this expression as in Section 2.1, separately for the two parts of the generating function. The constituents are all of the form

S_{a}(t) = \frac{1}{8}t^{4} − \frac{9}{2}t^{3} + a_{2}(t)t^{2} − a_{1}(t)t + a_{0}(t) − 72S_{7}(t),  (30)

where S_{7}(t) is a correction, to be defined in a moment, and

a_{2}(t) = \begin{cases}\frac{243}{4} , & if t ≡ 0 \\ \frac{218}{4} , & if t ≡ 1, 5 \\ \frac{227}{4} , & if t ≡ 2, 4 \\ \frac{234}{4} , & if t ≡ 3 \end{cases} (mod 6);

a_{1}(t) = \begin{cases}\frac{1968}{5} , & if t ≡ 0 \\ \frac{1158}{5} , & if t ≡ 1, 5 \\ \frac{1383}{5} , & if t ≡ 2, 10 \\ \frac{1653}{5} , & if t ≡ 3 \\ \frac{1428}{5} , & if t ≡ 4, 8 \\ \frac{1923}{5} , & if t ≡ 6 \\ \frac{1113}{5} , & if t ≡ 7, 11 \\ \frac{1698}{5} , & if t ≡ 9 \end{cases} (mod 12);

and a_{0}(t) is given in Table 6.

We call the constituents of the quasipolynomial

S_{a}(t) + 72S_{7}(t) = \frac{1}{8}t^{4} − \frac{9}{2}t^{3} + a_{2}(t)t^{2} − a_{1}(t)t + a_{0}(t)

the truncated constituents of S_{a}(t), since they correspond to the truncated generating function mentioned just above. The S_{7} term that undoes the truncation is

S_{7}(t) := ⌊\frac{t − 1}{21}⌋ + \begin{cases}1, & if t ≡ 10, 13, 16, 17, 19, 20 (mod 21); \\ 0, & otherwise \end{cases}

= \frac{t − t̄}{21} + s_{7}(t),

where t̄:= the least positive residue of t modulo 7 and

s_{7}(t) := \begin{cases}1, & if t ≡ 10, 13, 16, 17, 19, 20 (mod 21); \\ 0, & otherwise. \end{cases}

Note that t̄= 21 if t ≡ 0, so that S_{7}(0) = −1 and in general S_{7}(21k) = k − 1.

\begin{matrix}t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) \\ 0 & 1224 & 20 & 524 & 40 & 584 & 60 & 1188 & 80 & 560 & 100 & 548 \\ 1 & \frac{7259}{40} & 21 & \frac{31419}{40} & 41 & \frac{6299}{40} & 61 & \frac{8699}{40} & 81 & \frac{29979}{40} & 101 & \frac{7739}{40} \\ 2 & \frac{1801}{5} & 22 & \frac{1741}{5} & 42 & \frac{5121}{5} & 62 & \frac{1621}{5} & 82 & \frac{1921}{5} & 102 & \frac{4941}{5} \\ 3 & \frac{23067}{40} & 23 & \frac{827}{40} & 43 & \frac{347}{40} & 63 & \frac{24507}{40} & 83 & \frac{−613}{40} & 103 & \frac{1787}{40} \\ 4 & \frac{2452}{5} & 24 & \frac{5832}{5} & 44 & \frac{2332}{5} & 64 & \frac{2632}{5} & 84 & \frac{5652}{5} & 104 & \frac{2512}{5} \\ 5 & \frac{2239}{8} & 25 & \frac{2143}{8} & 45 & \frac{6975}{8} & 65 & \frac{1951}{8} & 85 & \frac{2431}{8} & 105 & \frac{6687}{8} \\ 6 & \frac{4653}{5} & 26 & \frac{1513}{5} & 46 & \frac{1453}{5} & 66 & \frac{4833}{5} & 86 & \frac{1333}{5} & 106 & \frac{1633}{5} \\ 7 & \frac{5243}{40} & 27 & \frac{26523}{40} & 47 & \frac{4283}{40} & 67 & \frac{3803}{40} & 87 & \frac{27963}{40} & 107 & \frac{2843}{40} \\ 8 & \frac{2224}{5} & 28 & \frac{2164}{5} & 48 & \frac{5544}{5} & 68 & \frac{2044}{5} & 88 & \frac{2344}{5} & 108 & \frac{5364}{5} \\ 9 & \frac{31131}{40} & 29 & \frac{8891}{40} & 49 & \frac{8411}{40} & 69 & \frac{32571}{40} & 89 & \frac{7451}{40} & 109 & \frac{9851}{40} \\ 10 & 413 & 30 & 1017 & 50 & 389 & 70 & 377 & 90 & 1053 & 110 & 353 \\ 11 & \frac{539}{40} & 31 & \frac{2939}{40} & 51 & \frac{24219}{40} & 71 & \frac{1979}{40} & 91 & \frac{1499}{40} & 111 & \frac{25659}{40} \\ 12 & \frac{5796}{5} & 32 & \frac{2656}{5} & 52 & \frac{2596}{5} & 72 & \frac{5976}{5} & 92 & \frac{2476}{5} & 112 & \frac{2776}{5} \\ 13 & \frac{7547}{40} & 33 & \frac{28827}{40} & 53 & \frac{6587}{40} & 73 & \frac{6107}{40} & 93 & \frac{30267}{40} & 113 & \frac{5147}{40} \\ 14 & \frac{1477}{5} & 34 & \frac{1777}{5} & 54 & \frac{4797}{5} & 74 & \frac{1657}{5} & 94 & \frac{1597}{5} & 114 & \frac{4977}{5} \\ 15 & \frac{5823}{8} & 35 & \frac{799}{8} & 55 & \frac{1279}{8} & 75 & \frac{5535}{8} & 95 & \frac{1087}{8} & 115 & \frac{991}{8} \\ 16 & \frac{2488}{5} & 36 & \frac{5508}{5} & 56 & \frac{2368}{5} & 76 & \frac{2308}{5} & 96 & \frac{5688}{5} & 116 & \frac{2188}{5} \\ 17 & \frac{8603}{40} & 37 & \frac{11003}{40} & 57 & \frac{32283}{40} & 77 & \frac{10043}{40} & 97 & \frac{9563}{40} & 117 & \frac{33723}{40} \\ 18 & \frac{4689}{5} & 38 & \frac{1189}{5} & 58 & \frac{1489}{5} & 78 & \frac{4509}{5} & 98 & \frac{1369}{5} & 118 & \frac{1309}{5} \\ 19 & \frac{2651}{40} & 39 & \frac{26811}{40} & 59 & \frac{1691}{40} & 79 & \frac{4091}{40} & 99 & \frac{25371}{40} & 119 & \frac{3131}{40} \end{matrix}

Table 6: Constant terms of the truncated constituents of S_{a}(t).

The period of the constant term of the truncated constituents is 120. It follows that S_{a}(t) has period 7 · 840, that is, 840.

The principal constituent of S_{a}(t) (that is, for t ≡ 0) is

\frac{1}{8}t^{4} − \frac{9}{2}t^{3} + \frac{243}{4}t^{2} − \frac{13896}{35}t + 1296.

(This incorporates the effect of the term −72S_{7}.) The constant term is the same as in the cubic count, as it is the number of order types of semimagic squares.

We give the first few nonzero values of S_{a}(t) in the following table. (This sequence is \underline{A173547} in the OEIS [15].) The third row is the number of normalized squares, or symmetry classes (sequence \underline{A173725}); this is S_{a}(t)/72. The last rows are the numbers of reduced squares (sequence \underline{A173728}) and of reduced, normalized squares (sequence \underline{A173726}) with magic sum t.

\begin{matrix}t & 12 & 13 & 14 & 15 & 16 & 17 & 18 & 19 & 20 & 21 & 22 & 23 & 24 \\ S_{a}(t) & 0 & 0 & 0 & 72 & 144 & 288 & 576 & 864 & 1440 & 2088 & 3024 & 3888 & 5904 \\ s_{a}(t) & 0 & 0 & 0 & 1 & 2 & 4 & 8 & 12 & 20 & 29 & 42 & 54 & 82 \\ R_{a}(t) & 72 & 144 & 288 & 504 & 720 & 1152 & 1512 & 2160 & 2448 & 3816 & 3960 & 5544 & 6264 \\ r_{a}(t) & 1 & 2 & 4 & 7 & 10 & 16 & 21 & 30 & 34 & 53 & 55 & 77 & 87 \end{matrix}

4.2.7 Alternative methods: Direct counting and direct computation

We verified our formulas by computing S_{a}(t) for t ≤ 100 through direct enumeration of normal squares. The results agree with those computed by expanding the generating function.

We also applied Proposition 1 to derive a formula, independent of all other methods, by which we calculated numbers (which we are not describing; see the “Six Little Squares” Web page [7]) that allowed us to find the 840 constituents by interpolation. These interpolated constituents fully agreed with the ones given above.

5 Magilatin squares of order 3

A magilatin square is like a semimagic square except that entries may be equal if they are in different rows and columns. The inside-out polytope is the same as with semimagic squares except that we omit those hyperplanes that prevent equality of entries in different rows and columns. Thus, in our count of reduced squares, we have to count the fractional lattice points in some of the faces of the polytope.

The reduced normal form of a magilatin square is the same as that of a semimagic square except that the restrictions are weaker. It might be thought that this would introduce ambiguity into the standard form because the minimum can occur in several cells, but it turns out that it does not.

Proposition 2. A reduced, normal 3 × 3 magilatin square has the form (13) with the restrictions

\begin{aligned}0 & < β, γ; \\ 0 & ≤ α; \\ 0 & ≤ δ ≤ β; \end{aligned}  (31)

and (15). Each reduced square with w in the upper right corner corresponds to exactly t − w − 1 different magilatin squares with entries in the range (0, t), for 0 < w < t. Each reduced square with magic sum s corresponds to one magilatin square with magic sum equal to t, if t ≡ s (mod 3), and none otherwise, for 0 < s < t.

Proof. The proof is similar to that for semimagic squares; we can arrange the square by permuting rows and columns and by reflection in the main diagonal so that x_{11} is the smallest entry, the first row and column are each increasing, and x_{21} ≥ x_{12}. We cannot say x_{21} > x_{12} because entries that do not share a row or column may be equal. Still, we obtain the form (13) with the bounds (31) and the same inequations (15) as in semimagic because all the latter depend on having no two equal values in the same line (row or column). ■

Each reduced, normal magilatin square gives rise to a family of true magilatin squares by adding a positive constant to each entry and by symmetries, which are generated by row and column permutations and reflection in the main diagonal. Call the set of symmetries 𝔊. As with semimagic, |𝔊| = 2(3!)^{2} = 72. Each normal, reduced square S gives rise to |𝔊/𝔊_{S}| = 72/|𝔊_{S}| squares via symmetries, where 𝔊_{S} is the stabilizer subgroup of S. If all entries are distinct, then the square is semimagic, 𝔊_{S} is trivial, and everything is as with semimagic squares. However, if α = 0 or δ = 0 or δ = β, the stabilizer is nontrivial. We consider each case in turn.

The case α = 0 < δ. Here δ < β because no line can repeat a value. To fix the square we cannot permute any rows or columns but we can reflect in the main diagonal, so |𝔊_{S}| = 2. Moreover, (15) reduces to

δ ≠ γ, \frac{β}{2} , \frac{β + γ}{2} .

We are in OBC, the x = 0 facet of Q, with the induced arrangement of three lines, ℐ^{x=0}. The number of reduced magilatin squares of this kind is 36r_{OBC}(t), where r_{OBC}(t) is the number of \frac{1}{t}-lattice points in the open facet and, equivalently, the number of symmetry types of reduced magilatin squares of this kind. We apply Equations (1)–(4) to the intersection poset ℒ(OBC^{∘}, ℐ^{x=0}), which is found in Figure 2. The Möbius function μ(OBC, u) equals (−1)^{codim u}.

The case δ = 0 < α. In this case a nontrivial member of 𝔊_{S} can only exchange the two zero positions. Such a symmetry that preserves the increase of the first row and column is unique (as one can easily see); thus |𝔊_{S}| = 2. Furthermore, (15) reduces to

α ≠ β.

We are in OAB, the z = 0 facet, with the induced arrangement ℐ^{z=0} of one line. The number of reduced magilatin squares of this kind is 36r_{OAB}(t), where r_{OAB}(t) is the number of \frac{1}{t}- lattice points in the open facet and, equivalently, the number of symmetry types of reduced squares. We apply Equations (1)–(4) to the intersection poset ℒ(OAB^{∘}, ℐ^{z=0}), shown in Figure 2. The Möbius function μ(OAB, u) equals (−1)^{codim u}.

The case δ = β. Here we must have α > 0. There are two zero positions in opposite corners. A symmetry that exchanges them and preserves increase in the first row and column is uniquely determined, so |𝔊_{S}| = 2. The inequations reduce to

β ≠ γ, α + γ.

We are in OAC, the facet where y = z, with the two-line induced arrangement ℐ^{y=z}. The number of reduced magilatin squares of this kind is 36r_{OAC}(t), where r_{OAC}(t) is the number of \frac{1}{t}-lattice points in the open facet, equally the number of reduced symmetry types. We apply Equations (1)–(4) to the intersection poset ℒ(OAC^{∘}, ℐ^{y=z}) in Figure 2. The Möbius function μ(OAC, u) equals (−1)^{codim u}.

The case α = 0 = δ. In these squares there are three zero positions and the whole square is a cyclic latin square. Any symmetry that fixes the zero positions also fixes the rest of the square. There are 3! symmetries that permute the zero positions, generated by row and column permutations. They all preserve the entire square. Therefore |𝔊_{S}| = 6. The inequations disappear. We are in the edge OB, which is the face where x = z = 0, with the empty arrangement, ℐ^{x=z=0} = ∅. The number of reduced magilatin squares of this kind is 12r_{OB}(t), where r_{OB}(t) is the number of \frac{1}{t}-lattice points in the open edge, also the number of reduced symmetry types. The intersection poset ℒ(OB^{∘}, ∅) consists of the one element OB, whose Möbius function μ(OB, OB) = 1.

To get the intersection posets we may examine Tables 1 and 3 to find the edges and vertices of (Q^{∘}, ℐ) in each closed facet. We also need to know which vertex is in which edge; this is easy. Although we do not need the fourth facet, ABC, we include it for the interest of its more complicated geometry.

Of course all the functions r_{s} and R_{ml} depend on whether we are counting cubically or affinely (thus, subscripted c or a); the two types will be treated separately. But the general conclusions hold that

R_{ml}(x) = 72r_{s}(x) + 36[r_{OAB}(x) + r_{OAC}(x) + r_{OBC}(x)] + 12r_{OB}(x),  (32)

and for the number of reduced symmetry types, r_{ml}(t) with generating function r_{ml}(x),

r_{ml}(x) = r_{s}(x) + r_{OAB}(x) + r_{OAC}(x) + r_{OBC}(x) + r_{OB}(x),  (33)

where r_{s}(x) is from semimagic and, by Equation (4) since |μ(X, Y)| = 1 for every lower interval in each facet poset (except facet ABC, which we do not use),

(−1)^{3}r_{OAB}(1/x) = E_{OAB}(x) + E_{OE}(x),
(−1)^{3}r_{OAC}(1/x) = E_{OAC}(x) + E_{AD}(x) + E_{DE′}(x),
(−1)^{3}r_{OBC}(1/x) = E_{OBC}(x) + E_{OE′′}(x) + E_{BD}(x) + E_{DE′′}(x) + E_{F} (x),
(−1)^{2}r_{OB}(1/x) = E_{OB}(x);

(34)

the sign and reciprocal on the left result from Equation (6).

There is also the generating function of the number of cubical or affine symmetry classes, l(t), whose generating function is l(x). This is obtained from r_{ml}(t) in the same way as L(t) is from R_{ml}(t), the exact way depending on whether the count is affine or cubic.

5.1 Magilatin squares: Cubical count (by upper bound)

The weak quasipolynomial is exactly as in the semimagic cubical problem.

PSfrag replacements

ℒ(OAB^{∘}, ℐ^{z=0}) ℒ(OAC^{∘}, ℐ^{y=z}) ℒ(OBC^{∘}, ℐ^{x=0})

ℒ(ABC^{∘}, ℐ^{x+y=1/2})

OAB OAC OBC ABC OE AD DE′ OE′′ BD DE′′ EE′′ BE′ AE′′ CE E′E′′

F

G G′ G′′

Figure 2: The four facet intersection posets of (Q^{∘}, ℐ).

5.1.1 Magilatin squares by upper bound

The number of 3 × 3 magilatin squares with strict upper bound t is L_{c}(t). We count them via R_{mlc}(w), the number of reduced magilatin squares with largest entry (which we know to be x_{13}) equal to w. The formula is

L_{c}(t) = ∑_{w=0}^{t−1}(t − 1 − w)R_{mlc}(w).  (35)

Equivalently, R_{mlc}(w) counts \frac{1}{w}-integral points in the interior and part of the boundary of the inside-out polytope (Q_{c}, ℐ_{c}) of Section 4.1, weighted variably by 72/|𝔊_{S}|.

Now we must calculate the closed Ehrhart generating function for each necessary face. This is done by LattE; here are the results. First, OAB_{c}:

\begin{matrix}(−1)^{3}r_{OAB_{c}}(1/x) = & E_{OAB_{c}}(x) + E_{OE_{c}}(x) \\ = & \frac{1}{(1 − x)^{2}(1 − x^{2})} + \frac{1}{(1 − x)(1 − x)^{3}} \\ = & \frac{x + 2}{(1 − x)(1 − x^{2})(1 − x^{3})} . \end{matrix}  (36)

Next is OAC_{c}:

(−1)^{3}r_{OAC_{c}}(1/x) = E_{OAC_{c}}(x) + E_{AD_{c}}(x) + E_{D_{c}E′_{c}}(x)
= \frac{1}{(1 − x)^{2}(1 − x^{2})} + \frac{1}{(1 − x^{2})^{2}} + \frac{1}{(1 − x^{2})(1 − x^{3})}
= \frac{x^{2} + 2x + 3}{(1 − x^{2})^{2}(1 − x^{3})} .

(37)

The last facet is OB_{c}C_{c}:

(−1)^{3}r_{OB_{c}C_{c}}(1/x) = E_{OB_{c}C_{c}}(x) + E_{OE′′_{c}}(x) + E_{B_{c}D_{c}}(x) + E_{D_{c}E′′_{c}}(x) + E_{F_{c}}(x)
= \frac{1}{(1 − x)^{3}} + \frac{1}{(1 − x)(1 − x^{2})} + \frac{1}{(1 − x)(1 − x^{2})}

+ \frac{2x^{2} + 1}{(1 − x^{2})^{2}} + \frac{1}{1 − x^{3}}

= \frac{−2x^{5} + 4x^{2} + 5x + 5}{(1 − x^{2})^{2}(1 − x^{3})} .

(38)

Finally, the edge OB_{c}:

(−1)^{2}r_{OB_{c}}(1/x) = E_{OB_{c}}(x) = \frac{1}{(1 − x)^{2}} .  (39)

Now R_{mlc}(x) results from (32), and then from (35) we see that

\begin{matrix}L_{c}(x) = \frac{x^{2}}{(1 − x)^{2}} R_{mlc}(x) = \\ \frac{12x^{4} \{ \begin{matrix}79x^{15} + 190x^{14} + 260x^{13} + 250x^{12} + 211x^{11} + 179x^{10} + 181x^{9} \\ + 198x^{8} + 210x^{7} + 181x^{6} + 125x^{5} + 61x^{4} + 22x^{3} + 8x^{2} + 4x + 1 \end{matrix}\}}{(1 − x^{4})(1 − x^{5})(1 − x^{3})^{2}(1 − x^{2})^{2}} . \end{matrix}  (40)

The constituents of L_{c} are extracted as described in Section 2.1, and here they are:

L_{c}(t) = \begin{cases}\frac{3}{10}t^{5} − \frac{51}{8}t^{4} + \frac{202}{3}t^{3} − \frac{3769}{10}t^{2} + c_{1}(t)t − c_{0}(t), & if t is even; \\ \frac{3}{10}t^{5} − \frac{51}{8}t^{4} + \frac{202}{3}t^{3} − \frac{7493}{20}t^{2} + c_{1}(t)t − c_{0}(t), & if t is odd; \end{cases}  (41)

where c_{1} varies with period 6, given by

c_{1}(t) = \begin{cases}994, & if t ≡ 0, 2 \\ 986, & if t ≡ 4 \\ \frac{1909}{2} , & if t ≡ 1 \\ \frac{1925}{2} , & if t ≡ 3, 5 \end{cases} (mod 6);

and c_{0}, given by Table 7, varies with period 60.

\begin{matrix}t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) & t & c_{0}(t) \\ 0 & 948 & 12 & 948 & 24 & \frac{4884}{5} & 36 & \frac{4452}{5} & 48 & \frac{4884}{5} \\ 1 & \frac{76933}{120} & 13 & \frac{87301}{120} & 25 & \frac{16769}{24} & 37 & \frac{16769}{24} & 49 & \frac{87301}{120} \\ 2 & \frac{2780}{3} & 14 & \frac{14332}{15} & 26 & \frac{13036}{15} & 38 & \frac{14332}{15} & 50 & \frac{2780}{3} \\ 3 & \frac{35607}{40} & 15 & \frac{6891}{8} & 27 & \frac{6891}{8} & 39 & \frac{35607}{40} & 51 & \frac{32151}{40} \\ 4 & \frac{13292}{15} & 16 & \frac{11996}{15} & 28 & \frac{13292}{15} & 40 & \frac{2572}{3} & 52 & \frac{2572}{3} \\ 5 & \frac{18433}{24} & 17 & \frac{18433}{24} & 29 & \frac{95621}{120} & 41 & \frac{85253}{120} & 53 & \frac{95621}{120} \\ 6 & \frac{4452}{5} & 18 & \frac{4884}{5} & 30 & 948 & 42 & 948 & 54 & \frac{4884}{5} \\ 7 & \frac{18497}{24} & 19 & \frac{95941}{120} & 31 & \frac{85573}{120} & 43 & \frac{95941}{120} & 55 & \frac{18497}{24} \\ 8 & \frac{14332}{15} & 20 & \frac{2780}{3} & 32 & \frac{2780}{3} & 44 & \frac{14332}{15} & 56 & \frac{13036}{15} \\ 9 & \frac{32727}{40} & 21 & \frac{29271}{40} & 33 & \frac{32727}{40} & 45 & \frac{6315}{8} & 57 & \frac{6315}{8} \\ 10 & \frac{2572}{3} & 22 & \frac{2572}{3} & 34 & \frac{13292}{15} & 46 & \frac{11996}{15} & 58 & \frac{13292}{15} \\ 11 & \frac{93893}{120} & 23 & \frac{104261}{120} & 35 & \frac{20161}{24} & 47 & \frac{20161}{24} & 59 & \frac{104261}{120} \end{matrix}

Table 7: Constant terms of the constituents of L_{c}(t), counting all magilatin squares by upper bound.

Thus the period of L_{c} turns out to be 60, just like that of S_{c} (not a surprise). However, again as with S_{c}, the even constant terms have half the period of the odd constant terms. That means L_{c}(2t) has period equal to half the denominator of the corresponding inside-out polytope (2P, ℋ). We have no explanation for this.

The principal constituent, that for t ≡ 0 (mod 60), is

\frac{3}{10}t^{5} − \frac{51}{8}t^{4} + \frac{202}{3}t^{3} − \frac{3769}{10}t^{2} + 994t − 948.

The constant term for magilatin squares does not have the simple interpretation as a number of linear orderings that it does for magic and semimagic squares, because the entries in the square need not all be different. (Still, there is an interpretation as a number of partial orderings of the nine cells; see Theorem 4.1 in our general magic and magilatin paper [6], and recall that an acyclic orientation of a graph can be represented by a partial ordering of the vertices.)

We confirmed the formulas by generating all magilatin squares and comparing the count with the coefficients of L_{c}(x) up to t = 91.

5.1.2 Symmetry types of magilatin squares, counted by upper bound

The number of symmetry types with strict bound t is l_{c}(t). We count them via r_{mlc}(w), given by Equation (33); then

l_{c}(t) = ∑_{w=0}^{t−1}(t − 1 − w) r_{mlc}(w).  (42)

Equivalently, r_{mlc}(w) counts the \frac{1}{w}-integral points in the interior and part of the boundary of the inside-out polytope (Q_{c}, ℐ_{c}) of Section 4.1.

We get r_{mlc}(x) from (33); then from (42) we see that

\begin{matrix}l_{c}(x) = \frac{x^{2}}{(1 − x)^{2}} r_{mlc}(x) = \\ \frac{x^{4} \left\{{9x^{15} + 20x^{14} + 23x^{13} + 16x^{12} + 10x^{11} + 13x^{10} + 27x^{9} + 43x^{8} \atop + 54x^{7} + 52x^{6} + 41x^{5} + 25x^{4} + 14x^{3} + 8x^{2} + 4x + 1}\right\}}{(1 − x^{2})^{2}(1 − x^{3})^{2}(1 − x^{4})(1 − x^{5})} . \end{matrix}  (43)

The constituents of l_{c} are:

l_{c}(t) = \begin{cases}\frac{1}{240}t^{5} − \frac{3}{64}t^{4} + \frac{97}{216}t^{3} − \frac{2029}{720}t^{2} + ĉ_{1}(t)t − ĉ_{0}(t), & if t is even; \\ \frac{1}{240}t^{5} − \frac{3}{64}t^{4} + \frac{97}{216}t^{3} − \frac{4013}{1440}t^{2} + ĉ_{1}(t)t − ĉ_{0}(t), & if t is odd; \end{cases}

where ĉ_{1} varies with period 6, given by

ĉ_{1}(t) = \begin{cases}\frac{17}{2} , & if t ≡ 0, 2 \\ \frac{151}{18} , & if t ≡ 4 \\ \frac{1163}{144} , & if t ≡ 1 \\ \frac{131}{16} , & if t ≡ 3, 5 \end{cases} (mod 6),

and ĉ_{0}, given by Table 8, varies with period 60.

\begin{matrix}t & ĉ_{0}(t) & t & ĉ_{0}(t) & t & ĉ_{0}(t) & t & ĉ_{0}(t) & t & ĉ_{0}(t) \\ 0 & 9 & 12 & 9 & 24 & \frac{47}{5} & 36 & \frac{41}{5} & 48 & \frac{47}{5} \\ 1 & \frac{49213}{8640} & 13 & \frac{59581}{8640} & 25 & \frac{11225}{1728} & 37 & \frac{11225}{1728} & 49 & \frac{59581}{8640} \\ 2 & \frac{235}{27} & 14 & \frac{1229}{135} & 26 & \frac{1067}{135} & 38 & \frac{1229}{135} & 50 & \frac{235}{27} \\ 3 & \frac{2823}{320} & 15 & \frac{539}{64} & 27 & \frac{539}{64} & 39 & \frac{2823}{320} & 51 & \frac{2439}{320} \\ 4 & \frac{1144}{135} & 16 & \frac{982}{135} & 28 & \frac{1144}{135} & 40 & \frac{218}{27} & 52 & \frac{218}{27} \\ 5 & \frac{12313}{1728} & 17 & \frac{12313}{1728} & 29 & \frac{65021}{8640} & 41 & \frac{54653}{8640} & 53 & \frac{65021}{8640} \\ 6 & \frac{41}{5} & 18 & \frac{47}{5} & 30 & 9 & 42 & 9 & 54 & \frac{47}{5} \\ 7 & \frac{12953}{1728} & 19 & \frac{68221}{8640} & 31 & \frac{57853}{8640} & 43 & \frac{68221}{8640} & 55 & \frac{12953}{1728} \\ 8 & \frac{1229}{135} & 20 & \frac{235}{27} & 32 & \frac{235}{27} & 44 & \frac{1229}{135} & 56 & \frac{1067}{135} \\ 9 & \frac{2503}{320} & 21 & \frac{2119}{320} & 33 & \frac{2503}{320} & 45 & \frac{475}{64} & 57 & \frac{475}{64} \\ 10 & \frac{218}{27} & 22 & \frac{218}{27} & 34 & \frac{1144}{135} & 46 & \frac{982}{135} & 58 & \frac{1144}{135} \\ 11 & \frac{63293}{8640} & 23 & \frac{73661}{8640} & 35 & \frac{14041}{1728} & 47 & \frac{14041}{1728} & 59 & \frac{73661}{8640} \end{matrix}

Table 8: Constant terms of the constituents of l_{c}(t), counting symmetry types of magilatin squares by upper bound.

Thus the period of l_{c} turns out to be 60. As with S_{c} and L_{c}, the period of the even constant terms is half that of the odd constant terms.

The principal constituent of l_{c} is

\frac{1}{240}t^{5} − \frac{3}{64}t^{4} + \frac{97}{216}t^{3} − \frac{2029}{720}t^{2} + \frac{17}{2}t − 9.

5.1.3 Some real numbers

For the first several nonzero values of the numbers of magilatin squares and of symmetry types, consult this table:

t 4 5 6 7  8  9  10 11  12  13  14  15
L_{c}(t) 12 48 120 384 1068 2472 4896 9072 15516 25608 40296 61608
l_{c}(t) 1 4 10 24 53 106 191 328 528 822 1230 1794
R_{mlc}(t) 12 24 36 192 420 720 1020 1752 2268 3648 4596 6624
r_{mlc}(t) 1 2 3 8 15 24 32 52  63  94  114 156

The third line contains the number of symmetry classes of 3 × 3 magilatin squares, counted by upper bound. The main numbers, L_{c}(t) and l_{c}(t), are sequences \underline{A173548} and \underline{A173729} in the OEIS [15]. The reduced numbers, R_{mlc}(t) and r_{mlc}(t), are sequences \underline{A174018} and \underline{A174019}. In contrast to the semimagic case, the number of squares is not a simple multiple of the number of symmetry types.

5.2 Magilatin squares: Affine count (by magic sum)

The last example is 3 × 3 magilatin squares, counted affinely. Let L_{a}(t) be the number of 3 × 3 magilatin squares with magic sum t > 0.

The weak quasipolynomial is the same as in affine semimagic.

5.2.1 Magilatin squares by magic sum

We compute L_{a}(t), the number of squares with magic sum t, via R_{mla}(s), the number of reduced squares with magic sum s. The formula is

L_{a}(t) = ∑_{\begin{matrix}0<s≤t−3 \\ s≡t (mod 3) \end{matrix}} R_{mla}(s)  if t > 0.  (44)

Equivalently, R_{mla}(s) counts \frac{1}{s}-integral points in the interior and part of the boundary of the inside-out polytope (Q_{a}, ℐ_{a}) of Section 4.2, each weighted by 72/|𝔊_{S}|.

Now we calculate (by LattE) the closed Ehrhart generating function for each necessary face. First, OAB_{a}:

\begin{matrix}(−1)^{3}r_{OAB_{a}}(1/x) = & E_{OAB_{a}}(x) + E_{OE_{a}}(x) \\ = & \frac{1}{(1 − x)(1 − x^{2})^{2}} + \frac{1}{(1 − x)(1 − x^{4})} \\ = & \frac{2}{(1 − x)(1 − x^{2})(1 − x^{4})} . \end{matrix}  (45)

Next is OAC_{a}:

(−1)^{3}r_{OAC_{a}}(1/x) = E_{OAC_{a}}(x) + E_{AD_{a}}(x) + E_{D_{a}E′_{a}}(x)
= \frac{1}{(1 − x)(1 − x^{2})^{2}} + \frac{1}{(1 − x^{2})(1 − x^{3})} + \frac{1}{(1 − x^{3})(1 − x^{4})}
= \frac{x^{3} + x^{2} + x + 3}{(1 − x^{2})(1 − x^{3})(1 − x^{4})} .

(46)

The last facet is OB_{a}C_{a}:

(−1)^{3}r_{OB_{a}C_{a}}(1/x) = E_{OB_{a}C_{a}}(x) + E_{OE′′_{a}}(x) + E_{B_{a}D_{a}}(x)

+ E_{D_{a}E′′_{a}}(x) + E_{F_{a}}(x)

= \frac{1}{(1 − x)(1 − x^{2})^{2}} + \frac{1}{(1 − x)(1 − x^{4})} + \frac{1}{(1 − x^{2})(1 − x^{3})}

+ \frac{1}{(1 − x^{3})(1 − x^{4})} + \frac{1}{1 − x^{5}}

= \frac{x^{5} + 4x^{4} + 6x^{3} + 7x^{2} + 7x + 5}{(1 − x^{5})(1 − x^{3})(1 − x^{4})(1 + x)} .

(47)

Finally, the edge OB_{a}:

(−1)^{2}r_{OB_{a}}(1/x) = E_{OB_{a}}(x) = \frac{1}{(1 − x)(1 − x^{2})} .  (48)

Now we get R_{mla}(x) from (32); then by (44) we deduce that

\begin{matrix}L_{a}(x) = \frac{x^{3}}{1 − x^{3}} R_{mla}(x) \\ = \frac{12x^{6} \{ \begin{matrix}1 + 3x + 7x^{2} + 15x^{3} + 33x^{4} + 65x^{5} + 128x^{6} + 208x^{7} + 316x^{8} \\ + 434x^{9} + 566x^{10} + 676x^{11} + 784x^{12} + 852x^{13} + 911x^{14} + 936x^{15} \\ + 967x^{16} + 967x^{17} + 1001x^{18} + 995x^{19} + 1000x^{20} + 955x^{21} \\ + 893x^{22} + 752x^{23} + 624x^{24} + 456x^{25} + 322x^{26} + 174x^{27} + 79x^{28} \end{matrix}\}}{(1 + x)(1 + x + x^{2})(1 + x^{2})(1 − x^{3})(1 − x^{5})(1 − x^{6})(1 − x^{7})(1 − x^{8})} . \end{matrix} (49)

The constituents of L_{a} are the following:

L_{a}(t) = \frac{1}{8}t^{4} − 3t^{3} + a_{2}(t)t^{2} − a_{1}(t)t + a_{0}(t) − 72S_{7}(t),  (50)

where a_{2} varies with period 6, given by

a_{2}(t) = \begin{cases}\frac{151}{4} , & if t ≡ 0 \\ \frac{135}{4} , & if t ≡ 2, 4 \\ \frac{63}{2} , & if t ≡ 1, 5 \\ \frac{71}{2} , & if t ≡ 3 \end{cases} (mod 6);

the linear coefficient varies with period 12, given by

a_{1}(t) = \begin{cases}\frac{1296}{5} , & if t ≡ 0 \\ \frac{1347}{10} , & if t ≡ 1, 5 \\ \frac{831}{5} , & if t ≡ 2, 10 \\ \frac{2097}{10} , & if t ≡ 3 \\ \frac{876}{5} , & if t ≡ 4, 8 \\ \frac{1251}{5} , & if t ≡ 6 \\ \frac{1257}{10} , & if t ≡ 7, 11 \\ \frac{2187}{10} , & if t ≡ 9 \end{cases} (mod 12);

the constant term a_{0}, given by Table 9, varies with period 120; and S_{7} is as in the affine semimagic count.

\begin{matrix}t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) & t & a_{0}(t) \\ 0 & 876 & 20 & 340 & 40 & 400 & 60 & 840 & 80 & 376 & 100 & 364 \\ 1 & \frac{4243}{40} & 21 & \frac{21843}{40} & 41 & \frac{3283}{40} & 61 & \frac{5683}{40} & 81 & \frac{20403}{40} & 101 & \frac{4723}{40} \\ 2 & \frac{1097}{5} & 22 & \frac{1037}{5} & 42 & \frac{3597}{5} & 62 & \frac{917}{5} & 82 & \frac{1217}{5} & 102 & \frac{3417}{5} \\ 3 & \frac{15219}{40} & 23 & \frac{−461}{40} & 43 & \frac{−941}{40} & 63 & \frac{16659}{40} & 83 & \frac{−1901}{40} & 103 & \frac{499}{40} \\ 4 & \frac{1604}{5} & 24 & \frac{4164}{5} & 44 & \frac{1484}{5} & 64 & \frac{1784}{5} & 84 & \frac{3984}{5} & 104 & \frac{1664}{5} \\ 5 & \frac{1463}{8} & 25 & \frac{1367}{8} & 45 & \frac{4887}{8} & 65 & \frac{1175}{8} & 85 & \frac{1655}{8} & 105 & \frac{4599}{8} \\ 6 & \frac{3201}{5} & 26 & \frac{881}{5} & 46 & \frac{821}{5} & 66 & \frac{3381}{5} & 86 & \frac{701}{5} & 106 & \frac{1001}{5} \\ 7 & \frac{3091}{40} & 27 & \frac{17811}{40} & 47 & \frac{2131}{40} & 67 & \frac{1651}{40} & 87 & \frac{19251}{40} & 107 & \frac{691}{40} \\ 8 & \frac{1448}{5} & 28 & \frac{1388}{5} & 48 & \frac{3948}{5} & 68 & \frac{1268}{5} & 88 & \frac{1568}{5} & 108 & \frac{3768}{5} \\ 9 & \frac{21267}{40} & 29 & \frac{5587}{40} & 49 & \frac{5107}{40} & 69 & \frac{22707}{40} & 89 & \frac{4147}{40} & 109 & \frac{6547}{40} \\ 10 & 265 & 30 & 705 & 50 & 241 & 70 & 229 & 90 & 741 & 110 & 205 \\ 11 & \frac{−1037}{40} & 31 & \frac{1363}{40} & 51 & \frac{16083}{40} & 71 & \frac{403}{40} & 91 & \frac{−77}{40} & 111 & \frac{17523}{40} \\ 12 & \frac{4092}{5} & 32 & \frac{1772}{5} & 52 & \frac{1712}{5} & 72 & \frac{4272}{5} & 92 & \frac{1592}{5} & 112 & \frac{1892}{5} \\ 13 & \frac{4819}{40} & 33 & \frac{19539}{40} & 53 & \frac{3859}{40} & 73 & \frac{3379}{40} & 93 & \frac{20979}{40} & 113 & \frac{2419}{40} \\ 14 & \frac{809}{5} & 34 & \frac{1109}{5} & 54 & \frac{3309}{5} & 74 & \frac{989}{5} & 94 & \frac{929}{5} & 114 & \frac{3489}{5} \\ 15 & \frac{4023}{8} & 35 & \frac{311}{8} & 55 & \frac{791}{8} & 75 & \frac{3735}{8} & 95 & \frac{599}{8} & 115 & \frac{503}{8} \\ 16 & \frac{1676}{5} & 36 & \frac{3876}{5} & 56 & \frac{1556}{5} & 76 & \frac{1496}{5} & 96 & \frac{4056}{5} & 116 & \frac{1376}{5} \\ 17 & \frac{5011}{40} & 37 & \frac{7411}{40} & 57 & \frac{22131}{40} & 77 & \frac{6451}{40} & 97 & \frac{5971}{40} & 117 & \frac{23571}{40} \\ 18 & \frac{3273}{5} & 38 & \frac{593}{5} & 58 & \frac{893}{5} & 78 & \frac{3093}{5} & 98 & \frac{773}{5} & 118 & \frac{713}{5} \\ 19 & \frac{787}{40} & 39 & \frac{18387}{40} & 59 & \frac{−173}{40} & 79 & \frac{2227}{40} & 99 & \frac{16947}{40} & 119 & \frac{1267}{40} \end{matrix}

Table 9: Constant terms of the truncated constituents of L_{a}(t), counting magilatin squares by magic sum.

The period of L_{a} turns out to be 840—the period of the constant terms, due to the combination of a_{0}(t) and the constant term of S_{7}(t). This is equal to the denominator.

The principal constituent of L_{a}, that is, for t ≡ 0 (mod 840), is

\frac{1}{8}t^{4} − 3t^{3} + \frac{151}{4}t^{2} − \frac{9192}{35}t + 948

(incorporating the −72S_{7} term). (As with the cubical magilatin count, there is an interpretation of the constant term 948 in terms of partial orderings; see Theorem 4.7 in our general paper [6].)

We verified the results by comparing an actual count of magilatin squares with magic sum t ≤ 100 to the coefficients in L_{a}.

5.2.2 Symmetry types of magilatin squares, counted by magic sum

We compute l_{a}(t), the number of symmetry types of squares with magic sum t, via r_{mla}(s), the number of reduced symmetry types with magic sum s. The formula is

l_{a}(t) = ∑_{\begin{matrix}0<s≤t−3 \\ s≡t (mod 3) \end{matrix}} r_{mla}(s)  if t > 0.  (51)

Equivalently, r_{mla}(s) counts \frac{1}{s}-integral points in the interior and part of the boundary of the inside-out polytope(Q_{a}, ℐ_{a}) of Section 4.2.

From (32) we get r_{mla}(x) and then from (51) we see that

\begin{matrix}l_{a}(x) = \frac{x^{3}}{1 − x^{3}} r_{mla}(x) \\ = \frac{x^{6} \{ \begin{matrix}1 + 3x + 7x^{2} + 13x^{3} + 23x^{4} + 37x^{5} + 60x^{6} + 86x^{7} + 118x^{8} + 149x^{9} \\ + 180x^{10} + 199x^{11} + 212x^{12} + 208x^{13} + 196x^{14} + 171x^{15} \\ + 145x^{16} + 115x^{17} + 96x^{18} + 79x^{19} + 72x^{20} + 67x^{21} \\ + 66x^{22} + 59x^{23} + 54x^{24} + 43x^{25} + 33x^{26} + 19x^{27} + 9x^{28} \end{matrix}\}}{(1 + x)(1 + x + x^{2})(1 + x^{2})(1 − x^{3})(1 − x^{5})(1 − x^{6})(1 − x^{7})(1 − x^{8})} . \end{matrix} (52)

The constituents of l_{a} are:

l_{a}(t) = \frac{1}{576}t^{4} − \frac{1}{48}t^{3} + â_{2}(t)t^{2} − â_{1}(t)t + â_{0}(t) − S_{7}(t),

where the quadratic term varies with period 6, given by

â_{2}(t) = \begin{cases}\frac{25}{96} , & if t ≡ 0 \\ \frac{25}{144} , & if t ≡ 1, 5 \\ \frac{59}{288} , & if t ≡ 2, 4 \\ \frac{11}{48} , & if t ≡ 3 \end{cases} (mod 6);

the linear term varies with period 12, given by

â_{1}(t) = \begin{cases}\frac{31}{15} , & if t ≡ 0 \\ \frac{103}{120} , & if t ≡ 1, 5 \\ \frac{133}{120} , & if t ≡ 2, 10 \\ \frac{47}{30} , & if t ≡ 3 \\ \frac{37}{30} , & if t ≡ 4, 8 \\ \frac{233}{120} , & if t ≡ 6 \\ \frac{11}{15} , & if t ≡ 7, 11 \\ \frac{203}{120} , & if t ≡ 9 \end{cases} (mod 12);

and â_{0}, given by Table 10, varies with period 120.

\begin{matrix}t & â_{0}(t) & t & â_{0}(t) & t & â_{0}(t) & t & â_{0}(t) & t & â_{0}(t) & t & â_{0}(t) \\ 0 & 8 & 20 & \frac{47}{18} & 40 & \frac{31}{9} & 60 & \frac{15}{2} & 80 & \frac{28}{9} & 100 & \frac{53}{18} \\ 1 & \frac{2027}{2880} & 21 & \frac{1523}{320} & 41 & \frac{1067}{2880} & 61 & \frac{3467}{2880} & 81 & \frac{1363}{320} & 101 & \frac{2507}{2880} \\ 2 & \frac{553}{360} & 22 & \frac{493}{360} & 42 & \frac{257}{40} & 62 & \frac{373}{360} & 82 & \frac{673}{360} & 102 & \frac{237}{40} \\ 3 & \frac{979}{320} & 23 & \frac{−949}{2880} & 43 & \frac{−1429}{2880} & 63 & \frac{1139}{320} & 83 & \frac{−2389}{2880} & 103 & \frac{11}{2880} \\ 4 & \frac{229}{90} & 24 & \frac{38}{5} & 44 & \frac{199}{90} & 64 & \frac{137}{45} & 84 & \frac{71}{10} & 104 & \frac{122}{45} \\ 5 & \frac{847}{576} & 25 & \frac{751}{576} & 45 & \frac{343}{64} & 65 & \frac{559}{576} & 85 & \frac{1039}{576} & 105 & \frac{311}{64} \\ 6 & \frac{221}{40} & 26 & \frac{409}{360} & 46 & \frac{349}{360} & 66 & \frac{241}{40} & 86 & \frac{229}{360} & 106 & \frac{529}{360} \\ 7 & \frac{1739}{2880} & 27 & \frac{1171}{320} & 47 & \frac{779}{2880} & 67 & \frac{299}{2880} & 87 & \frac{1331}{320} & 107 & \frac{−661}{2880} \\ 8 & \frac{104}{45} & 28 & \frac{193}{90} & 48 & \frac{36}{5} & 68 & \frac{163}{90} & 88 & \frac{119}{45} & 108 & \frac{67}{10} \\ 9 & \frac{1427}{320} & 29 & \frac{3083}{2880} & 49 & \frac{2603}{2880} & 69 & \frac{1587}{320} & 89 & \frac{1643}{2880} & 109 & \frac{4043}{2880} \\ 10 & \frac{149}{72} & 30 & \frac{49}{8} & 50 & \frac{125}{72} & 70 & \frac{113}{72} & 90 & \frac{53}{8} & 110 & \frac{89}{72} \\ 11 & \frac{−1813}{2880} & 31 & \frac{587}{2880} & 51 & \frac{1043}{320} & 71 & \frac{−373}{2880} & 91 & \frac{−853}{2880} & 111 & \frac{1203}{320} \\ 12 & \frac{73}{10} & 32 & \frac{131}{45} & 52 & \frac{247}{90} & 72 & \frac{39}{5} & 92 & \frac{217}{90} & 112 & \frac{146}{45} \\ 13 & \frac{2891}{2880} & 33 & \frac{1299}{320} & 53 & \frac{1931}{2880} & 73 & \frac{1451}{2880} & 93 & \frac{1459}{320} & 113 & \frac{491}{2880} \\ 14 & \frac{301}{360} & 34 & \frac{601}{360} & 54 & \frac{229}{40} & 74 & \frac{481}{360} & 94 & \frac{421}{360} & 114 & \frac{249}{40} \\ 15 & \frac{279}{64} & 35 & \frac{−17}{576} & 55 & \frac{463}{576} & 75 & \frac{247}{64} & 95 & \frac{271}{576} & 115 & \frac{175}{576} \\ 16 & \frac{128}{45} & 36 & \frac{69}{10} & 56 & \frac{113}{45} & 76 & \frac{211}{90} & 96 & \frac{37}{5} & 116 & \frac{181}{90} \\ 17 & \frac{2219}{2880} & 37 & \frac{4619}{2880} & 57 & \frac{1491}{320} & 77 & \frac{3659}{2880} & 97 & \frac{3179}{2880} & 117 & \frac{1651}{320} \\ 18 & \frac{233}{40} & 38 & \frac{157}{360} & 58 & \frac{457}{360} & 78 & \frac{213}{40} & 98 & \frac{337}{360} & 118 & \frac{277}{360} \\ 19 & \frac{−277}{2880} & 39 & \frac{1267}{320} & 59 & \frac{−1237}{2880} & 79 & \frac{1163}{2880} & 99 & \frac{1107}{320} & 119 & \frac{203}{2880} \end{matrix}

Table 10: Constant terms of the truncated constituents of l_{a}(t), the number of symmetry types of magilatin squares with given magic sum.

The period of l_{a} is 840. The principal constituent, that for t ≡ 0 (mod 840), is

\frac{1}{576}t^{4} − \frac{1}{48}t^{3} + \frac{25}{96}t^{2} − \frac{74}{35}t + 9.

(This incorporates the −S_{7} term.)

5.2.3 Some numbers

The first several nonzero values are given in the table. The third line gives the number of symmetry classes of squares.

t 6 7 8 9 10 11 12 13 14 15 16 17 18 19
L_{a}(t) 12 12 24 72 156 240 552 600 1020 1548 2004 2568 4008 4644
l_{a}(t) 1 1 2 4 7 10 20 22 35 50 63 78 116 131
R_{mla}(t) 12 12 24 60 144 216 480 444 780 996 1404 1548 2640 3696
r_{mla}(t) 1 1 2 3 6 8 16 15 25 30 41 43 66 68

The sequences L_{a}(t) and l_{a}(t) are \underline{A173549} and \underline{A173730} in the OEIS [15]. The numbers R_{mla}(t) of reduced magilatin squares and r_{mla}(t) of normalized, reduced squares with largest value t are sequences \underline{A174020} and \underline{A174021}.

6 Observations and conjectures

A remarkable fact is that the period of every one of our strong Ehrhart quasipolynomials equals the denominator, when it could be much smaller.

For some small values of t we calculated M_{a}(t) = E_{P^{∘},ℋ}^{∘} (t) by hand, which is feasible because the problem is 2-dimensional. The process of counting lattice points in a diagram drew our attention to some remarkable phenomena that apply to the semimagic and magilatin problems as well. Let δ := dim s; let c_{k} be the coefficient in the quasipolynomial E_{P^{∘},ℋ}^{∘} (t) and let c_{k}^{w} be that in the Ehrhart quasipolynomial of P, and let p_{k}, p_{k}^{w} be their periods. We observe that the variation in c_{δ−1} is exactly the same as that in c_{δ−1}^{w}, i.e.,

c_{δ−1}(t) − c_{δ−1}(t − 1) = c_{δ−1}^{w}(t) − c_{δ−1}^{w}(t − 1),

but that is not so for most lower coefficients, especially c_{0}. The reason is that adding each new excluded hyperplane results in a constant deduction in degree δ − 1 (as we discussed at Equation (4.9) in our first article [5]) but a more irregular one in lower terms. We observe that p_{k} increases—that is, there is longer-term variation in c_{k}—as k decreases in every case. Thus we propose some daring conjectures.

Conjecture 3. In an inside-out counting problem, let δ := dim P.

(a) p_{k}^{w} | p_{k} for 0 ≤ k ≤ δ. (We know that p_{δ−1} = p_{δ−1}^{w} because c_{δ−1} and c_{δ−1}^{w} have the same variation.)

(b) If p_{j} = p_{j}^{w} for all j ≥ k, then the variation in c_{k} is the same as that in c_{k}^{w}.

(c) The period ratios increase by a multiplicative factor as k decreases:

\frac{p_{k}}{p_{k}^{w}} | \frac{p_{k−1}}{p_{k−1}^{w}} for 0 ≤ k ≤ δ.

We do not suggest p_{k} | p_{k+1} because that is false in general in ordinary Ehrhart theory, according to McAllister and Woods [13]. However, it might be true for the kinds of insideout polytopes that arise in cubical and affine counting.

7 Acknowledgments

The authors are grateful to the anonymous referee for several helpful suggestions.

References

[1] Maya Ahmed, Jesús De Loera, and Raymond Hemmecke, Polyhedral cones of magic cubes and squares, in Boris Aronov et al., eds., Discrete and Computational Geometry: The Goodman–Pollack Festschrift, Algorithms Combin., Vol. 25, Springer-Verlag, 2003, pp. 25–41. MR 2004m:05016. Zbl 1077.52506.

[2] Matthias Beck, Moshe Cohen, Jessica Cuomo, and Paul Gribelyuk, The number of “magic” squares and hypercubes, Amer. Math. Monthly 110 (2003), 707–717. MR 2004k:05009. Zbl 1043.05501.

[3] Matthias Beck and Andrew van Herick, Enumeration of 4 × 4 magic squares, Math. Comp., to appear.

[4] Matthias Beck and Dennis Pixton, The Ehrhart polynomial of the Birkhoff polytope, Discrete Comput. Geom. 30 (2003), 623–637. MR 2004g:52015. Zbl 1065.52007.

[5] Matthias Beck and Thomas Zaslavsky, Inside-out polytopes, Adv. Math. 205 (1) (2006), 134–162. MR 2007e:52017. Zbl 1107.52009.

[6] Matthias Beck and Thomas Zaslavsky, An enumerative geometry for magic and magilatin labellings, Ann. Combin. 10 (2006), 395–413. MR 2007m:05010. Zbl 1116.05071.

[7] Matthias Beck and Thomas Zaslavsky, “Six Little Squares and How their Numbers Grow” Web Site: With a detailed version (as of January 25, 2007) of this paper, Maple worksheets, and supporting documentation.

URL http://www.math.binghamton.edu/zaslav/Tmath/SLSfiles/

[8] Garrett Birkhoff, Three observations on linear algebra (in Spanish), Rev. Univ. Nac. Tucumán, Ser. A 5 (1946), 147–151. MR 8, 561a. Zbl 60, 79f (e: 060.07906).

[9] Jesús A. De Loera, David Haws, Raymond Hemmecke, Peter Huggins, Jeremy Tauzer, and Ruriko Yoshida, A User’s Guide for LattE v1.1, Univ. of California at Davis, 2003,

URL http://www.math.ecdavis.edu/~ latte/

[10] Eugène Ehrhart, Sur un problème de géometrie diophantienne linéaire. I: Polyèdres et réseaux. II: Systèmes diophantiens linéaires, J. reine angew. Math. 226 (1967), 1–29; 227 (1967), 25–49. Correction, ibid. 231 (1968), 220. MR 35 #4184, 36 #105. Zbl 155.37503, 164.05304.

[11] Eugène Ehrhart, Sur les carrés magiques, C. R. Acad. Sci. Paris Sér. A–B 277 (1973), A651–A654. MR 48 #10859. Zbl 267.05014.

[12] Percy A. MacMahon, Combinatory Analysis, Vols. I and II, Cambridge University Press, 1915–1916. Reprinted in one volume by Chelsea, 1960. MR 25 #5003. Zbl 101, 251 (e: 101.25102).

[13] Tyrrell B. McAllister and Kevin M. Woods, The minimum period of the Ehrhart quasipolynomial of a rational polytope, J. Combin. Theory Ser. A 109 (2005), 345–352. MR 2005i:52018. Zbl 1063.52006.

[14] Jochi Shigeru, The dawn of wasan (Japanese mathematics), in Helaine Selin and Ubiratan D’Ambrosio, eds., Mathematics Across Cultures: The History of Non-Western Mathematics, Science Across Cultures: The History of Non-Western Science, Vol. 2, Kluwer, 2000, pp. 423–454. MR (book) 1805670 (2002a:01001). Zbl 981.01005.

[15] N. J. A. Sloane, The On-Line Encyclopedia of Integer Sequences,

    http://www.research.att.com/~njas/sequences/

[16] Richard P. Stanley, Linear homogeneous Diophantine equations and magic labelings of graphs, Duke Math. J. 40 (1973), 607–632. MR 47 #6519. Zbl 269.05109.

[17] Richard P. Stanley, Enumerative Combinatorics, Vol. I, Wadsworth & Brooks/Cole, 1986. MR 87j:05003. Zbl 608.05001. Corrected reprint, Cambridge Stud. Adv. Math., Vol. 49, Cambridge University Press, 1997. MR 98a:05001. Zbl 889.05001, 945.05006.

[18] Guoce Xin, Constructing all magic squares of order three, Discrete Math. 308 (2008), 3393–3398. MR 2423421 (2009e:05039). Zbl 1145.05012.

2010 Mathematics Subject Classification: Primary 05B15; Secondary 05A15, 52B20, 52C35. Keywords: magic square, semimagic square, magic graph, latin square, magilatin square, lattice-point counting, rational convex polytope, arrangement of hyperplanes.

(Concerned with sequences \underline{A108235}, \underline{A108236}, \underline{A108576}, \underline{A108577}, \underline{A108578}, \underline{A108579}, \underline{A173546}, \underline{A173547}, \underline{A173548}, \underline{A173549}, \underline{A173723}, \underline{A173724}, \underline{A173725}, \underline{A173726}, \underline{A173727}, \underline{A173728}, \underline{A173729}, \underline{A173730}, \underline{A174018}, \underline{A174019}, \underline{A174020}, \underline{A174021}, \underline{A174256}, and \underline{A174257}.)

Received March 9 2010; revised version received June 1 2010. Published in Journal of Integer Sequences, June 2 2010. Revised, June 8 2010.

Return to Journal of Integer Sequences home page.