polymake/h_star_vector.cc

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