add isl.py to distribution
[barvinok.git] / edomain.cc
blobc10af44c50ce635470227cac3754acfd34ead4a3
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 int i;
313 if (e->x.p->type != fractional)
314 continue;
315 if (find_floor(&e->x.p->arr[0]) == -1)
316 ++fract;
319 int rows = D->NbConstraints+2*fract+1;
320 int cols = 2+D->Dimension+fract;
321 ge->M = Matrix_Alloc(rows, cols);
322 for (int i = 0; i < D->NbConstraints; ++i) {
323 Vector_Copy(D->Constraint[i], ge->M->p[i], 1+D->Dimension);
324 value_assign(ge->M->p[i][1+D->Dimension+fract],
325 D->Constraint[i][1+D->Dimension]);
327 value_set_si(ge->M->p[rows-1][0], 1);
328 fract = 0;
329 evalue *e;
330 for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
331 if (e->x.p->type == fractional) {
332 int i, pos;
334 i = find_floor(&e->x.p->arr[0]);
335 if (i >= 0)
336 pos = D->Dimension-floors.size()+i;
337 else
338 pos = D->Dimension+fract;
340 add_fract(e, ge->M->p[rows-1], cols, 1+pos);
342 if (pos < D->Dimension)
343 continue;
345 EDomain_floor *new_floor;
346 new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
348 /* constraints for the new floor */
349 int row = D->NbConstraints+2*fract;
350 Vector_Copy(new_floor->v->p+1, ge->M->p[row]+1, dimension());
351 value_assign(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]);
352 value_oppose(ge->M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
353 value_set_si(ge->M->p[row][0], 1);
354 assert(value_eq(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]));
355 assert(Vector_Equal(new_floor->v->p+1, ge->M->p[row]+1, dimension()));
357 Vector_Scale(ge->M->p[row]+1, ge->M->p[row+1]+1, mone.x.n, cols-1);
358 value_set_si(ge->M->p[row+1][0], 1);
359 value_addto(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1],
360 ge->M->p[row+1][1+D->Dimension+fract]);
361 value_decrement(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1]);
363 ge->new_floors.push_back(new_floor);
365 ++fract;
366 } else {
367 assert(e->x.p->type == polynomial);
368 assert(e->x.p->size == 2);
369 add_coeff(ge->M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
372 add_coeff(ge->M->p[rows-1], cols, e, cols-1);
373 ge->simplified = ConstraintSimplify(ge->M->p[rows-1], ge->M->p[rows-1],
374 cols, &ge->M->p[rows-1][0]);
375 value_set_si(ge->M->p[rows-1][0], 1);
376 free_evalue_refs(&mone);
377 return ge;
380 /* "align" matrix to have nrows by inserting
381 * the necessary number of rows and an equal number of columns at the end
382 * right before the constant row/column
384 static Matrix *align_matrix_initial(Matrix *M, int nrows)
386 int i;
387 int newrows = nrows - M->NbRows;
388 Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
389 for (i = 0; i < M->NbRows-1; ++i) {
390 Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
391 value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
393 for (i = 0; i <= newrows; ++i)
394 value_assign(M2->p[M->NbRows-1+i][M->NbColumns-1+i],
395 M->p[M->NbRows-1][M->NbColumns-1]);
396 return M2;
399 static Matrix *InsertColumns(Matrix *M, int pos, int n)
401 Matrix *R;
402 int i;
404 R = Matrix_Alloc(M->NbRows, M->NbColumns+n);
405 for (i = 0; i < M->NbRows; ++i) {
406 Vector_Copy(M->p[i], R->p[i], pos);
407 Vector_Copy(M->p[i]+pos, R->p[i]+pos+n, M->NbColumns-pos);
409 return R;
412 void EDomain_floor::substitute(evalue **subs, Matrix *T)
414 /* This is a hack. The EDomain_floor elements are possibly shared
415 * by many EDomain structures and we want to perform the substitution
416 * only once.
418 if (substituted)
419 return;
420 substituted = true;
422 evalue_substitute(e, subs);
424 assert(T->NbRows == v->Size-1);
425 Vector *tmp = Vector_Alloc(1+T->NbColumns);
426 Vector_Matrix_Product(v->p+1, T, tmp->p+1);
427 value_multiply(tmp->p[0], v->p[0], T->p[T->NbRows-1][T->NbColumns-1]);
428 Vector_Free(v);
429 v = tmp;
432 /* T is a homogeneous matrix that maps the original variables to the new variables.
433 * this has constraints in the new variables and this method
434 * transforms this to constraints in the original variables.
436 void EDomain::substitute(evalue **subs, Matrix *T, Matrix *Eq, unsigned MaxRays)
438 int nexist = floors.size();
439 Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
440 Polyhedron *new_D = DomainPreimage(D, M, MaxRays);
441 Polyhedron_Free(D);
442 D = new_D;
443 Matrix_Free(M);
445 M = nexist ? InsertColumns(Eq, 1+dimension(), nexist) : Eq;
446 new_D = DomainAddConstraints(D, M, MaxRays);
447 Polyhedron_Free(D);
448 D = new_D;
449 if (nexist)
450 Matrix_Free(M);
451 for (int i = 0; i < floors.size(); ++i)
452 floors[i]->substitute(subs, T);
455 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
457 Matrix M;
458 Polyhedron_Matrix_View(*P, &M, (*P)->NbEq);
460 Matrix *T = compress_variables(&M, nparam);
462 if (!T) {
463 *P = NULL;
464 return NULL;
466 if (isIdentity(T)) {
467 Matrix_Free(T);
468 T = NULL;
469 } else
470 *P = Polyhedron_Preimage(*P, T, MaxRays);
472 return T;
475 bool EDomain::not_empty(lexmin_options *options)
477 Polyhedron *P = D;
478 Polyhedron *Porig = P;
480 POL_ENSURE_VERTICES(P);
481 if (emptyQ2(P))
482 return false;
484 for (int i = 0; i < P->NbRays; ++i)
485 if (value_one_p(P->Ray[i][1+P->Dimension])) {
486 sample = Vector_Alloc(P->Dimension + 1);
487 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
488 return true;
491 if (options->emptiness_check == BV_LEXMIN_EMPTINESS_CHECK_COUNT) {
492 bool notzero;
493 Value cb;
494 value_init(cb);
495 barvinok_count_with_options(P, &cb, options->verify->barvinok);
496 notzero = value_notzero_p(cb);
497 value_clear(cb);
498 return notzero;
501 Matrix *T = NULL;
502 while (P && !emptyQ2(P) && P->NbEq > 0) {
503 Polyhedron *Q = P;
504 Matrix *T2 = remove_equalities(&P, 0, options->verify->barvinok->MaxRays);
505 if (!T)
506 T = T2;
507 else {
508 if (T2) {
509 Matrix *T3 = Matrix_Alloc(T->NbRows, T2->NbColumns);
510 Matrix_Product(T, T2, T3);
511 Matrix_Free(T);
512 Matrix_Free(T2);
513 T = T3;
515 if (Q != Porig)
516 Polyhedron_Free(Q);
519 if (P)
520 sample = Polyhedron_Sample(P, options->verify->barvinok);
521 if (sample) {
522 if (T) {
523 Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
524 Matrix_Vector_Product(T, sample->p, P_sample->p);
525 Vector_Free(sample);
526 sample = P_sample;
529 if (T) {
530 Polyhedron_Free(P);
531 Matrix_Free(T);
534 if (sample)
535 assert(in_domain(Porig, sample->p));
537 return sample != NULL;