vertex_cone.cc

   1 #include <assert.h>
   2 #include <NTL/mat_ZZ.h>
   3 #include <barvinok/evalue.h>
   4 #include "bernoulli.h"
   5 #include "conversion.h"
   6 #include "lattice_point.h"
   7 #include "vertex_cone.h"
   8
   9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  10
  11 /* Compares the position of the first non-zero element
  12  * in both vectors and the second non-zero element
  13  * if the position of the first non-zero element is the same.
  14  */
  15 static int pos_cmp(const void *a, const void *b)
  16 {
  17     int pa, pb;
  18     const Vector *va = *(const Vector **)a;
  19     const Vector *vb = *(const Vector **)b;
  20
  21     pa = First_Non_Zero(va->p, va->Size);
  22     pb = First_Non_Zero(vb->p, vb->Size);
  23
  24     if (pa == pb) {
  25         pa = First_Non_Zero(va->p+pa+1, va->Size-pa-1);
  26         pb = First_Non_Zero(vb->p+pa+1, vb->Size-pa-1);
  27     }
  28
  29     return pa - pb;
  30 }
  31
  32 vertex_cone::vertex_cone(unsigned dim) : dim(dim)
  33 {
  34     E_vertex = new evalue *[dim];
  35     bernoulli_t = new evalue **[dim];
  36
  37     assert(dim > 0);
  38
  39     coeff = ALLOCN(Vector *, dim);
  40     for (int i = 0; i < dim; ++i)
  41         coeff[i] = Vector_Alloc(dim);
  42
  43     Rays.NbRows = Rays.NbColumns = dim;
  44     Rays.p = ALLOCN(Value *, dim);
  45
  46     coeff_power = new struct power **[dim];
  47     for (int i = 0; i < dim; ++i)
  48         coeff_power[i] = new struct power *[dim];
  49 }
  50
  51 void vertex_cone::init(const mat_ZZ &rays, Param_Vertices *V,
  52                         unsigned max_power)
  53 {
  54     unsigned nparam = V->Vertex->NbColumns - 2;
  55     this->max_power = max_power;
  56
  57     for (int i = 0; i < dim; ++i)
  58         zz2values(rays[i], coeff[i]->p);
  59     qsort(coeff, dim, sizeof(Vector *), pos_cmp);
  60
  61     for (int i = 0; i < dim; ++i) {
  62         for (int j = 0; j < dim; ++j) {
  63             if (value_notzero_p(coeff[i]->p[j]))
  64                 coeff_power[i][j] = new struct power(coeff[i]->p[j], max_power);
  65             else
  66                 coeff_power[i][j] = NULL;
  67         }
  68     }
  69
  70     for (int i = 0; i < dim; ++i)
  71         Rays.p[i] = coeff[i]->p;
  72     Matrix *T = Transpose(&Rays);
  73     Matrix *L = relative_coordinates(V, T);
  74     Matrix_Free(T);
  75
  76     for (int i = 0; i < dim; ++i)
  77         E_vertex[i] = ceiling(L->p[i], V->Vertex->p[0][nparam+1], nparam, NULL);
  78     Matrix_Free(L);
  79
  80     for (int j = 0; j < dim; ++j) {
  81         bernoulli_t[j] = new evalue *[max_power];
  82         for (int k = 0; k < max_power; ++k)
  83             bernoulli_t[j][k] = NULL;
  84     }
  85 }
  86
  87 /*
  88  * Returns b(n, E_vertex[i])
  89  */
  90 const evalue *vertex_cone::get_bernoulli(int n, int i)
  91 {
  92     assert(n > 0);
  93     if (!bernoulli_t[i][n-1]) {
  94         struct poly_list *bernoulli = bernoulli_compute(n);
  95         bernoulli_t[i][n-1] = evalue_polynomial(bernoulli->poly[n],
  96                                                 E_vertex[i]);
  97     }
  98     return bernoulli_t[i][n-1];
  99 }
 100
 101 void vertex_cone::clear()
 102 {
 103     for (int i = 0; i < dim; ++i)
 104         if (E_vertex[i])
 105             evalue_free(E_vertex[i]);
 106
 107     for (int i = 0; i < dim; ++i) {
 108         for (int j = 0; j < max_power; ++j) {
 109             if (bernoulli_t[i][j])
 110                 evalue_free(bernoulli_t[i][j]);
 111         }
 112         delete [] bernoulli_t[i];
 113     }
 114
 115     for (int i = 0; i < dim; ++i) {
 116         for (int j = 0; j < dim; ++j)
 117             if (coeff_power[i][j])
 118                 delete coeff_power[i][j];
 119     }
 120 }