The Ehrhart polynomial of the Birkhoff polytope ^{1}

Matthias Beck and Dennis Pixton

Abstract: The n^{th} Birkhoff polytope is the set of all doubly stochastic n × n matrices, that is, those matrices with nonnegative real coefficients in which every row and column sums to one. A wide open problem concerns the volumes of these polytopes, which have been known for n ≤ 8. We present a new, complex-analytic way to compute the Ehrhart polynomial of the Birkhoff polytope, that is, the function counting the integer points in the dilated polytope. The leading term of the Ehrhart polynomial is—up to a trivial factor—the volume of the polytope, which is one reason why we are interested in this counting function. We implemented our methods in form of a computer program, which yielded the Ehrhart polynomial (and hence the volume) of the ninth Birkhoff polytope.

1 Introduction

One of the most intriguing objects of combinatorial geometry is the n^{th} Birkhoff polytope

B_{n} = \{( \begin{matrix}x_{11} & ⋯ & x_{1n} \\ ⋮ & & ⋮ \\ x_{n1} & … & x_{nn} \end{matrix} ) ∈ ℝ^{n^{2}} : x_{jk} ≥ 0, \begin{matrix}∑_{j} x_{jk} = 1 for all 1 ≤ k ≤ n \\ ∑_{k} x_{jk} = 1 for all 1 ≤ j ≤ n \end{matrix} \} ,

often described as the set of all n × n doubly stochastic matrices. B_{n} is a convex polytope with integer vertices. It possesses fascinating combinatorial properties [4, 5, 6, 8, 22] and relates to many mathematical areas [10, 14]. A long-standing open problem is the determination of the relative volume of B_{n}, which has been known only up to n = 8 [7]. In this paper, we propose a new method of calculating this volume and use it to compute vol B_{9}.

One of the recent attempts to compute vol B_{n} relies on the theory of counting functions for the integer points in polytopes. Ehrhart proved [11] that for a polytope P ⊂ ℝ^{d} with integral vertices, the number

L_{P}(t) := # (tP ∩ ℤ^{d})

is a polynomial in the positive integer variable t. He showed various other properties for this counting function (in fact, in the more general setting of P having rational vertices), of which we mention three here:

• The degree of L_{P} is the dimension of P.

• The leading term of L_{P} is the relative volume of P, normalized with respect to the sublattice of ℤ^{d} on the affine subspace spanned by P.

• Since L_{P} is a polynomial, we can evaluate it at nonpositive integers. These evaluations yield

L_{P}(0) = χ(P) ,  (1)
L_{P}(−t) = (−1)^{dim P}L_{P^{∘}}(t) .  (2)

1Version of 2/22/2002

Here χ(P) denotes the Euler characteristic, P^{∘} the relative interior of P. The reciprocity law (2) was in its full generality proved by Macdonald [15].

The application of this theory to the Birkhoff polytope B_{n} incorporates the nice interpretation of the number of integral points in tB_{n} as the number of semi-magic squares, namely, square matrices whose nonnegative integral coefficients sum up to the same integer t along each row and column.

We will denote the Ehrhart polynomial of B_{n} by

H_{n}(t) := L_{B_{n}}(t) .

It is not hard to see that dim B_{n} = (n − 1)^{2}, hence H_{n} is a polynomial in t of degree (n − 1)^{2}. The first two of these polynomials are trivial:

H_{1}(t) = 1 H_{2}(t) = t + 1 ,

the first nontrivial case was computed by MacMahon [16] as

H_{3}(t) = 3\binom{t + 3}{4} + \binom{t + 2}{2} .  (3)

The structural properties of H_{n} were first studied in [12, 18, 19]. It is a nice exercise to deduce from (2) that

H_{n}(−n − t) = (−1)^{n−1} H_{n}(t)  (4)

and

H_{n}(−1) = H_{n}(−2) = ⋯ = H_{n}(−n + 1) = 0 .

This allows the following strategy of computing H_{n}, and therefore, the volume of B_{n}: compute the first \binom{n−1}{2} values of H_{n}, use the above symmetry and trivial values of H_{n}, and calculate the polynomial H_{n} by interpolation. In fact, as far as we are aware of, the volume of B_{8} was computed using essentially this method, combined with some nice computational tricks [7].

We propose a new, completely different approach of computing H_{n} (and hence vol B_{n}). It is based on an analytic method by the first author of computing the Ehrhart polynomial of a polytope [2]. We will introduce the application of this method to the Birkhoff polytope in the following section.

Some recent refreshing approaches—of a more algebraic-geometric/topological flavor—to the problem of computing vol B_{n} and H_{n} can be found in [1, 9, 21].

2 An integral counting integers

A convex polytope P ⊂ ℝ^{d} is an intersection of halfspaces. This allows the compact description

P = {x ∈ ℝ^{d} : Ax ≤ b} ,

