1 | initial version |
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.
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.