polymake/h_star_vector.cc

   1 #include <gmp.h>
   2 #include <cstdlib>
   3 #include <Rational.h>
   4 #include <Poly.h>
   5 #include <Matrix.h>
   6 #include <Vector.h>
   7 #include <barvinok/barvinok.h>
   8 #include <barvinok/util.h>
   9 #include "convert.h"
  10
  11 #define MAXRAYS     0
  12
  13 namespace polymake { namespace polytope {
  14
  15 /*
  16  * The h^*-vector and the Ehrhart polynomial f(t) of a lattice polytope are
  17  * related through
  18  *
  19  *              f(t) = \sum_{i=0}^d h^*_i { t+d-i \choose d }
  20  *
  21  * I.e.,
  22  *              c = M h^*
  23  *
  24  * with c and h^* column vectors,
  25  *              f(t) = \sum_{i=0}^d c_i t^i,
  26  * and column i of M the coefficients of the polynomial { t+d-i \choose d }
  27  *
  28  * Since it is actually easier for us to compute the Ehrhart polynomial,
  29  * we compute the h^*-vector from this Ehrhart polynomial as
  30  *
  31  *              h^* = M^{-1} c
  32  *
  33  * All coefficients are ordered from index 0 to d.
  34  */
  35 void h_star_vector(Poly& p)
  36 {
  37     Matrix<Rational> F = p.give("FACETS | INEQUALITIES");
  38     ::Matrix *M = polymake_constraints2polylib(F);
  39     Polyhedron *A = Constraints2Polyhedron(M, 0);
  40     Matrix_Free(M);
  41     for (unsigned i = 0; i < A->NbRays; ++i) {
  42         if (value_zero_p(A->Ray[i][0]) || value_zero_p(A->Ray[i][1+A->Dimension])) {
  43             Polyhedron_Free(A);
  44             throw std::runtime_error("polyhedron must be bounded");
  45         }
  46         if (value_notone_p(A->Ray[i][1+A->Dimension])) {
  47             Polyhedron_Free(A);
  48             throw std::runtime_error("polytope must have integer vertices");
  49         }
  50     }
  51     Polyhedron *C = Cone_over_Polyhedron(A);
  52     Polyhedron_Free(A);
  53     Polyhedron *U = Universe_Polyhedron(1);
  54     evalue *EP = barvinok_enumerate_ev(C, U, MAXRAYS);
  55     Polyhedron_Free(C);
  56     Polyhedron_Free(U);
  57
  58     assert(value_zero_p(EP->d));
  59     assert(EP->x.p->type == partition);
  60     assert(EP->x.p->size == 2);
  61     evalue *poly = &EP->x.p->arr[1];
  62     if (value_notzero_p(poly->d)) {
  63         assert(value_one_p(poly->d));
  64         p.take("H_STAR_VECTOR") << Integer(poly->x.n);
  65     } else {
  66         assert(poly->x.p->type == polynomial);
  67         int d = poly->x.p->size - 1;
  68         ::Matrix *M = Matrix_Alloc(d+2, d+2);
  69         ::Vector *V = Vector_Alloc(d+2);
  70         ::Vector *h = Vector_Alloc(d+2);
  71
  72         Value one;
  73         value_init(one);
  74         value_set_si(one, 1);
  75         Value k;
  76         value_init(k);
  77         for (int i = 0; i <= d; ++i) {
  78             Vector_Set(V->p, 0, d+2);
  79             value_set_si(V->p[d], 1);
  80             for (int j = d; j > 0; --j) {
  81                 value_set_si(k, j-i);
  82                 Vector_Combine(V->p + j-1, V->p + j, V->p + j-1, k, one, d-j+2);
  83             }
  84             for (int j = 0; j <= d; ++j)
  85                 value_assign(M->p[j][i], V->p[d-j]);
  86         }
  87         value_clear(one);
  88         value_clear(k);
  89
  90         Factorial(d, &M->p[d+1][d+1]);
  91         ::Matrix *inv = Matrix_Alloc(d+2, d+2);
  92         int ok = Matrix_Inverse(M, inv);
  93         Matrix_Free(M);
  94         assert(ok);
  95         for (int i = 0; i <= d; ++i) {
  96             assert(value_one_p(poly->x.p->arr[i].d));
  97             value_assign(V->p[i], poly->x.p->arr[i].x.n);
  98         }
  99         value_set_si(V->p[d+1], 1);
 100         Matrix_Vector_Product(inv, V->p, h->p);
 101         Vector_Free(V);
 102         Matrix_Free(inv);
 103         assert(value_one_p(h->p[d+1]));
 104
 105         Vector<Integer> h_star(d+1);
 106         for (int i = 0; i <= d; ++i)
 107             h_star[i] = Integer(h->p[i]);
 108         Vector_Free(h);
 109         p.take("H_STAR_VECTOR") << h_star;
 110     }
 111     free_evalue_refs(EP);
 112     free(EP);
 113 }
 114
 115 } }
 116
 117 using namespace polymake;
 118
 119 int main(int argc, const char *argv[])
 120 {
 121     if (argc != 2) {
 122         return 1;
 123     }
 124     try {
 125         Poly p(argv[1], ios::in | ios::out);
 126         polytope::h_star_vector(p);
 127     } catch (const std::exception& e) {
 128         cerr << e.what() << endl;
 129         return 1;
 130    }
 131     return 0;
 132 }