for some (m × d)-matrix A and m-dimensional vector b. Here the inequality is understood componentwise. In fact, we may convert these inequalities into equalities by introducing ‘slack variables.’ If P has rational vertices (those polytopes are called rational), we can choose A and b in such a way that all their entries are integers. In summary, we may assume that a convex rational polytope P is given by

P = {x ∈ ℝ_{≥0}^{d} : Ax = b} ,  (5)

where A ∈ M_{m×d}(ℤ) and b ∈ ℤ^{m}. (If we are interested in counting the integer points in P, we may assume that P is in the nonnegative orthant, i.e., the points in P have nonnegative coordinates, as translation by an integer vector does not change the lattice-point count.) The following straightforward theorem can be found in [2]. We use the standard multivariate notation v^{w} := v_{1}^{w_{1}} ⋯ v_{n}^{w_{n}}.

Theorem 1 [2, Theorem 8] Suppose the convex rational polytope P is given by (5), and denote the columns of A by c_{1}, …, c_{d}. Then

L_{P}(t) = \frac{1}{(2πi)^{m}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{m}|=ϵ_{m}} \frac{z_{1}^{−tb_{1}−1} ⋯ z_{m}^{−tb_{m}−1}}{(1 − z^{c_{1}}) ⋯ (1 − z^{c_{d}})} dz .

Here 0 < ϵ_{1}, …, ϵ_{m} < 1 are different numbers such that we can expand all the \frac{1}{1−z^{c_{k}}} into power series about 0.

We can view the Birkhoff polytope B_{n} as given in the form of (5), where

A = ( \begin{matrix}1 & ⋯ & 1 \\ & & & 1 & ⋯ & 1 \\ . & . & . \\ & & & & & & & 1 & ⋯ & 1 \\ 1 & & & 1 & & & & 1 \\ & ... & & & ... & & ⋯ & & ... \\ & & 1 & & & 1 & & & & 1 \end{matrix} )

is an (2n × n^{2})-matrix and b = (1, …, 1) ∈ ℤ^{2n}. Hence Theorem 1 gives for this special case

H_{n}(t) = \frac{1}{(2πi)^{2n}} ∫ ⋯ ∫ \frac{(z_{1} ⋯ z_{2n})^{−t−1}}{(1 − z_{1}z_{n+1})(1 − z_{1}z_{n+2}) ⋯ (1 − z_{n}z_{2n})} dz .

Here it is understood that each integral is over a circle with radius < 1 centered at 0; all appearing radii should be different. We can separate, say, the last n variables and obtain

H_{n}(t) = \frac{1}{(2πi)^{n}} ∫ ⋯ ∫ (z_{1} ⋯ z_{n})^{−t−1} (\frac{1}{2πi} ∫ \frac{z^{−t−1}}{(1 − z_{1}z) ⋯ (1 − z_{n}z)} dz)^{n} dz_{n} ⋯ dz_{1} .

The innermost integral is easy to compute: It is equal to the residue at 0 of

\frac{1}{z^{t+1}(1 − z_{1}z) ⋯ (1 − z_{n}z)}

and, by the residue theorem, equal to the negative of the sum of the residues at \frac{1}{z_{1}} , …, \frac{1}{z_{n}}. (Note that here we use the fact that t > 0.) The residues at these simple poles are easily computed, and we obtain the starting point for our computations.

Theorem 2 For any pairwise different 0 < ϵ_{1}, …, ϵ_{n} < 1,

H_{n}(t) = \frac{1}{(2πi)^{n}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}}(z_{1} ⋯ z_{n})^{−t−1} (∑_{k=1}^{n} \frac{z_{k}^{t+n−1}}{∏_{j≠k}(z_{k} − z_{j})})^{n} dz_{n} ⋯ dz_{1} .

Remark. It is not hard to conclude from the form of the integrand that H_{n} is indeed a polynomial in t. In fact, one can make this property more apparent by noticing that the expression in parenthesis is actually a polynomial; namely

H_{n}(t) = \frac{1}{(2πi)^{n}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}}(z_{1} ⋯ z_{n})^{−t−1} ( ∑_{m_{1}+···+m_{n}=t} z_{1}^{m_{1}} ⋯ z_{n}^{m_{n}})^{n} dz_{n} ⋯ dz_{1} ,

where the sum is over all ordered partitions of t. This formula can also be proved “more directly” combinatorially.^{2}

3 Small n do not require a computer

We will now illustrate the computation of H_{n} (and hence vol B_{n}) by means of Theorem 2 for n = 3 and 4. These calculations “by hand” give an idea what computational tricks one might use in tackling larger n with the aid of a computer.

By the theorem,

H_{3}(t) = \frac{1}{(2πi)^{3}} ∫ (z_{1}z_{2}z_{3})^{−t−1} (\frac{z_{1}^{t+2}}{(z_{1} − z_{2})(z_{1} − z_{3})} + \frac{z_{2}^{t+2}}{(z_{2} − z_{1})(z_{2} − z_{3})} + \frac{z_{3}^{t+2}}{(z_{3} − z_{1})(z_{3} − z_{2})})^{3} dz.

