iscc: add "source" operation
[barvinok.git] / vertex_cone.cc
blob74f6bcdcd719e86962c1f22e117398cbabc7e964
1 #include "bernoulli.h"
2 #include "conversion.h"
3 #include "lattice_point.h"
4 #include "vertex_cone.h"
6 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
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.
12 static int pos_cmp(const void *a, const void *b)
14 int pa, pb;
15 const Vector *va = *(const Vector **)a;
16 const Vector *vb = *(const Vector **)b;
18 pa = First_Non_Zero(va->p, va->Size);
19 pb = First_Non_Zero(vb->p, vb->Size);
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);
26 return pa - pb;
29 vertex_cone::vertex_cone(unsigned dim) : dim(dim)
31 E_vertex = new evalue *[dim];
32 bernoulli_t = new evalue **[dim];
34 coeff = ALLOCN(Vector *, dim);
35 for (int i = 0; i < dim; ++i)
36 coeff[i] = Vector_Alloc(dim);
38 Rays.NbRows = Rays.NbColumns = dim;
39 Rays.p = ALLOCN(Value *, dim);
41 coeff_power = new struct power **[dim];
42 for (int i = 0; i < dim; ++i)
43 coeff_power[i] = new struct power *[dim];
46 void vertex_cone::init(const mat_ZZ &rays, Param_Vertices *V,
47 unsigned max_power)
49 unsigned nparam = V->Vertex->NbColumns - 2;
50 this->max_power = max_power;
52 for (int i = 0; i < dim; ++i)
53 zz2values(rays[i], coeff[i]->p);
54 qsort(coeff, dim, sizeof(Vector *), pos_cmp);
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;
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);
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);
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;
83 * Returns b(n, E_vertex[i])
85 const evalue *vertex_cone::get_bernoulli(int n, int i)
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]);
93 return bernoulli_t[i][n-1];
96 void vertex_cone::clear()
98 for (int i = 0; i < dim; ++i)
99 if (E_vertex[i])
100 evalue_free(E_vertex[i]);
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]);
107 delete [] bernoulli_t[i];
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];