1 Introduction↩︎

Let \(\lambda/\mu\) be a skew shape and \(\mathbf{w}= (w_1, \dotsc, w_k)\) a weight vector with \(|\mathbf{w}| = |\lambda/\mu|\). The Gelfand–Tsetlin polytopeGelfand–Tsetlin polytope \(\mathrm{GT}(\lambda/\mu, \mathbf{w})\) is the set of all real-valued GT patterns of shape \(\lambda/\mu\) and weight \(\mathbf{w}\). Its integer points are in bijection with semistandard Young tableaux of shape \(\lambda/\mu\) and weight \(\mathbf{w}\), so \[\label{eq:kostka} K_{\lambda/\mu, \mathbf{w}} = \#(\mathrm{GT}(\lambda/\mu, \mathbf{w}) \cap \mathbb{Z}^d),\tag{1}\] where \(d = \dim \mathrm{GT}(\lambda/\mu, \mathbf{w})\).

By a theorem of E. Rassart [1], the function \(n \mapsto K_{n\lambda/n\mu, n\mathbf{w}}\) is a polynomial in \(n\): the Ehrhart polynomialEhrhart polynomial \(L_P(n)\) of the polytope \(P = \mathrm{GT}(\lambda/\mu, \mathbf{w})\). Computing this polynomial is important for understanding the structure of Kostka coefficients and the geometry of GT polytopes.

The polynomiality of stretched Kostka and Littlewood–Richardson coefficients, and the structure of the resulting Ehrhart polynomials, have been studied by R. King, C. Tollu, and F. Toumazet [2], [3]. Their later hive-model paper [3] computes several explicit stretched LR polynomials as illustrative examples: the degree-8 polynomial for \((\lambda,\mu,\nu)=((4,3,3,2,1),(4,3,2,2,1),(7,4,4,4,3,2,1))\); a factorization example with \((\lambda,\mu,\nu)=((9,7,6,2,0),(13,5,3,1,0),(14,12,11,5,4))\); and degree-bound examples including \((\lambda,\mu,\nu)=((7,6,5,4,0,0),(7,7,7,4,0,0),(12,8,8,7,6,4,2))\) and \((\lambda,\mu,\nu)=((9,7,3,0,0),(9,9,3,2,0),(10,9,9,8,6))\). These examples support conjectures on factorization, degree bounds, and linear factors; they are not a systematic benchmark computation of the family of GT/Kostka polytopes considered here. Nzeutchap, Toumazet, and Butelle [4] later developed a distributed method for stretched Kostka and Littlewood–Richardson coefficients based on the hive model. Their approach computes values at positive dilation factors across a local network and interpolates the resulting stretched polynomial; it also uses hive inequalities to obtain an improved upper bound, denoted maxDeg, on the degree before the computation begins. Their reported LR examples are arbitrary selected triples, with maxDeg between 4 and 14; in several cases the real degree is strictly smaller, including examples of real degree 0. In contrast, for the GT/Kostka family considered here we compute the dimension of the relevant GT polytope directly, so the exact degree of the stretched Kostka polynomial is known before interpolation. The reciprocity step then reduces the number and cost of the needed evaluations, since many negative evaluations are strict pattern counts and often vanish.

Kostka coefficients arise as special cases of Littlewood–Richardson (LR) coefficients. By the identity \(K_{\lambda/\nu,\mu} = \langle s_{\lambda/\nu}, h_\mu \rangle = \langle s_\lambda, s_\nu h_\mu \rangle\), we construct a skew shape \(\rho/\sigma\) that is a disjoint union (no shared rows or columns) of \(\nu\) and single rows of lengths \(\mu_1, \mu_2, \dotsc\), giving \(s_{\rho/\sigma} = s_\nu h_\mu\) and hence \(K_{\lambda/\nu,\mu} = c^{\rho}_{\lambda, \sigma}\). Thus software for LR coefficients such as A. Buch’s lrcalc [5] can be used for individual Kostka evaluations. However, the dilated partitions \(n\lambda\), \(n\mathbf{w}\) grow with \(n\), making LR computation expensive at large dilations. Even lrcalc (highly optimized C code) takes over 8 minutes for the 11 evaluations needed for \(\mathrm{GT}((3,3,3), (1^9))\) at dimension 10, compared to 1 ms for our dynamic programming (DP) approach with reciprocity (Table 1). Our horizontal-strip DP exploits the specific structure of Kostka coefficients directly, without the overhead of the LR reduction.

The standard approach evaluates \(L_P(n)\) at \(n = 1, 2, \dotsc, d+1\) using a DP algorithm for Kostka coefficients, then recovers the polynomial by Lagrange interpolation. Each evaluation requires a DP that can be expensive for large dilations.

Our approach uses Ehrhart–Macdonald reciprocityEhrhart–Macdonald reciprocity to access the negative side of the polynomial for free or cheaply. The same idea also applies in other settings where one can count both lattice points and relative interior lattice points efficiently; in particular, we will later discuss order polytopes and an application to permutation posets.

2 Background↩︎

2.1 Ehrhart–Macdonald reciprocity↩︎

For a rational polytope \(P \subset \mathbb{R}^d\) of dimension \(d\) (one whose vertices have rational coordinates), the function \(L_P(n) = \#(nP \cap \mathbb{Z}^d)\) is a quasi-polynomial in \(n\), and a polynomial when \(P\) is a lattice polytope. GT polytopes are rational but in general not lattice polytopes, so one would a priori expect \(L_P\) to be only a quasi-polynomial. However, by Rassart’s theorem [1], \(L_P\) is a genuine polynomial—a case of so-called period collapse. Ehrhart–Macdonald reciprocity [6], [7] gives: \[\label{eq:reciprocity} (-1)^d L_P(-n) = \relax\ifmmode\textcolor{blue}{L^{*}_P(n)}\else \textcolor{blue}{\emph{L^{*}_P(n)}}\fi \mathrel{\vcenter{:}}= \#(\mathop{\mathrm{int}}(nP) \cap \mathbb{Z}^d),\tag{2}\] where \(\mathop{\mathrm{int}}(nP)\)\(\mathop{\mathrm{int}}(nP)\) denotes the relative interior of \(nP\).

2.2 Gelfand–Tsetlin polytopes↩︎

Recall that an SSYT of skew shape \(\lambda/\mu\) and content \(\mathbf{w}= (w_1, \dotsc, w_k)\) corresponds bijectively to a chain of partitions \[\mu = \alpha^{(0)} \subseteq \alpha^{(1)} \subseteq \dotsb \subseteq \alpha^{(k)} = \lambda,\] where each \(\alpha^{(i)} / \alpha^{(i-1)}\) is a horizontal strip of size \(w_i\). The Gelfand–Tsetlin polytopeGelfand–Tsetlin polytope \(\mathrm{GT}(\lambda/\mu, \mathbf{w}) \subset \mathbb{R}^m\) is the convex polytope obtained by relaxing the integrality constraint: the entries of the intermediate partitions \(\alpha^{(1)}, \dotsc, \alpha^{(k-1)}\) are allowed to be real, subject to the interlacing inequalities \[\alpha^{(i-1)}_j \leq \alpha^{(i)}_j \leq \alpha^{(i-1)}_{j-1}\] and the row-sum constraints \(|\alpha^{(i)}| - |\alpha^{(i-1)}| = w_i\). The lattice points of \(\mathrm{GT}(\lambda/\mu, \mathbf{w})\) are precisely the SSYT, so \(K_{\lambda/\mu, \mathbf{w}} = \#(\mathrm{GT}(\lambda/\mu, \mathbf{w}) \cap \mathbb{Z}^m)\).

2.3 Strict GT patterns↩︎

The interior lattice points of \(n \cdot \mathrm{GT}(\lambda/\mu, \mathbf{w})\) are precisely the strict GT patternsstrict GT patterns: those where every interlacing inequality that is not forced as an equality in the affine hull of the polytope is strict. Equivalently, one imposes strictness only on the non-globally-tight inequalities. Thus \[\label{eq:strict-kostka} L^{*}_P(n) = K^{\mathrm{strict}}_{n\lambda/n\mu, n\mathbf{w}},\tag{3}\] where \(K^{\mathrm{strict}}\) counts strict semistandard tableaux (equivalently, GT patterns with strictness on all non-forced internal inequalities). This is computed by the same DP as the ordinary Kostka coefficient, but with strict inequality constraints.

3 The adaptive reciprocity algorithm↩︎

To compute \(L_P(n)\) for a polytope of dimension \(d\), we need \(d+1\) evaluations at distinct integers. We always have \(L_P(0) = 1\) for free, and by 2 and 3 , each negative evaluation \((-1)^d L_P(-n) = K^{\mathrm{strict}}_{n\lambda, n\mu, n\mathbf{w}}\) reduces to a strict Kostka computation. The key observation is that strict Kostka coefficients are often zero for small \(n\): if any part of \(\mathbf{w}\) is 1, then strict GT patterns require the corresponding row to have distinct entries in a range that may be empty for small dilations. Even when the strict count is nonzero, its DP state space is typically much smaller than the ordinary Kostka DP for the same dilation.

This motivates an adaptive strategy:

The “cost” heuristic uses the most recent Kostka/strict-Kostka value as a proxy for the DP state-space size at the next dilation. The negative side typically wins for the first several points because strict Kostka equals zero.

4 Implementation↩︎

Our implementation is in Rust. The specialized research code used for the computations in this article is available as the kostka crate.¹ It contains the GT/Kostka dynamic programs, the Ehrhart interpolation and reciprocity routines, and the strict-count variants used below. The order-polytope experiments discussed later use the same reciprocity/interpolation viewpoint, but are conceptually separate from the lrcalc compatibility layer. At its core is a DP algorithm for (skew) Kostka coefficients \(K_{\lambda/\mu, \mathbf{w}}\) based on horizontal-strip removal, using BigUint arithmetic for exact computation. A variant with strict inequality constraints in the DP transitions handles the strict Kostka coefficients \(K^{\mathrm{strict}}_{\lambda/\mu, \mathbf{w}}\). The dimension of \(\mathrm{GT}(\lambda/\mu, \mathbf{w})\) is determined by counting free entries in the GT pattern with propagation of forced values. This gives the exact degree of the stretched Kostka polynomial, rather than only an a priori upper bound; see J. A. De Loera and T. B. McAllister [8] and Nzeutchap, Toumazet, and Butelle [4] for related degree computations. The Ehrhart polynomial is then recovered from evaluation points via Lagrange interpolation over \(\mathbb{Q}\).

The kostka CLI supports:

  kostka ehrhart --lambda 3,2,1 -w 1,1,1,1,1,1
  kostka ehrhart --lambda 4,3 --mu 2,1 -w 2,2,2
  kostka hstar --lambda 4,3,2,1 -w 1,1,1,1,1,1,1,1,1,1
  kostka hstar --lambda 4,3,1 --mu 2,1 -w 2,2,1

The tool also computes Kostka coefficients associated with flagged Schur functions, Littlewood–Richardson coefficients, and batch computations over all weights for a given shape. We refer to the Git repository for full documentation.

There is also a separate Rust implementation aimed at the classical lrcalc interface: lrcalc-rs.² That project is not the specialized Ehrhart engine used in Tables ¿tbl:tab:benchmarks? and 1; rather, it is a drop-in-oriented implementation of the liblrcalc ABI and command-line tools, with native code for Littlewood–Richardson coefficients, Schur products, skew Schur expansions, Schubert products, and fast Kostka special cases. The speedups in the present article come from exploiting the GT/Kostka and reciprocity structure directly, while lrcalc-rs packages related optimizations in a broader LR/Schur-function compatibility library.

5 Benchmarks↩︎

We compare our method against two alternatives: the \(\mathtt{OSCAR}\)/\(\mathtt{polymake}\) pipeline [9], [10] and A. Buch’s lrcalc [5], which reduces Kostka evaluations to LR coefficient computations. (\(\mathtt{OSCAR}\) is a Julia framework that calls \(\mathtt{polymake}\) for polyhedral computations, which in turn delegates lattice-point enumeration and \(h^*\)-vector computation to \(\mathtt{Normaliz}\) [11]; thus “\(\mathtt{OSCAR}\),” “\(\mathtt{polymake}\),” and “\(\mathtt{Normaliz}\)” refer to different layers of the same computational stack.) Tables ¿tbl:tab:benchmarks? and 1 summarize the results.

The \(\mathtt{OSCAR}\)/\(\mathtt{polymake}\) approach requires constructing the full vertex representation of the GT polytope—expensive in high dimensions—and does not exploit Ehrhart–Macdonald reciprocity. Our method avoids vertex enumeration entirely: each evaluation of \(L_P(n)\) or \(L^{*}_P(n)\) reduces to a Kostka coefficient DP whose complexity depends on the dilation \(n\) (not the number of vertices), and reciprocity provides many free evaluation points where strict Kostka equals zero.

6 The \(h^{*}\)-vector↩︎

For a rational polytope \(P\) with denominator \(q\) (the smallest positive integer such that \(qP\) is a lattice polytope), the lattice-point counting function \(L_P(n) = |nP \cap \mathbb{Z}^d|\) is a quasi-polynomialquasi-polynomial of degree \(d\) and period dividing \(q\). The Ehrhart series then satisfies \[\label{eq:hstar-rational} \sum_{n \geq 0} L_P(n) \, z^n = \frac{h^*_q(P;z)}{(1-z^q)^{d+1}},\tag{4}\] where \(h^*_q(P;z)\) is a polynomial of degree at most \(q(d+1)-1\) with non-negative integer coefficients, see [12] and [13]. Since GT polytopes have period collapse by Rassart’s theorem, we define the \(h^{*}\)-vector\(h^{*}\)-vector \((h^*_0, \dotsc, h^*_d)\) as in the lattice polytope case (\(q=1\)): \[\label{eq:hstar-lattice} \sum_{n \geq 0} L_P(n) \, z^n = \frac{h^*_0 + h^*_1 z + \dotsb + h^*_d z^d}{(1-z)^{d+1}}.\tag{5}\] Comparing 4 and 5 via the factorization \((1-z^q)^{d+1} = (1-z)^{d+1}(1+z+\dotsb+z^{q-1})^{d+1}\), we see that \(h^*_q(P;z) = (h^*_0 + \dotsb + h^*_d z^d) \cdot (1+z+\dotsb+z^{q-1})^{d+1}\). In particular, non-negativity of the \(h^*_i\) implies non-negativity of the coefficients of \(h^*_q(P;z)\). R. King, C. Tollu, and F. Toumazet [2] conjectured that the Ehrhart polynomial \(L_P(n)\) of a GT polytope has non-negative coefficients (in the standard monomial basis). Note that this is an independent property from \(h^*\)-positivity: neither implies the other in general. Our computations verify both properties for all straight-shape GT polytopes \(\mathrm{GT}(\lambda,(1^{|\lambda|}))\) with \(|\lambda| \leq 15\)—over 54,000 cases—providing further evidence for the conjecture. See also [14], [15] for related positivity conjectures for Ehrhart polynomials of polytopes related to key polynomials, flagged Schur functions, and cylindric Schur functions. For recent progress on Ehrhart positivity for marked order polytopes, relevant to skew GT polytopes without a weight condition, see [16].

As observed by C. Haase and T. B. McAllister [17], there exist rational polytopes with period collapse whose Ehrhart polynomial is not the Ehrhart polynomial of any lattice polytope. It is natural to ask whether the same phenomenon occurs for GT polytopes: is \(L_{\mathrm{GT}(\lambda/\mu, \mathbf{w})}(n)\) always realizable as the Ehrhart polynomial of a lattice polytope, or are GT polytopes genuinely different in this respect?

7 Order polytopes↩︎

The reciprocity approach extends naturally to order polytopesorder polytopes. For a finite poset \(P\) on \(\{0,1,\dotsc,n-1\}\), the order polytopeorder polytope [18] \[\mathcal{O}(P) = \bigl\{x \in [0,1]^n : x_a \leq x_b \text{ for all } a <_P b \bigr\}\] has the property that its lattice points at dilation \(t\) are weakly order-preserving maps \(P \to \{0,1,\dotsc,t\}\), so \(L_{\mathcal{O}(P)}(t) = \Omega(P, t+1)\), where \(\Omega(P,k)\) is Stanley’s order polynomial. Interior lattice points at dilation \(t\) correspond to strictly order-preserving maps \(P \to \{1,\dotsc,t-1\}\), giving Ehrhart–Macdonald reciprocity in the form \[\label{eq:order-reciprocity} L_{\mathcal{O}(P)}(-k) = (-1)^n \, \bar{\Omega}(P, k-1),\tag{6}\] where \(\bar{\Omega}(P,m)\) counts strict order-preserving maps \(P \to \{1,\dotsc,m\}\).

7.1 Frontier DP for order-preserving maps↩︎

The standard approach to evaluating \(\Omega(P,k)\) enumerates all order-preserving maps by backtracking over vertices, giving \(O(k^n)\) worst-case time. We introduce a frontier DPfrontier DP that processes vertices in a natural labeling order, maintaining as state only the values of live vertices—those with at least one unprocessed child. When a vertex \(v\) has no children among the remaining vertices, its value does not affect future constraints, so all valid choices are collapsed into a single multiplier \((k - \ell_v + 1)\), where \(\ell_v\) is the lower bound from parent constraints. States that agree on all live vertex values are merged.

The complexity is \(O(n \cdot k^{w})\), where \(w\) is the maximum frontier width (the largest number of simultaneously live vertices). For chains and zigzag (fence) posets, \(w = 1\); for antichains, \(w = 0\) (each vertex is independent); for Young diagram posets of shape \((\lambda_1, \dotsc, \lambda_r)\), \(w \leq \lambda_1\).

7.2 Order polytope benchmarks↩︎

We benchmark three methods for computing the \(h^*\)-vector of order polytopes:

1. Rust DP: our frontier DP with Ehrhart–Macdonald reciprocity and BigRational interpolation;

2. linext: Stanley’s formula—enumerate all linear extensions \(\mathcal{L}(P)\) and tabulate descents (streaming, \(O(n)\) memory);

3. \(\mathtt{polymake}\): the \(\mathtt{polymake}\)/\(\mathtt{Normaliz}\) backend for general lattice-point enumeration.

By Stanley’s theorem [19], for a naturally labeled poset \(P\), the \(h^*\)-vector of \(\mathcal{O}(P)\) equals the \(P\)-Eulerian polynomial: \[h^*_i = \#\{\sigma \in \mathcal{L}(P) : \mathrm{des}(\sigma) = i\},\] where \(\mathcal{L}(P)\) is the set of linear extensions and \(\mathrm{des}(\sigma)\) counts descents. Counting linear extensions is \(\#\mathsf{P}\)-complete in general [20], and even for structured posets the number \(|\mathcal{L}(P)|\) can be enormous—e.g., \(|\mathcal{L}(P)| = n!\) for the antichain on \(n\) elements. Our Ehrhart/reciprocity approach avoids linear extension enumeration entirely: it requires only \(\sim n/2\) evaluations of the order polynomial via the frontier DP, each costing \(O(n \cdot k^w)\).

7.3 Non-real-rooted \(h^*\)-vectors from permutation posets↩︎

As an application of the fast \(h^*\)-computation, we investigate real-rootedness of \(h^*\)-polynomials for permutation posetspermutation posets. Given \(w \in S_n\), the permutation poset \(P_w\) is the partial order on \(\{1,\dotsc,n\}\) with \(i <_{P_w} j\) iff \(i < j\) and \(w_i < w_j\). Since \(P_w\) is always naturally labeled, \(h^*(\mathcal{O}(P_w); t)\) equals the \(P\)-Eulerian polynomial by Stanley’s theorem.

Starting from a \(321\)-avoiding permutation of Stembridge, \[w_0 = (2,4,6,8,10,1,12,3,15,5,17,7,9,11,13,14,16) \in S_{17},\] whose permutation poset has \(h^*(t) = 1 + 32t + 336t^2 + 1420t^3 + 2534t^4 + 1946t^5 + 658t^6 + 86t^7 + 3t^8\) (ultra-log-concave but not real-rooted), we searched over 34,226 nearby \(4321\)-avoiding permutations (within 3 arbitrary transpositions) and found exactly 22 with non-real-rooted \(h^*\). All 22 counterexamples are ultra-log-concave, hence also log-concave.

8 Further directions↩︎

Our method applies whenever one has a combinatorial interpretation of both \(L_P(n)\) and the interior count \(L^{*}_P(n)\) that admits a DP or transfer-matrix algorithm, together with a regime where \(L^{*}_P(n) = 0\) for small \(n\). GT polytopes are closely related to flow polytopesflow polytopes on directed acyclic graphs [21].

Another natural test case is the Birkhoff polytope \(B_\ell\). Here \(L_{B_\ell}(n)\) and \(L^{*}_{B_\ell}(n)\) count nonnegative and positive integer magic squares with line sum \(n\), so the reciprocity setup again applies. For Birkhoff polytopes, however, the symmetry-based interpolation trick is classical rather than new. Beck and Pixton [22] record the standard consequence of Ehrhart–Macdonald reciprocity: from reciprocity one gets \(L_{B_\ell}(-\ell-t)=(-1)^{\ell-1}L_{B_\ell}(t)\) and \(L_{B_\ell}(-1)=\cdots=L_{B_\ell}(-\ell+1)=0\). They explicitly describe the resulting strategy: compute the first \(\binom{\ell-1}{2}\) values, use the symmetry and trivial zeros, and interpolate the polynomial. They also note that the volume of \(B_8\) had been computed by essentially this method, combined with additional computational tricks. Beck–Pixton’s own main contribution in that paper is a different residue-based method for computing \(H_\ell(t)\) and the volumes of larger Birkhoff polytopes. For context, our implementation computes the Ehrhart polynomial in under 10 ms for \(B_3\) and \(B_4\), in 0.28 s for \(B_5\), and in about 265 s for \(B_6\). Moreover, \(L^{*}_{B_\ell}(n)=0\) for \(1 \leq n < \ell\), so the first \(\ell-1\) negative evaluation points are free. Thus the new contribution in this direction would be the GT-DP implementation and its adaptation of that classical Birkhoff symmetry trick, not the trick itself. This suggests that the method is viable for smaller Birkhoff polytopes, although the growth by \(\ell=6\) is already steep.

Matroid polytopes are another natural family where the reciprocity approach may be effective.

Acknowledgements↩︎

Some Rust code, installation work, and benchmark scripts were developed with assistance from OpenAI Codex.

References↩︎