We have to order the radii of the integration paths for each variable; we choose 0 < ϵ_{3} < ϵ_{2} < ϵ_{1} < 1. We heavily use this fact after multiply out the cubic: integrating, for example, the term

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

with respect to z_{3} gives 0, as this function is analytic at the z_{3}-origin and |z_{1}| , |z_{2}| > ϵ_{3}. After exploiting this observation for all the terms stemming from the cubic, the only integrals surviving are

\frac{1}{(2πi)^{3}} ∫ \frac{z_{1}^{2t+5}z_{2}^{−t−1}z_{3}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{3}} dz

and

−\frac{3}{(2πi)^{3}} ∫ \frac{z_{1}^{t+3}z_{2}z_{3}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{2}(z_{2} − z_{3})} dz .

2The authors thank Sinai Robins and Frank Sottile for their help in the proof of this equivalence and its combinatorial interpretation.

The first integral factors and yields, again by residue calculus,

\frac{1}{(2πi)^{3}} ∫ \frac{z_{1}^{2t+5}z_{2}^{−t−1}z_{3}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{3}} dz = \frac{1}{(2πi)^{3}} ∫ z_{1}^{2t+5} (∫ \frac{z^{−t−1}}{(z_{1} − z)^{3}(z_{1} − z)^{3}} dz)^{2} dz_{1} = \binom{t + 2}{2}^{2} .

For the second integral, it is most efficient to integrate with respect to z_{2} first:

−\frac{3}{(2πi)^{3}} ∫ \frac{z_{1}^{t+3}z_{2}z_{3}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{2}(z_{2} − z_{3})} dz = −\frac{3}{(2πi)^{2}} ∫ \frac{z_{1}^{t+3}z_{3}^{−t}}{(z_{1} − z_{3})^{5}} dz_{3} dz_{1} = −3\binom{t + 3}{4} .

Adding up the last two lines gives finally

H_{3}(t) = \binom{t + 2}{2}^{2} − 3\binom{t + 3}{4} = \frac{1}{8}t^{4} + \frac{3}{4}t^{3} + \frac{15}{8}t^{2} + \frac{9}{4}t + 1 ,

which is equal to (3). To obtain the volume, the leading term of H_{3} has to be multiplied by the relative volume of the fundamental domain of the sublattice of ℤ^{9} spanned by B_{3}, which is 9; hence

vol B_{3} = \frac{9}{8} .

In general, it is not hard to convince oneself that the relative volume of the fundamental domain of the sublattice of ℤ^{n^{2}} spanned by B_{n} is n^{n−1}.

The number of integrals we have to evaluate to compute H_{4} is only slightly higher. By Theorem 2,

H_{4}(t) = \frac{1}{(2πi)^{4}} ∫_{|z_{1}|=ϵ_{1}} ∫_{|z_{2}|=ϵ_{2}} ∫_{|z_{3}|=ϵ_{3}} ∫_{|z_{4}|=ϵ_{4}}(z_{1}z_{2}z_{3}z_{4})^{−t−1} (∑_{k=1}^{4} \frac{z_{k}^{t+3}}{∏_{j≠k}(z_{k} − z_{j})})^{4} dz .

Again we have a choice of ordering the radii; we use 0 < ϵ_{4} < ϵ_{3} < ϵ_{2} < ϵ_{1} < 1. After multiplying out the quartic, we have to calculate five integrals; their evaluation—again straightforward by means of the residue theorem—is as follows. As before, we can ‘save’ computation effort by choosing a particular order with which we integrate.

\frac{1}{(2πi)^{4}} ∫ \frac{z_{1}^{3t+11}z_{2}^{−t−1}z_{3}^{−t−1}z_{4}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{3}(z_{1} − z_{4})^{3}} dz = \frac{1}{(2πi)^{4}} ∫ z_{1}^{3t+11} (∫ \frac{z^{−t−1}}{(z_{1} − z)^{3}} dz)^{3} dz_{1}

= \binom{t + 3}{3}^{3},

− \frac{4}{(2πi)^{4}} ∫ \frac{z_{1}^{2t+8}z_{2}^{2}z_{3}^{−t−1}z_{4}^{−t−1}}{(z_{1} − z_{2})^{4}(z_{1} − z_{3})^{3}(z_{1} − z_{4})^{3}(z_{2} − z_{3})(z_{2} − z_{4})} dz

= −\frac{4}{(2πi)^{4}} ∫ \frac{z_{1}^{2t+8}z_{2}^{2}}{(z_{1} − z_{2})^{4}} (∫ \frac{z^{−t−1}}{(z_{1} − z)^{3}(z_{2} − z)} dz)^{2} dz_{1} dz_{2}

