Ask Your Question

Revision history [back]

Sage doesn't have a direct 'print out the h* vector' command. It's not too difficult to manipulate the Ehrhart polynomial to get it, though I'm not sure my way could be cleanly implemented into Sage proper.

Let's say P is your polytope, of dimension n=P.dimension(), with Ehrhart polynomial p=P.ehrhart_polynomial() . Also, let t=p.parent().gen() . We need a polynomial ring to work over, and this let's us work over the same polynomial ring as p. Now, compute the first n+2 terms of the Ehrhart series (or more, if you want to be safe) X=0 for i in range(0,n+10): X+=p(i)*t^i print(X)

Now, look at X(1-t)^(n+1) . If X were the entire infinite Ehrhart series, we would get exactly the h vector. But since we took sufficiently many terms of the Ehrhart series, we'll get lower order terms corresponding to the h* vector, and higher order error terms that can be discarded.

click to hide/show revision 2
No.2 Revision

Sage doesn't have a direct 'print out the h* vector' command. It's not too difficult to manipulate the Ehrhart polynomial to get it, though I'm not sure my way could be cleanly implemented into Sage proper.

Let's say P P is your polytope, of dimension n=P.dimension(), n=P.dimension(), with Ehrhart polynomial p=P.ehrhart_polynomial() p=P.ehrhart_polynomial() . Also, let t=p.parent().gen() t=p.parent().gen() . We need a polynomial ring to work over, and this let's us work over the same polynomial ring as p. Now, compute the first n+2 terms of the Ehrhart series (or more, if you want to be safe) safe)

X=0
for i in range(0,n+10):
    X+=p(i)*t^i
print(X)

print(X)

Now, look at X(1-t)^(n+1) X*(1-t)^(n+1) . If X were the entire infinite Ehrhart series, we would get exactly the h h* vector. But since we took sufficiently many terms of the Ehrhart series, we'll get lower order terms corresponding to the h* vector, and higher order error terms that can be discarded.