update pet for support for recent versions of clang
[barvinok.git] / vertex_cone.cc
blob0c79f3205f65b63fe599e0ddf44f7f8426635996
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"
9 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
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.
15 static int pos_cmp(const void *a, const void *b)
17 int pa, pb;
18 const Vector *va = *(const Vector **)a;
19 const Vector *vb = *(const Vector **)b;
21 pa = First_Non_Zero(va->p, va->Size);
22 pb = First_Non_Zero(vb->p, vb->Size);
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);
29 return pa - pb;
32 vertex_cone::vertex_cone(unsigned dim) : dim(dim)
34 E_vertex = new evalue *[dim];
35 bernoulli_t = new evalue **[dim];
37 assert(dim > 0);
39 coeff = ALLOCN(Vector *, dim);
40 for (int i = 0; i < dim; ++i)
41 coeff[i] = Vector_Alloc(dim);
43 Rays.NbRows = Rays.NbColumns = dim;
44 Rays.p = ALLOCN(Value *, dim);
46 coeff_power = new struct power **[dim];
47 for (int i = 0; i < dim; ++i)
48 coeff_power[i] = new struct power *[dim];
51 void vertex_cone::init(const mat_ZZ &rays, Param_Vertices *V,
52 unsigned max_power)
54 unsigned nparam = V->Vertex->NbColumns - 2;
55 this->max_power = max_power;
57 for (int i = 0; i < dim; ++i)
58 zz2values(rays[i], coeff[i]->p);
59 qsort(coeff, dim, sizeof(Vector *), pos_cmp);
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;
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);
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);
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;
88 * Returns b(n, E_vertex[i])
90 const evalue *vertex_cone::get_bernoulli(int n, int i)
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]);
98 return bernoulli_t[i][n-1];
101 void vertex_cone::clear()
103 for (int i = 0; i < dim; ++i)
104 if (E_vertex[i])
105 evalue_free(E_vertex[i]);
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]);
112 delete [] bernoulli_t[i];
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];