= \frac{4}{(2πi)^{2}} ∫ \frac{z_{1}^{2t+8}z_{2}^{−t+1}}{(z_{1} − z_{2})^{7}} (2\binom{t + 2}{2}\frac{z_{1}^{−t−3}}{z_{1} − z_{2}} + 2(t + 1)\frac{z_{1}^{−t−2}}{(z_{1} − z_{2})^{2}} + 2\frac{z_{1}^{−t−1}}{(z_{1} − z_{2})^{3}} − \frac{z_{2}^{−t−1}}{(z_{1} − z_{2})^{3}} ) dz

= 8\binom{t + 2}{2}\binom{t + 5}{7} + 8(t + 1)\binom{t + 6}{8} + 8\binom{t + 7}{9} − 4\binom{2t + 8}{9} ,

\frac{4}{(2πi)^{4}} ∫ \frac{z_{1}^{2t+8}z_{2}^{−t−1}z_{3}^{2}z_{4}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{4}(z_{1} − z_{4})^{3}(z_{2} − z_{3})(z_{3} − z_{4})} dz

= \frac{4}{(2πi)^{3}} ∫ \frac{z_{1}^{2t+8}z_{2}^{−t+1}z_{4}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{4})^{7}(z_{2} − z_{4})} dz

= \frac{4}{(2πi)^{2}} ∫ \frac{z_{1}^{2t+8}z_{4}^{−t−1}}{(z_{1} − z_{4})^{7}} (\binom{t + 2}{2}\frac{z_{1}^{−t−3}}{z_{1} − z_{4}} + (t + 1)\frac{z_{1}^{−t−2}}{(z_{1} − z_{4})^{2}} + \frac{z_{1}^{−t−1}}{(z_{1} − z_{4})^{3}} ) dz

= 4 (\binom{t + 2}{2}\binom{t + 5}{7} + (t + 1)\binom{t + 6}{8} + \binom{t + 7}{9}) ,

\frac{6}{(2πi)^{4}} ∫ \frac{z_{1}^{t+5}z_{2}^{t+5}z_{3}^{−t−1}z_{4}^{−t−1}}{(z_{1} − z_{2})^{4}(z_{1} − z_{3})^{2}(z_{1} − z_{4})^{2}(z_{2} − z_{3})^{2}(z_{2} − z_{4})^{2}} dz

= \frac{6}{(2πi)^{4}} ∫ \frac{z_{1}^{t+5}z_{2}^{t+5}}{(z_{1} − z_{2})^{4}} (∫ \frac{z^{−t−1}}{(z_{1} − z)^{2}(z_{2} − z)^{2}} dz)^{2} dz_{1} dz_{2}

= \frac{6}{(2πi)^{2}} ∫ \frac{z_{1}^{t+5}z_{2}^{t+5}}{(z_{1} − z_{2})^{4}} ((t + 1)^{2} \frac{z_{2}^{−2t−4}}{(z_{1} − z_{2})^{4}} − 4(t + 1)\frac{z_{2}^{−2t−3}}{(z_{1} − z_{2})^{5}} + 4\frac{z_{2}^{−2t−2}}{(z_{1} − z_{2})^{6}} ) dz

= 6(t + 1)^{2}\binom{t + 5}{7} − 24(t + 1)\binom{t + 5}{8} + 24\binom{t + 5}{9} ,

− \frac{12}{(2πi)^{4}} ∫ \frac{z_{1}^{t+5}z_{2}^{2}z_{3}^{2}z_{4}^{−t−1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{3})^{3}(z_{1} − z_{4})^{2}(z_{2} − z_{3})^{2}(z_{2} − z_{4})(z_{3} − z_{4})} dz

= −\frac{12}{(2πi)^{3}} ∫ \frac{z_{1}^{t+5}z_{2}^{2}z_{4}^{−t+1}}{(z_{1} − z_{2})^{3}(z_{1} − z_{4})^{5}(z_{2} − z_{4})^{3}} dz

= −\frac{12}{(2πi)^{2}} ∫ \frac{z_{1}^{t+5}z_{4}^{−t+1}}{(z_{1} − z_{4})^{5}} (\frac{1}{(z_{1} − z_{4})^{3}} + 6\frac{z_{4}}{(z_{1} − z_{4})^{4}} + 6\frac{z_{4}^{2}}{(z_{1} − z_{4})^{5}} ) dz

= −12\binom{t + 5}{7} − 72\binom{t + 5}{8} − 72\binom{t + 5}{9} .

Adding them all up gives

H_{4}(t) = \binom{t + 3}{3}^{3} + 6 (2t^{2} + 5t + 1) \binom{t + 5}{7} − 24 (t + 4) \binom{t + 5}{8} + 12(t + 1)\binom{t + 6}{8}

− 4\binom{2t + 8}{9} − 48\binom{t + 5}{9} + 12\binom{t + 7}{9}

