vertex_cone.cc

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