barvinok 0.41.8
[barvinok.git] / edomain.cc
blob1852f7936645f9db86f4625ab853a2abac54566e
1 #include <assert.h>
2 #include <sstream>
3 //#include "fdstream.h"
4 #include <barvinok/util.h>
5 #include <barvinok/sample.h>
6 #include <barvinok/barvinok.h>
7 #include "edomain.h"
8 #include "evalue_util.h"
10 using std::vector;
11 using std::endl;
12 using std::ostream;
14 EDomain *EDomain::new_from_ge_constraint(ge_constraint *ge, int sign,
15 barvinok_options *options)
17 if (ge->simplified && sign == 0)
18 return NULL;
20 EDomain *ED;
21 Matrix *M2;
22 M2 = Matrix_Copy(ge->M);
23 if (sign == 1) {
24 if (!ge->simplified)
25 value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
26 M2->p[M2->NbRows-1][M2->NbColumns-1]);
27 } else if (sign == -1) {
28 Value mone;
29 value_init(mone);
30 value_set_si(mone, -1);
31 Vector_Scale(M2->p[M2->NbRows-1]+1, M2->p[M2->NbRows-1]+1,
32 mone, M2->NbColumns-1);
33 value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
34 M2->p[M2->NbRows-1][M2->NbColumns-1]);
35 value_clear(mone);
36 } else {
37 value_set_si(M2->p[M2->NbRows-1][0], 0);
39 Polyhedron *D2 = Constraints2Polyhedron(M2, options->MaxRays);
40 ED = new EDomain(D2, ge->D, ge->new_floors);
41 Polyhedron_Free(D2);
42 Matrix_Free(M2);
43 return ED;
46 static void print_term(ostream& os, Value v, int pos, int dim,
47 char **names, int *first)
49 if (value_zero_p(v)) {
50 if (first && *first && pos >= dim)
51 os << "0";
52 return;
55 if (first) {
56 if (!*first && value_pos_p(v))
57 os << "+";
58 *first = 0;
60 if (pos < dim) {
61 if (value_mone_p(v)) {
62 os << "-";
63 } else if (!value_one_p(v))
64 os << VALUE_TO_INT(v);
65 os << names[pos];
66 } else
67 os << VALUE_TO_INT(v);
70 void EDomain_floor::print(ostream& os, char **p) const
72 int first = 1;
73 os << "[";
74 os << "(";
75 for (int i = 0; i < v->Size-2; ++i)
76 print_term(os, v->p[1+i], i, v->Size-2, p, &first);
77 os << ")";
78 os << "/";
79 print_term(os, v->p[0], v->Size-2, v->Size-2, p, NULL);
80 os << "]";
83 void EDomain::print_constraints(ostream& os, char **p,
84 barvinok_options *options) const
86 Value tmp;
87 value_init(tmp);
89 Matrix *M = Matrix_Alloc(2*floors.size(), 1+D->Dimension+1);
90 value_set_si(tmp, -1);
91 for (int i = 0; i < floors.size(); ++i) {
92 value_set_si(M->p[2*i][0], 1);
93 Vector_Copy(floors[i]->v->p+1, M->p[2*i]+1, dimension());
94 value_assign(M->p[2*i][1+D->Dimension], floors[i]->v->p[1+dimension()]);
95 value_oppose(M->p[2*i][1+dimension()+i], floors[i]->v->p[0]);
97 Vector_Scale(M->p[2*i]+1, M->p[2*i+1]+1, tmp, D->Dimension+1);
98 value_addto(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension],
99 M->p[2*i+1][1+dimension()+i]);
100 value_decrement(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension]);
101 value_set_si(M->p[2*i+1][0], 1);
103 Polyhedron *E = Constraints2Polyhedron(M, options->MaxRays);
104 Matrix_Free(M);
105 Polyhedron *SD = DomainSimplify(D, E, options->MaxRays);
106 Polyhedron_Free(E);
108 char **names = p;
109 unsigned dim = dimension();
110 if (dim < SD->Dimension) {
111 names = new char * [SD->Dimension];
112 int i;
113 for (i = 0; i < dim; ++i)
114 names[i] = p[i];
115 for ( ; i < SD->Dimension; ++i) {
116 std::ostringstream strm;
117 floors[i-dim]->print(strm, p);
118 names[i] = strdup(strm.str().c_str());
122 for (int i = 0; i < SD->NbConstraints; ++i) {
123 int first = 1;
124 int v = First_Non_Zero(SD->Constraint[i]+1, SD->Dimension);
125 if (v == -1)
126 continue;
127 if (i)
128 os << " && ";
129 if (value_pos_p(SD->Constraint[i][v+1])) {
130 print_term(os, SD->Constraint[i][v+1], v, SD->Dimension,
131 names, NULL);
132 if (value_zero_p(SD->Constraint[i][0]))
133 os << " = ";
134 else
135 os << " >= ";
136 for (int j = v+1; j <= SD->Dimension; ++j) {
137 value_oppose(tmp, SD->Constraint[i][1+j]);
138 print_term(os, tmp, j, SD->Dimension,
139 names, &first);
141 } else {
142 value_oppose(tmp, SD->Constraint[i][1+v]);
143 print_term(os, tmp, v, SD->Dimension,
144 names, NULL);
145 if (value_zero_p(SD->Constraint[i][0]))
146 os << " = ";
147 else
148 os << " <= ";
149 for (int j = v+1; j <= SD->Dimension; ++j)
150 print_term(os, SD->Constraint[i][1+j], j, SD->Dimension,
151 names, &first);
154 value_clear(tmp);
155 Domain_Free(SD);
157 if (dim < D->Dimension) {
158 for (int i = dim; i < D->Dimension; ++i)
159 free(names[i]);
160 delete [] names;
165 void EDomain::print(FILE *out, char **p)
167 fdostream os(dup(fileno(out)));
168 for (int i = 0; i < floors.size(); ++i) {
169 os << "floor " << i << ": [";
170 evalue_print(os, floors[i]->e, p);
171 os << "]" << endl;
173 Polyhedron_Print(out, P_VALUE_FMT, D);
177 static int type_offset(enode *p)
179 return p->type == fractional ? 1 :
180 p->type == flooring ? 1 : 0;
183 static void add_coeff(Value *cons, int len, evalue *coeff, int pos)
185 Value tmp;
187 assert(value_notzero_p(coeff->d));
189 value_init(tmp);
191 value_lcm(tmp, cons[0], coeff->d);
192 value_division(tmp, tmp, cons[0]);
193 Vector_Scale(cons, cons, tmp, len);
194 value_division(tmp, cons[0], coeff->d);
195 value_addmul(cons[pos], tmp, coeff->x.n);
197 value_clear(tmp);
200 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len);
202 static void add_fract(evalue *e, Value *cons, int len, int pos)
204 evalue mone;
205 value_init(mone.d);
206 evalue_set_si(&mone, -1, 1);
208 /* contribution of alpha * fract(X) is
209 * alpha * X ...
211 assert(e->x.p->size == 3);
212 evalue argument;
213 value_init(argument.d);
214 evalue_copy(&argument, &e->x.p->arr[0]);
215 emul(&e->x.p->arr[2], &argument);
216 evalue2constraint_r(NULL, &argument, cons, len);
217 free_evalue_refs(&argument);
219 /* - alpha * floor(X) */
220 emul(&mone, &e->x.p->arr[2]);
221 add_coeff(cons, len, &e->x.p->arr[2], pos);
222 emul(&mone, &e->x.p->arr[2]);
224 free_evalue_refs(&mone);
227 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len)
229 int r = 0;
230 if (value_zero_p(E->d) && E->x.p->type == fractional) {
231 int i;
232 assert(E->x.p->size == 3);
233 r = evalue2constraint_r(D, &E->x.p->arr[1], cons, len);
234 assert(value_notzero_p(E->x.p->arr[2].d));
235 if (D && (i = D->find_floor(&E->x.p->arr[0])) >= 0) {
236 add_fract(E, cons, len, 1+D->D->Dimension-D->floors.size()+i);
237 } else {
238 if (value_pos_p(E->x.p->arr[2].x.n)) {
239 evalue coeff;
240 value_init(coeff.d);
241 value_init(coeff.x.n);
242 value_set_si(coeff.d, 1);
243 evalue_denom(&E->x.p->arr[0], &coeff.d);
244 value_decrement(coeff.x.n, coeff.d);
245 emul(&E->x.p->arr[2], &coeff);
246 add_coeff(cons, len, &coeff, len-1);
247 free_evalue_refs(&coeff);
249 r = 1;
251 } else if (value_zero_p(E->d)) {
252 assert(E->x.p->type == polynomial);
253 assert(E->x.p->size == 2);
254 r = evalue2constraint_r(D, &E->x.p->arr[0], cons, len) || r;
255 add_coeff(cons, len, &E->x.p->arr[1], E->x.p->pos);
256 } else {
257 add_coeff(cons, len, E, len-1);
259 return r;
262 EDomain_floor::EDomain_floor(const evalue *f, int dim)
264 e = new evalue;
265 value_init(e->d);
266 evalue_copy(e, f);
267 v = Vector_Alloc(1+dim+1);
268 value_set_si(v->p[0], 1);
269 evalue2constraint_r(NULL, e, v->p, v->Size);
270 refcount = 1;
271 substituted = false;
274 void EDomain_floor::eval(Value *values, Value *res) const
276 Inner_Product(v->p+1, values, v->Size-2, res);
277 value_addto(*res, *res, v->p[v->Size-1]);
278 value_pdivision(*res, *res, v->p[0]);
281 int evalue2constraint(EDomain *D, evalue *E, Value *cons, int len)
283 Vector_Set(cons, 0, len);
284 value_set_si(cons[0], 1);
285 return evalue2constraint_r(D, E, cons, len);
288 bool EDomain::contains(Value *point, int len) const
290 assert(len <= D->Dimension);
291 if (len == D->Dimension)
292 return in_domain(D, point);
294 Vector *all_val = Vector_Alloc(D->Dimension);
295 Vector_Copy(point, all_val->p, len);
296 for (int i = len-dimension(); i < floors.size(); ++i)
297 floors[i]->eval(all_val->p, &all_val->p[dimension()+i]);
298 bool in = in_domain(D, all_val->p);
299 Vector_Free(all_val);
300 return in;
303 ge_constraint *EDomain::compute_ge_constraint(evalue *constraint) const
305 ge_constraint *ge = new ge_constraint(this);
306 evalue mone;
307 value_init(mone.d);
308 evalue_set_si(&mone, -1, 1);
309 int fract = 0;
310 for (evalue *e = constraint; value_zero_p(e->d);
311 e = &e->x.p->arr[type_offset(e->x.p)]) {
312 if (e->x.p->type != fractional)
313 continue;
314 if (find_floor(&e->x.p->arr[0]) == -1)
315 ++fract;
318 int rows = D->NbConstraints+2*fract+1;
319 int cols = 2+D->Dimension+fract;
320 ge->M = Matrix_Alloc(rows, cols);
321 for (int i = 0; i < D->NbConstraints; ++i) {
322 Vector_Copy(D->Constraint[i], ge->M->p[i], 1+D->Dimension);
323 value_assign(ge->M->p[i][1+D->Dimension+fract],
324 D->Constraint[i][1+D->Dimension]);
326 value_set_si(ge->M->p[rows-1][0], 1);
327 fract = 0;
328 evalue *e;
329 for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
330 if (e->x.p->type == fractional) {
331 int i, pos;
333 i = find_floor(&e->x.p->arr[0]);
334 if (i >= 0)
335 pos = D->Dimension-floors.size()+i;
336 else
337 pos = D->Dimension+fract;
339 add_fract(e, ge->M->p[rows-1], cols, 1+pos);
341 if (pos < D->Dimension)
342 continue;
344 EDomain_floor *new_floor;
345 new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
347 /* constraints for the new floor */
348 int row = D->NbConstraints+2*fract;
349 Vector_Copy(new_floor->v->p+1, ge->M->p[row]+1, dimension());
350 value_assign(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]);
351 value_oppose(ge->M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
352 value_set_si(ge->M->p[row][0], 1);
353 assert(value_eq(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]));
354 assert(Vector_Equal(new_floor->v->p+1, ge->M->p[row]+1, dimension()));
356 Vector_Scale(ge->M->p[row]+1, ge->M->p[row+1]+1, mone.x.n, cols-1);
357 value_set_si(ge->M->p[row+1][0], 1);
358 value_addto(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1],
359 ge->M->p[row+1][1+D->Dimension+fract]);
360 value_decrement(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1]);
362 ge->new_floors.push_back(new_floor);
364 ++fract;
365 } else {
366 assert(e->x.p->type == polynomial);
367 assert(e->x.p->size == 2);
368 add_coeff(ge->M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
371 add_coeff(ge->M->p[rows-1], cols, e, cols-1);
372 ge->simplified = ConstraintSimplify(ge->M->p[rows-1], ge->M->p[rows-1],
373 cols, &ge->M->p[rows-1][0]);
374 value_set_si(ge->M->p[rows-1][0], 1);
375 free_evalue_refs(&mone);
376 return ge;
379 /* "align" matrix to have nrows by inserting
380 * the necessary number of rows and an equal number of columns at the end
381 * right before the constant row/column
383 static Matrix *align_matrix_initial(Matrix *M, int nrows)
385 int i;
386 int newrows = nrows - M->NbRows;
387 Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
388 for (i = 0; i < M->NbRows-1; ++i) {
389 Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
390 value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
392 for (i = 0; i <= newrows; ++i)
393 value_assign(M2->p[M->NbRows-1+i][M->NbColumns-1+i],
394 M->p[M->NbRows-1][M->NbColumns-1]);
395 return M2;
398 static Matrix *InsertColumns(Matrix *M, int pos, int n)
400 Matrix *R;
401 int i;
403 R = Matrix_Alloc(M->NbRows, M->NbColumns+n);
404 for (i = 0; i < M->NbRows; ++i) {
405 Vector_Copy(M->p[i], R->p[i], pos);
406 Vector_Copy(M->p[i]+pos, R->p[i]+pos+n, M->NbColumns-pos);
408 return R;
411 void EDomain_floor::substitute(evalue **subs, Matrix *T)
413 /* This is a hack. The EDomain_floor elements are possibly shared
414 * by many EDomain structures and we want to perform the substitution
415 * only once.
417 if (substituted)
418 return;
419 substituted = true;
421 evalue_substitute(e, subs);
423 assert(T->NbRows == v->Size-1);
424 Vector *tmp = Vector_Alloc(1+T->NbColumns);
425 Vector_Matrix_Product(v->p+1, T, tmp->p+1);
426 value_multiply(tmp->p[0], v->p[0], T->p[T->NbRows-1][T->NbColumns-1]);
427 Vector_Free(v);
428 v = tmp;
431 /* T is a homogeneous matrix that maps the original variables to the new variables.
432 * this has constraints in the new variables and this method
433 * transforms this to constraints in the original variables.
435 void EDomain::substitute(evalue **subs, Matrix *T, Matrix *Eq, unsigned MaxRays)
437 int nexist = floors.size();
438 Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
439 Polyhedron *new_D = DomainPreimage(D, M, MaxRays);
440 Polyhedron_Free(D);
441 D = new_D;
442 Matrix_Free(M);
444 M = nexist ? InsertColumns(Eq, 1+dimension(), nexist) : Eq;
445 new_D = DomainAddConstraints(D, M, MaxRays);
446 Polyhedron_Free(D);
447 D = new_D;
448 if (nexist)
449 Matrix_Free(M);
450 for (int i = 0; i < floors.size(); ++i)
451 floors[i]->substitute(subs, T);
454 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
456 Matrix M;
457 Polyhedron_Matrix_View(*P, &M, (*P)->NbEq);
459 Matrix *T = compress_variables(&M, nparam);
461 if (!T) {
462 *P = NULL;
463 return NULL;
465 if (isIdentity(T)) {
466 Matrix_Free(T);
467 T = NULL;
468 } else
469 *P = Polyhedron_Preimage(*P, T, MaxRays);
471 return T;
474 bool EDomain::not_empty(lexmin_options *options)
476 Polyhedron *P = D;
477 Polyhedron *Porig = P;
479 POL_ENSURE_VERTICES(P);
480 if (emptyQ2(P))
481 return false;
483 for (int i = 0; i < P->NbRays; ++i)
484 if (value_one_p(P->Ray[i][1+P->Dimension])) {
485 sample = Vector_Alloc(P->Dimension + 1);
486 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
487 return true;
490 if (options->emptiness_check == BV_LEXMIN_EMPTINESS_CHECK_COUNT) {
491 bool notzero;
492 Value cb;
493 value_init(cb);
494 barvinok_count_with_options(P, &cb, options->verify->barvinok);
495 notzero = value_notzero_p(cb);
496 value_clear(cb);
497 return notzero;
500 Matrix *T = NULL;
501 while (P && !emptyQ2(P) && P->NbEq > 0) {
502 Polyhedron *Q = P;
503 Matrix *T2 = remove_equalities(&P, 0, options->verify->barvinok->MaxRays);
504 if (!T)
505 T = T2;
506 else {
507 if (T2) {
508 Matrix *T3 = Matrix_Alloc(T->NbRows, T2->NbColumns);
509 Matrix_Product(T, T2, T3);
510 Matrix_Free(T);
511 Matrix_Free(T2);
512 T = T3;
514 if (Q != Porig)
515 Polyhedron_Free(Q);
518 if (P)
519 sample = Polyhedron_Sample(P, options->verify->barvinok);
520 if (sample) {
521 if (T) {
522 Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
523 Matrix_Vector_Product(T, sample->p, P_sample->p);
524 Vector_Free(sample);
525 sample = P_sample;
528 if (T) {
529 Polyhedron_Free(P);
530 Matrix_Free(T);
533 if (sample)
534 assert(in_domain(Porig, sample->p));
536 return sample != NULL;