= \frac{11}{11340} t^{9} + \frac{11}{630} t^{8} + \frac{19}{135} t^{7} + \frac{2}{3} t^{6} + \frac{1109}{540} t^{5} + \frac{43}{10} t^{4} + \frac{35117}{5670} t^{3} + \frac{379}{63} t^{2} + \frac{65}{18} t + 1

and hence

vol B_{4} = 4^{3} · \frac{11}{11340} = \frac{176}{2835} .

4 Larger n do

As we have seen in the examples, after multiplying out the integrand of Theorem 2, many of the terms do not contribute to the integral. The following proposition will provide us with a general statement to that extend.

For a rational function f in n variables z_{j} we use the notation d_{r}(f) for the degree of f in the variables z_{1}, …, z_{r}.

Proposition 3 Suppose p_{1}, …, p_{n} are integers, q_{jk} are nonnegative integers (1 ≤ j < k ≤ n), 1 > ϵ_{1} > ⋯ > ϵ_{n} > 0, and

f(z_{1}, …, z_{n}) = \frac{∏_{1≤j≤n} z_{j}^{p_{j}}}{∏_{1≤j<k≤n}(z_{j} − z_{k})^{q_{jk}}} .

If 1 ≤ r ≤ n and d_{r}(f) < −r then

∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}} f(z) dz = 0 .

Proof. This is clear if r = 1 since f, considered as a function of z_{1}, is analytic outside the circle |z_{1}| = ϵ_{1} and has zero residue at infinity since its degree is less than −1.

We continue by induction, so suppose r > 1 and the proposition is true for smaller r. Suppose that d_{r}(f) < −r; we may assume d_{r−1} ≥ −(r − 1). Note that d_{r}(f) ≥ d_{r−1}(f) + d where d is the degree of f in the single variable z_{r} (the discrepancy is the sum of the exponents q_{jr} for j < r). Hence d ≤ d_{r}(f) − d_{r−1}(f) < −1. We consider f as a function of z_{r} and apply the residue theorem to the region outside the circle |z_{r}| = ϵ_{r}. As above f has zero residue at infinity, so we only need to consider the residues at the poles z_{j} for j < r. Evaluating these residues converts the integral of f into a (possibly huge) linear combination of integrals of functions of the same form as f, but in the n − 1 variables z_{1}, …, z_{r−1}, z_{r+1}, …, z_{n}. If g is any one of these functions then we easily calculate d_{r−1}(g) = d_{r}(f) + 1 < −(r − 1), and, by the induction hypothesis, the integral of g is zero. □

From this proposition we obtain the starting point for our ‘algorithm.’

Corollary 4 For 1 > ϵ_{1} > ⋯ > ϵ_{n} > 0 and t ≥ 0, H_{n}(t) =

\frac{1}{(2πi)^{n}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}}(z_{1} ⋯ z_{n})^{−t−1} ∑_{m_{1}+···+m_{n}=n}^{∗}  \binom{n}{m_{1}, …, m_{n}} ∏_{k=1}^{n} (\frac{z_{k}^{t+n−1}}{∏_{j≠k}(z_{k} − z_{j})})^{m_{k}} dz ,

where ∑^{∗} denotes that we only sum over those n-tuples of non-negative integers satisfying m_{1} + ⋯ + m_{n} = n and m_{1} + ⋯ + m_{r} > r if 1 ≤ r < n.

Remark. The condition on m_{1}, …, m_{n} can be visualized through lattice paths from (0, 0) to (n, n) using the steps (1, m_{1}), (1, m_{2}), …, (1, m_{n}). The condition means that these paths stay strictly above the diagonal (except at the start and end).^{3}

Proof. By Theorem 2,

H_{n}(t) = \frac{1}{(2πi)^{n}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}}(z_{1} ⋯ z_{n})^{−t−1} (∑_{k=1}^{n} \frac{z_{k}^{t+n−1}}{∏_{j≠k}(z_{k} − z_{j})})^{n} dz

= \frac{1}{(2πi)^{n}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{n}|=ϵ_{n}}(z_{1} ⋯ z_{n})^{−t−1} ∑_{m_{1}+···+m_{n}=n} \binom{n}{m_{1}, …, m_{n}} ∏_{k=1}^{n} (\frac{z_{k}^{t+n−1}}{∏_{j≠k}(z_{k} − z_{j})})^{m_{k}} dz .

We select a partition m_{1} + ⋯ + m_{n} = n and rewrite the corresponding integrand in the language of Proposition 3:

p_{j} = (n − 1)m_{j} + t(m_{j} − 1) − 1 , q_{jk} = m_{j} + m_{k} .

Now suppose 1 ≤ r < n. The degree of the denominator in z_{1}, …, z_{r} is

∑_{j=1}^{r} ∑_{k=j+1}^{n}(m_{j} + m_{k}) = (n − 1) ∑_{j=1}^{r} m_{j} + r ∑_{j=r+1}^{n} m_{j} = (n − 1) ∑_{j=1}^{r} m_{j} + nr − r ∑_{j=1}^{r} m_{j}

