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     assert(dim > 0);
  35
  36     coeff = ALLOCN(Vector *, dim);
  37     for (int i = 0; i < dim; ++i)
  38         coeff[i] = Vector_Alloc(dim);
  39
  40     Rays.NbRows = Rays.NbColumns = dim;
  41     Rays.p = ALLOCN(Value *, dim);
  42
  43     coeff_power = new struct power **[dim];
  44     for (int i = 0; i < dim; ++i)
  45         coeff_power[i] = new struct power *[dim];
  46 }
  47
  48 void vertex_cone::init(const mat_ZZ &rays, Param_Vertices *V,
  49                         unsigned max_power)
  50 {
  51     unsigned nparam = V->Vertex->NbColumns - 2;
  52     this->max_power = max_power;
  53
  54     for (int i = 0; i < dim; ++i)
  55         zz2values(rays[i], coeff[i]->p);
  56     qsort(coeff, dim, sizeof(Vector *), pos_cmp);
  57
  58     for (int i = 0; i < dim; ++i) {
  59         for (int j = 0; j < dim; ++j) {
  60             if (value_notzero_p(coeff[i]->p[j]))
  61                 coeff_power[i][j] = new struct power(coeff[i]->p[j], max_power);
  62             else
  63                 coeff_power[i][j] = NULL;
  64         }
  65     }
  66
  67     for (int i = 0; i < dim; ++i)
  68         Rays.p[i] = coeff[i]->p;
  69     Matrix *T = Transpose(&Rays);
  70     Matrix *L = relative_coordinates(V, T);
  71     Matrix_Free(T);
  72
  73     for (int i = 0; i < dim; ++i)
  74         E_vertex[i] = ceiling(L->p[i], V->Vertex->p[0][nparam+1], nparam, NULL);
  75     Matrix_Free(L);
  76
  77     for (int j = 0; j < dim; ++j) {
  78         bernoulli_t[j] = new evalue *[max_power];
  79         for (int k = 0; k < max_power; ++k)
  80             bernoulli_t[j][k] = NULL;
  81     }
  82 }
  83
  84 /*
  85  * Returns b(n, E_vertex[i])
  86  */
  87 const evalue *vertex_cone::get_bernoulli(int n, int i)
  88 {
  89     assert(n > 0);
  90     if (!bernoulli_t[i][n-1]) {
  91         struct poly_list *bernoulli = bernoulli_compute(n);
  92         bernoulli_t[i][n-1] = evalue_polynomial(bernoulli->poly[n],
  93                                                 E_vertex[i]);
  94     }
  95     return bernoulli_t[i][n-1];
  96 }
  97
  98 void vertex_cone::clear()
  99 {
 100     for (int i = 0; i < dim; ++i)
 101         if (E_vertex[i])
 102             evalue_free(E_vertex[i]);
 103
 104     for (int i = 0; i < dim; ++i) {
 105         for (int j = 0; j < max_power; ++j) {
 106             if (bernoulli_t[i][j])
 107                 evalue_free(bernoulli_t[i][j]);
 108         }
 109         delete [] bernoulli_t[i];
 110     }
 111
 112     for (int i = 0; i < dim; ++i) {
 113         for (int j = 0; j < dim; ++j)
 114             if (coeff_power[i][j])
 115                 delete coeff_power[i][j];
 116     }
 117 }