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