= (n − 1) ∑_{j=1}^{r} m_{j} + nr − r ∑_{j=1}^{r}(m_{j} − 1) − r^{2} .

We subtract this from ∑_{j=1}^{r} p_{j} to get

d_{r}(f) = (t + r) ∑_{j=1}^{r}(m_{j} − 1) − r(n − r) − r .

Here t is non-negative and r(n − r) is positive, so Proposition 3 implies that the integral is zero unless ∑_{j=1}^{r}(m_{j} − 1) > 0. □

Theoretically, Corollary 4 tells us what we have to do to compute H_{n}. For practical purposes, its statement is almost worthless. For once, the number of terms in the sum equals the n^{th} Catalan number

\frac{1}{n + 1}\binom{2n}{n} = \frac{(2n)!}{(n + 1)!n!} ,

which grows exponentially with n. Another slippery point is the evaluation of each integral. As we have seen in the examples, and as can be easily seen for the general case, we can compute each integral step by step one variable at a time. However, this means at each step we ‘convert’ a rational function into a sum of rational functions (of one variable less) by means of the residue theorem. Again, this means that the number of (single-variable) integrals we have to compute grows immensely as n increases. In fact, if we just ‘feed’ the statement of Corollary 4 into a computer and tell it to integrate each summand, say, starting with z_{1}, then z_{2}, and so on, the computation time explodes once one tries n = 7 or 8. We feel that computationally this is as involved as calculating a sufficient number of values of H_{n} and then interpolating this polynomial. However, complex analysis allows us some shortcuts which turn out to speed up the computation by a huge factor and which make Corollary 4 valuable, even from a computational perspective. These ‘tricks’ all showed up already in the examples and include

3The authors thank Lou Billera for pointing out this lattice-path interpretation.

1. realizing when a function is analytic at the z_{k}-origin,

2. choosing the most efficient order of variables to integrate, and

3. factoring the integral if some of the variables appear in a symmetric fashion.

Finally, if we are only interested in the volume of B_{n}, we may also be

4. dispensing a particular integral if it does not contribute to the leading term of H_{n}.

It is worth noting that every of these computational ‘speed-ups’ decreases the total computation time substantially. By applying them to Corollary 4, we implemented a C-program for the specific functions we have to integrate to compute H_{n}.^{4} We were able to verify all previous values (n ≤ 8) and to compute vol B_{9} and H_{9}. The 65 coefficients of this polynomial can be found at www.math.binghamton.edu/dennis/Birkhoff. The following table gives the volumes together with the respective computing time (on a 1GHz PC running under Linux). Note that computing the full polynomial H_{n} takes longer, as we cannot make use of shortcut #4. In fact, the computation of H_{9} took more than 300 days. (We divided the process onto several computers.)

n  vol B_{n}  time
1  1  < .01 sec
2  2  < .01 sec
3  \frac{9}{8}  < .01 sec
4  \frac{176}{2835}  < .01 sec
5  \frac{23590375}{167382319104}  < .01 sec
6  \frac{9700106723}{1319281996032·10^{6}}  .18 sec
7  \frac{77436678274508929033}{137302963682235238399868928·10^{8}}  15 sec
8  \frac{5562533838576105333259507434329}{12589036260095477950081480942693339803308928·10^{10}}  54 min
9 \frac{559498129702796022246895686372766052475496691}{92692623409952636498965146712806984296051951329202419606108477153345536·10^{14}} 317 hr

4The program can be found at www.math.binghamton.edu/dennis/Birkhoff.

5 An outlook towards transportation polytopes

The Birkhoff polytopes are special cases of transportation polytopes, which are defined below. The study of this class of polytopes, which are naturally at least as fascinating as the Birkhoff polytopes, was motivated by problems in linear programming; for combinatorial properties see, for example, [10, 13]. The goal of this section is to show how our methods can be applied in this more general setting.

Fix positive real numbers a_{1}, …, a_{m}, b_{1}, …, b_{n} such that a_{1} + ⋯ + a_{m} = b_{1} + ⋯ + b_{n}. Let a = (a_{1}, …, a_{m}), b = (b_{1}, …, b_{n}).

T_{a,b} = \{( \begin{matrix}x_{11} & ⋯ & x_{1n} \\ ⋮ & & ⋮ \\ x_{m1} & … & x_{mn} \end{matrix} ) ∈ ℝ^{mn} : x_{jk} ≥ 0, \begin{matrix}∑_{k} x_{jk} = a_{j} for all 1 ≤ j ≤ m, \\ ∑_{j} x_{jk} = b_{k} for all 1 ≤ k ≤ n \end{matrix} \}

is the set of solutions to the transportation problem with parameters a, b. It is a convex (m − 1)(n − 1)-dimensional polytope in ℝ^{mn}—hence we refer to the T_{a,b} as transportation polytopes. Note that B_{n} = T_{1,1} where 1 = (1, …, 1) ∈ ℤ^{n}. We will be exclusively interested in the case when a_{1}, …, a_{m}, b_{1}, …, b_{n} are integers. One reason for this is that the vertices of T_{a,b} are then integral. Let

T (a, b) = T (a_{1}, …, a_{m}, b_{1}, …, b_{n}) = # (T_{a,b} ∩ ℤ^{mn})

denote the number of integer points in the transportation polytope with integral parameters a, b. Geometrically, each of these parameters determines the position of a hyperplane bounding the polytope T_{a,b}. It is well known (see, for example, [20]) that T (a, b) is a piecewise-defined polynomial in a_{1}, …, a_{m}, b_{1}, …, b_{n}, namely, it is a (different) polynomial for each region in the (m + n)- dimensional space of the parameters a_{1}, …, a_{m}, b_{1}, …, b_{n} in which the combinatorial type (that is, the face structure) of T_{a,b} does not change.

In fact, one can derive results for the counting function T (a, b) which are ‘higher-dimensional’ extensions of (4), by application of the main theorem in [3], which in turn is a generalization of the Ehrhart-Macdonald reciprocity theorem (2). To state this theorem, we need to define the following integer-point counting function. Suppose P is a rational polytope given in the form (5); denote by P(t) = {x ∈ ℝ^{d} : Ax = t} a polytope which we obtain from P = P(b) by translating (some of) its bounding hyperplanes. (Classical Ehrhart dilation is the special case t = tb.) Let

L_{P}(t) = # (P(t) ∩ ℤ^{d}) .

[3, Theorem 4] states that L_{P} and L_{P^{∘}} are piecewise defined (multivariable) quasipolynomials satisfying

L_{P}(−t) = (−1)^{dim P}L_{P^{∘}}(t) .

As easy as we deduced before (4) from (2), this reciprocity theorem yields a symmetry result for the transporation counting function. Denote the d-dimensional vector all of whose entries are one by 1_{d}. Then

T (−a − n1_{m}, −b − m1_{n}) = (−1)^{mn−1} T (a, b) .

We finally turn to the problem of writing T (a, b) in form of an integral. As with B_{n}, we can view T_{a,b} as given in the form of (5) and apply Theorem 1 to obtain

T (a, b) = \frac{1}{(2πi)^{m+n}} ∫ \frac{z_{1}^{−a_{1}−1} ⋯ z_{m}^{−a_{m}−1}w_{1}^{−b_{1}−1} ⋯ w_{n}^{−b_{n}−1}}{∏_{{1≤j≤m \atop 1≤k≤n}} (1 − z_{j}w_{k})} dz dw .

Again the integral with respect to each (one-dimensional) variable is over a circle with radius < 1 centered at 0, and all appearing radii are different. As in the Birkhoff case, we can separate, say, the w-variables and obtain

T (a, b) = \frac{1}{(2πi)^{m+n}} ∫ z_{1}^{−a_{1}−1} ⋯ z_{m}^{−a_{m}−1} ∏_{k=1}^{n} ∫ \frac{w_{k}^{−b_{k}−1}}{∏_{j=1}^{m} (1 − z_{j}w_{k})} dw_{k} dz .

And as with the Birkhoff polytope, the innermost integrals are easy to compute by means of the residue theorem. This yields the following statement which generalizes Theorem 2.

Theorem 5 For any pairwise different 0 < ϵ_{1}, …, ϵ_{m} < 1,

T (a, b) = \frac{1}{(2πi)^{m}} ∫_{|z_{1}|=ϵ_{1}} ⋯ ∫_{|z_{m}|=ϵ_{m}} z_{1}^{−a_{1}−1} ⋯ z_{m}^{−a_{m}−1} ∏_{k=1}^{n} ∑_{j=1}^{m} \frac{z_{j}^{b_{k}+m−1}}{∏_{l≠j}(z_{j} − z_{l})} dz .

Remark. It is not hard to conclude from the form of the integrand that T (a, b) is indeed a piecewisedefined polynomial in a_{1}, …, a_{m}, b_{1}, …, b_{n}.

Theoretically we could now use this theorem to produce formulas for T (a, b) just as we did for H_{n}(t). There is one major difference: T (a, b) is only a piecewise-defined polynomial. In fact, we can see this from the form of the integral in Theorem 5: whether we will get a nonzero contribution at a certain step in the computation depends heavily on the relationship between a_{1}, …, a_{m}, b_{1}, …, b_{n}.

The fact that the counting function T (a, b) is of a somewhat more delicate nature naturally has computational consequences. We feel that providing any general results on this function would go beyond the scope of this article and will hopefully be subject of a future project. On the other hand, we adjusted our algorithm to examples of values of T (a, b) for three (fixed) pairs (a, b), which have been previously computed by Mount [17] and DeLoera and Sturmfels [9]. The first example is

T ((3046, 5173, 6116, 10928), (182, 778, 3635, 9558, 11110)) =

23196436596128897574829611531938753

Its computation took 20/10 minutes with Mount’s/DeLoera–Sturmfels’s algorithm, respectively [9, 17]. We computed this number in 0.2 seconds, based on Theorem 5. A similar phenomenon happens with

T ((338106, 574203, 678876, 1213008), (20202, 142746, 410755, 1007773, 1222717)) =

316052820930116909459822049052149787748004963058022997262397 and

T ((30201, 59791, 70017, 41731, 58270), (81016, 68993, 47000, 43001, 20000)) = 24640538268151981086397018033422264050757251133401758112509495633028

which took 35 minutes/10 days with the DeLoera–Sturmfels algorithm [9], 0.3/2.9 seconds with ours.

References

[1] Welleda Baldoni-Silva and Michèle Vergne, Residues formulae for volumes and Ehrhart polynomials of convex polytopes, Preprint (arXiv:math.CO/0103097), 2001.

[2] Matthias Beck, Counting lattice points by means of the residue theorem, Ramanujan J. 4 (2000), no. 3, 299–310. MR 2001k:11191

[3] Matthias Beck, Multidimensional Ehrhart reciprocity, J. Combin. Theory Ser. A 97 (2002), no. 1, 187–194.

[4] Louis J. Billera and A. Sarangarajan, The combinatorics of permutation polytopes, Formal power series and algebraic combinatorics (New Brunswick, NJ, 1994), Amer. Math. Soc., Providence, RI, 1996, pp. 1–23. MR 96m:52014

[5] Richard A. Brualdi and Peter M. Gibson, Convex polyhedra of doubly stochastic matrices. I. Applications of the permanent function, J. Combinatorial Theory Ser. A 22 (1977), no. 2, 194–230. MR 55 #10486

[6] , Convex polyhedra of doubly stochastic matrices. III. Affine and combinatorial properties of U_{n}, J. Combinatorial Theory Ser. A 22 (1977), no. 3, 338–351. MR 55 #10488

[7] Clara S. Chan and David P. Robbins, On the volume of the polytope of doubly stochastic matrices, Experiment. Math. 8 (1999), no. 3, 291–300. MR 2000k:15048

[8] Clara S. Chan, David P. Robbins, and David S. Yuen, On the volume of a certain polytope, Experiment. Math. 9 (2000), no. 1, 91–99. MR 2001d:52022

[9] Jesús DeLoera and Bernd Sturmfels, Algebraic unimodular counting, Preprint

(arXiv:math.CO/0104286), 2001.

[10] Persi Diaconis and Anil Gangolli, Rectangular arrays with fixed margins, Discrete probability and algorithms (Minneapolis, MN, 1993), Springer, New York, 1995, pp. 15–41. MR 97e:05013

[11] Eugène Ehrhart, Sur un problème de géométrie diophantienne linéaire. II. Systèmes diophantiens linéaires, J. Reine Angew. Math. 227 (1967), 25–49. MR 36 #105

[12] , Sur les carrés magiques, C. R. Acad. Sci. Paris Sér. A-B 277 (1973), A651–A654. MR 48 #10859 [13] Victor Klee and Christoph Witzgall, Facets and vertices of transportation polytopes, Mathe-

matics of the Decision Sciences, Part I (Seminar, Stanford, Calif., 1967), Amer. Math. Soc., Providence, R.I., 1968, pp. 257–282. MR 38 #4134

[14] Donald E. Knuth, Permutations, matrices, and generalized Young tableaux, Pacific J. Math. 34 (1970), 709–727. MR 42 #7535

[15] I. G. Macdonald, Polynomials associated with finite cell-complexes, J. London Math. Soc. (2) 4 (1971), 181–192. MR 45 #7594

[16] Percy A. MacMahon, Combinatory analysis, Chelsea Publishing Co., New York, 1960. MR 25 #5003

[17] John Mount, Fast unimodular counting, Combin. Probab. Comput. 9 (2000), no. 3, 277–285. MR 2002a:05016

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

[19] M. L. Stein and P. R. Stein, Enumeration of stochastic matrices with integer elements, Report No. LA-4434, Los Alamos Scientific Laboratory, University of California, Los Alamos, New Mexico (1970).

[20] Bernd Sturmfels, On vector partition functions, J. Combin. Theory Ser. A 72 (1995), no. 2, 302–309. MR 97b:52014

[21] , Equations defining toric varieties, Algebraic geometry—Santa Cruz 1995, Amer. Math. Soc., Providence, RI, 1997, pp. 437–449. MR 99b:14058

[22] Doron Zeilberger, Proof of a conjecture of Chan, Robbins, and Yuen, Electron. Trans. Numer. Anal. 9 (1999), 147–148 (electronic), Orthogonal polynomials: numerical and symbolic algorithms (Leganés, 1998). MR 2001b:05015

Department of Mathematical Sciences State University of New York Binghamton, NY 13902-6000

matthias@math.binghamton.edu
dennis@math.binghamton.edu