evalue.c: evalue_denom: declare variable at start of code block.
[barvinok.git] / lexmin.cc
blobc3e68dc1f37df3e8faef4a7c76a120dda8823e43
1 #include <iostream>
2 #include <vector>
3 #include <map>
4 #include <set>
5 #include <gmp.h>
6 #include <NTL/vec_ZZ.h>
7 #include <NTL/mat_ZZ.h>
8 extern "C" {
9 #include <polylib/polylibgmp.h>
11 #include <barvinok/barvinok.h>
12 #include <barvinok/evalue.h>
13 #include <barvinok/options.h>
14 #include <barvinok/util.h>
15 #include "conversion.h"
16 #include "decomposer.h"
17 #include "lattice_point.h"
18 #include "reduce_domain.h"
19 #include "mat_util.h"
20 #include "combine.h"
21 #include "edomain.h"
22 #include "evalue_util.h"
23 #include "sample.h"
24 #include "config.h"
26 #ifdef NTL_STD_CXX
27 using namespace NTL;
28 #endif
30 using std::vector;
31 using std::map;
32 using std::cerr;
33 using std::cout;
34 using std::endl;
35 using std::ostream;
37 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
38 #define RANGE 50
40 /* SRANGE : small range for evalutations */
41 #define SRANGE 15
43 /* if dimension >= BIDDIM, use SRANGE */
44 #define BIGDIM 5
46 /* VSRANGE : very small range for evalutations */
47 #define VSRANGE 5
49 /* if dimension >= VBIDDIM, use VSRANGE */
50 #define VBIGDIM 8
52 #ifndef HAVE_GETOPT_H
53 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
54 #else
55 #include <getopt.h>
56 #define NO_EMPTINESS_CHECK 256
57 struct option lexmin_options[] = {
58 { "verify", no_argument, 0, 'T' },
59 { "print-all", no_argument, 0, 'A' },
60 { "no-emptiness-check", no_argument, 0, NO_EMPTINESS_CHECK },
61 { "min", required_argument, 0, 'm' },
62 { "max", required_argument, 0, 'M' },
63 { "range", required_argument, 0, 'r' },
64 { "version", no_argument, 0, 'V' },
65 { 0, 0, 0, 0 }
67 #endif
69 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
71 static int type_offset(enode *p)
73 return p->type == fractional ? 1 :
74 p->type == flooring ? 1 : 0;
77 struct indicator_term {
78 int sign;
79 int pos;
80 mat_ZZ den;
81 evalue **vertex;
83 indicator_term(unsigned dim, int pos) {
84 den.SetDims(0, dim);
85 vertex = new evalue* [dim];
86 this->pos = pos;
87 sign = 0;
89 indicator_term(unsigned dim) {
90 den.SetDims(dim, dim);
91 vertex = new evalue* [dim];
92 pos = -1;
94 indicator_term(const indicator_term& src) {
95 sign = src.sign;
96 pos = src.pos;
97 den = src.den;
98 unsigned dim = den.NumCols();
99 vertex = new evalue* [dim];
100 for (int i = 0; i < dim; ++i) {
101 vertex[i] = new evalue();
102 value_init(vertex[i]->d);
103 evalue_copy(vertex[i], src.vertex[i]);
106 ~indicator_term() {
107 unsigned dim = den.NumCols();
108 for (int i = 0; i < dim; ++i) {
109 free_evalue_refs(vertex[i]);
110 delete vertex[i];
112 delete [] vertex;
114 void print(ostream& os, char **p) const;
115 void substitute(Matrix *T);
116 void normalize();
117 void substitute(evalue *fract, evalue *val);
118 void substitute(int pos, evalue *val);
119 void reduce_in_domain(Polyhedron *D);
120 bool is_opposite(indicator_term *neg);
123 bool indicator_term::is_opposite(indicator_term *neg)
125 if (sign + neg->sign != 0)
126 return false;
127 if (den != neg->den)
128 return false;
129 for (int k = 0; k < den.NumCols(); ++k)
130 if (!eequal(vertex[k], neg->vertex[k]))
131 return false;
132 return true;
135 void indicator_term::reduce_in_domain(Polyhedron *D)
137 for (int k = 0; k < den.NumCols(); ++k) {
138 reduce_evalue_in_domain(vertex[k], D);
139 if (evalue_range_reduction_in_domain(vertex[k], D))
140 reduce_evalue(vertex[k]);
144 void indicator_term::print(ostream& os, char **p) const
146 unsigned dim = den.NumCols();
147 unsigned factors = den.NumRows();
148 if (sign == 0)
149 os << " s ";
150 else if (sign > 0)
151 os << " + ";
152 else
153 os << " - ";
154 os << '[';
155 for (int i = 0; i < dim; ++i) {
156 if (i)
157 os << ',';
158 evalue_print(os, vertex[i], p);
160 os << ']';
161 for (int i = 0; i < factors; ++i) {
162 os << " + t" << i << "*[";
163 for (int j = 0; j < dim; ++j) {
164 if (j)
165 os << ',';
166 os << den[i][j];
168 os << "]";
170 if (sign == 0)
171 os << " (" << pos << ")";
174 /* Perform the substitution specified by T on the variables.
175 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
176 * The substitution is performed as in gen_fun::substitute
178 void indicator_term::substitute(Matrix *T)
180 unsigned dim = den.NumCols();
181 unsigned nparam = T->NbColumns - dim - 1;
182 unsigned newdim = T->NbRows - nparam - 1;
183 evalue **newvertex;
184 mat_ZZ trans;
185 matrix2zz(T, trans, newdim, dim);
186 trans = transpose(trans);
187 den *= trans;
188 newvertex = new evalue* [newdim];
190 vec_ZZ v;
191 v.SetLength(nparam+1);
193 evalue factor;
194 value_init(factor.d);
195 value_set_si(factor.d, 1);
196 value_init(factor.x.n);
197 for (int i = 0; i < newdim; ++i) {
198 values2zz(T->p[i]+dim, v, nparam+1);
199 newvertex[i] = multi_monom(v);
201 for (int j = 0; j < dim; ++j) {
202 if (value_zero_p(T->p[i][j]))
203 continue;
204 evalue term;
205 value_init(term.d);
206 evalue_copy(&term, vertex[j]);
207 value_assign(factor.x.n, T->p[i][j]);
208 emul(&factor, &term);
209 eadd(&term, newvertex[i]);
210 free_evalue_refs(&term);
213 free_evalue_refs(&factor);
214 for (int i = 0; i < dim; ++i) {
215 free_evalue_refs(vertex[i]);
216 delete vertex[i];
218 delete [] vertex;
219 vertex = newvertex;
222 static void evalue_add_constant(evalue *e, ZZ v)
224 Value tmp;
225 value_init(tmp);
227 /* go down to constant term */
228 while (value_zero_p(e->d))
229 e = &e->x.p->arr[type_offset(e->x.p)];
230 /* and add v */
231 zz2value(v, tmp);
232 value_multiply(tmp, tmp, e->d);
233 value_addto(e->x.n, e->x.n, tmp);
235 value_clear(tmp);
238 /* Make all powers in denominator lexico-positive */
239 void indicator_term::normalize()
241 vec_ZZ extra_vertex;
242 extra_vertex.SetLength(den.NumCols());
243 for (int r = 0; r < den.NumRows(); ++r) {
244 for (int k = 0; k < den.NumCols(); ++k) {
245 if (den[r][k] == 0)
246 continue;
247 if (den[r][k] > 0)
248 break;
249 sign = -sign;
250 den[r] = -den[r];
251 extra_vertex += den[r];
252 break;
255 for (int k = 0; k < extra_vertex.length(); ++k)
256 if (extra_vertex[k] != 0)
257 evalue_add_constant(vertex[k], extra_vertex[k]);
260 static void substitute(evalue *e, evalue *fract, evalue *val)
262 evalue *t;
264 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
265 if (t->x.p->type == fractional && eequal(&t->x.p->arr[0], fract))
266 break;
268 if (value_notzero_p(t->d))
269 return;
271 free_evalue_refs(&t->x.p->arr[0]);
272 evalue *term = &t->x.p->arr[2];
273 enode *p = t->x.p;
274 value_clear(t->d);
275 *t = t->x.p->arr[1];
277 emul(val, term);
278 eadd(term, e);
279 free_evalue_refs(term);
280 free(p);
282 reduce_evalue(e);
285 void indicator_term::substitute(evalue *fract, evalue *val)
287 unsigned dim = den.NumCols();
288 for (int i = 0; i < dim; ++i) {
289 ::substitute(vertex[i], fract, val);
293 static void substitute(evalue *e, int pos, evalue *val)
295 evalue *t;
297 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
298 if (t->x.p->type == polynomial && t->x.p->pos == pos)
299 break;
301 if (value_notzero_p(t->d))
302 return;
304 evalue *term = &t->x.p->arr[1];
305 enode *p = t->x.p;
306 value_clear(t->d);
307 *t = t->x.p->arr[0];
309 emul(val, term);
310 eadd(term, e);
311 free_evalue_refs(term);
312 free(p);
314 reduce_evalue(e);
317 void indicator_term::substitute(int pos, evalue *val)
319 unsigned dim = den.NumCols();
320 for (int i = 0; i < dim; ++i) {
321 ::substitute(vertex[i], pos, val);
325 struct indicator_constructor : public polar_decomposer, public vertex_decomposer {
326 vec_ZZ vertex;
327 vector<indicator_term*> *terms;
328 Matrix *T; /* Transformation to original space */
329 Param_Polyhedron *PP;
331 indicator_constructor(Polyhedron *P, unsigned dim, Param_Polyhedron *PP,
332 Matrix *T) :
333 vertex_decomposer(P, PP->nbV, this), T(T), PP(PP) {
334 vertex.SetLength(dim);
335 terms = new vector<indicator_term*>[nbV];
337 ~indicator_constructor() {
338 for (int i = 0; i < nbV; ++i)
339 for (int j = 0; j < terms[i].size(); ++j)
340 delete terms[i][j];
341 delete [] terms;
343 void substitute(Matrix *T);
344 void normalize();
345 void print(ostream& os, char **p);
347 virtual void handle_polar(Polyhedron *P, int sign);
350 void indicator_constructor::handle_polar(Polyhedron *C, int s)
352 unsigned dim = vertex.length();
354 assert(C->NbRays-1 == dim);
356 indicator_term *term = new indicator_term(dim);
357 term->sign = s;
358 terms[vert].push_back(term);
360 lattice_point(V, C, vertex, term->vertex);
362 for (int r = 0; r < dim; ++r) {
363 values2zz(C->Ray[r]+1, term->den[r], dim);
364 for (int j = 0; j < dim; ++j) {
365 if (term->den[r][j] == 0)
366 continue;
367 if (term->den[r][j] > 0)
368 break;
369 term->sign = -term->sign;
370 term->den[r] = -term->den[r];
371 vertex += term->den[r];
372 break;
376 for (int i = 0; i < dim; ++i) {
377 if (!term->vertex[i]) {
378 term->vertex[i] = new evalue();
379 value_init(term->vertex[i]->d);
380 value_init(term->vertex[i]->x.n);
381 zz2value(vertex[i], term->vertex[i]->x.n);
382 value_set_si(term->vertex[i]->d, 1);
383 continue;
385 if (vertex[i] == 0)
386 continue;
387 evalue_add_constant(term->vertex[i], vertex[i]);
390 if (T) {
391 term->substitute(T);
392 term->normalize();
395 lex_order_rows(term->den);
398 void indicator_constructor::print(ostream& os, char **p)
400 for (int i = 0; i < nbV; ++i)
401 for (int j = 0; j < terms[i].size(); ++j) {
402 os << "i: " << i << ", j: " << j << endl;
403 terms[i][j]->print(os, p);
404 os << endl;
408 void indicator_constructor::normalize()
410 for (int i = 0; i < nbV; ++i)
411 for (int j = 0; j < terms[i].size(); ++j) {
412 vec_ZZ vertex;
413 vertex.SetLength(terms[i][j]->den.NumCols());
414 for (int r = 0; r < terms[i][j]->den.NumRows(); ++r) {
415 for (int k = 0; k < terms[i][j]->den.NumCols(); ++k) {
416 if (terms[i][j]->den[r][k] == 0)
417 continue;
418 if (terms[i][j]->den[r][k] > 0)
419 break;
420 terms[i][j]->sign = -terms[i][j]->sign;
421 terms[i][j]->den[r] = -terms[i][j]->den[r];
422 vertex += terms[i][j]->den[r];
423 break;
426 lex_order_rows(terms[i][j]->den);
427 for (int k = 0; k < vertex.length(); ++k)
428 if (vertex[k] != 0)
429 evalue_add_constant(terms[i][j]->vertex[k], vertex[k]);
433 struct indicator;
435 enum order_sign { order_lt, order_le, order_eq, order_ge, order_gt, order_unknown };
437 struct partial_order {
438 indicator *ind;
440 map<indicator_term *, int > pred;
441 map<indicator_term *, vector<indicator_term * > > lt;
442 map<indicator_term *, vector<indicator_term * > > le;
443 map<indicator_term *, vector<indicator_term * > > eq;
445 map<indicator_term *, vector<indicator_term * > > pending;
447 partial_order(indicator *ind) : ind(ind) {}
448 void copy(const partial_order& order,
449 map< indicator_term *, indicator_term * > old2new);
451 order_sign compare(indicator_term *a, indicator_term *b);
452 void set_equal(indicator_term *a, indicator_term *b);
453 void unset_le(indicator_term *a, indicator_term *b);
455 bool compared(indicator_term* a, indicator_term* b);
456 void add(indicator_term* it, std::set<indicator_term *> *filter);
457 void remove(indicator_term* it);
459 void print(ostream& os, char **p);
462 void partial_order::unset_le(indicator_term *a, indicator_term *b)
464 vector<indicator_term *>::iterator i;
465 i = find(le[a].begin(), le[a].end(), b);
466 le[a].erase(i);
467 pred[b]--;
468 i = find(pending[a].begin(), pending[a].end(), b);
469 if (i != pending[a].end())
470 pending[a].erase(i);
473 void partial_order::set_equal(indicator_term *a, indicator_term *b)
475 if (eq[a].size() == 0)
476 eq[a].push_back(a);
477 if (eq[b].size() == 0)
478 eq[b].push_back(b);
479 a = eq[a][0];
480 b = eq[b][0];
481 assert(a != b);
482 if (b < a) {
483 indicator_term *c = a;
484 a = b;
485 b = c;
488 indicator_term *base = a;
490 map<indicator_term *, vector<indicator_term * > >::iterator i;
492 for (int j = 0; j < eq[b].size(); ++j) {
493 eq[base].push_back(eq[b][j]);
494 eq[eq[b][j]][0] = base;
496 eq[b].resize(1);
498 i = lt.find(b);
499 if (i != lt.end()) {
500 for (int j = 0; j < lt[b].size(); ++j) {
501 if (find(eq[base].begin(), eq[base].end(), lt[b][j]) != eq[base].end())
502 pred[lt[b][j]]--;
503 else if (find(lt[base].begin(), lt[base].end(), lt[b][j])
504 != lt[base].end())
505 pred[lt[b][j]]--;
506 else
507 lt[base].push_back(lt[b][j]);
509 lt.erase(i);
512 i = le.find(b);
513 if (i != le.end()) {
514 for (int j = 0; j < le[b].size(); ++j) {
515 if (find(eq[base].begin(), eq[base].end(), le[b][j]) != eq[base].end())
516 pred[le[b][j]]--;
517 else if (find(le[base].begin(), le[base].end(), le[b][j])
518 != le[base].end())
519 pred[le[b][j]]--;
520 else
521 le[base].push_back(le[b][j]);
523 le.erase(i);
526 i = pending.find(base);
527 if (i != pending.end()) {
528 vector<indicator_term * > old = pending[base];
529 pending[base].clear();
530 for (int j = 0; j < old.size(); ++j) {
531 if (find(eq[base].begin(), eq[base].end(), old[j]) == eq[base].end())
532 pending[base].push_back(old[j]);
536 i = pending.find(b);
537 if (i != pending.end()) {
538 for (int j = 0; j < pending[b].size(); ++j) {
539 if (find(eq[base].begin(), eq[base].end(), pending[b][j]) == eq[base].end())
540 pending[base].push_back(pending[b][j]);
542 pending.erase(i);
546 void partial_order::copy(const partial_order& order,
547 map< indicator_term *, indicator_term * > old2new)
549 map<indicator_term *, vector<indicator_term * > >::const_iterator i;
550 map<indicator_term *, int >::const_iterator j;
552 for (j = order.pred.begin(); j != order.pred.end(); ++j)
553 pred[old2new[(*j).first]] = (*j).second;
555 for (i = order.lt.begin(); i != order.lt.end(); ++i) {
556 for (int j = 0; j < (*i).second.size(); ++j)
557 lt[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
559 for (i = order.le.begin(); i != order.le.end(); ++i) {
560 for (int j = 0; j < (*i).second.size(); ++j)
561 le[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
563 for (i = order.eq.begin(); i != order.eq.end(); ++i) {
564 for (int j = 0; j < (*i).second.size(); ++j)
565 eq[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
567 for (i = order.pending.begin(); i != order.pending.end(); ++i) {
568 for (int j = 0; j < (*i).second.size(); ++j)
569 pending[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
573 static void interval_minmax(Polyhedron *I, int *min, int *max)
575 assert(I->Dimension == 1);
576 *min = 1;
577 *max = -1;
578 POL_ENSURE_VERTICES(I);
579 for (int i = 0; i < I->NbRays; ++i) {
580 if (value_pos_p(I->Ray[i][1]))
581 *max = 1;
582 else if (value_neg_p(I->Ray[i][1]))
583 *min = -1;
584 else {
585 if (*max < 0)
586 *max = 0;
587 if (*min > 0)
588 *min = 0;
593 static void interval_minmax(Polyhedron *D, Matrix *T, int *min, int *max,
594 unsigned MaxRays)
596 Polyhedron *I = Polyhedron_Image(D, T, MaxRays);
597 if (MaxRays & POL_INTEGER)
598 I = DomainConstraintSimplify(I, MaxRays);
599 if (emptyQ2(I)) {
600 Polyhedron_Free(I);
601 I = Polyhedron_Image(D, T, MaxRays);
603 interval_minmax(I, min, max);
604 Polyhedron_Free(I);
607 struct max_term {
608 unsigned dim;
609 EDomain *domain;
610 vector<evalue *> max;
612 void print(ostream& os, char **p, barvinok_options *options) const;
613 void substitute(Matrix *T, unsigned MaxRays);
614 Vector *eval(Value *val, unsigned MaxRays) const;
616 ~max_term() {
617 for (int i = 0; i < max.size(); ++i) {
618 free_evalue_refs(max[i]);
619 delete max[i];
621 delete domain;
626 * Project on first dim dimensions
628 Polyhedron* Polyhedron_Project_Initial(Polyhedron *P, int dim)
630 int i;
631 Matrix *T;
632 Polyhedron *I;
634 if (P->Dimension == dim)
635 return Polyhedron_Copy(P);
637 T = Matrix_Alloc(dim+1, P->Dimension+1);
638 for (i = 0; i < dim; ++i)
639 value_set_si(T->p[i][i], 1);
640 value_set_si(T->p[dim][P->Dimension], 1);
641 I = Polyhedron_Image(P, T, P->NbConstraints);
642 Matrix_Free(T);
643 return I;
646 struct indicator {
647 vector<indicator_term*> term;
648 indicator_constructor& ic;
649 partial_order order;
650 EDomain *D;
651 Polyhedron *P;
652 Param_Domain *PD;
653 barvinok_options *options;
655 indicator(indicator_constructor& ic, Param_Domain *PD, EDomain *D,
656 barvinok_options *options) :
657 ic(ic), PD(PD), D(D), order(this), options(options), P(NULL) {}
658 indicator(const indicator& ind, EDomain *D) :
659 ic(ind.ic), PD(ind.PD), D(NULL), order(this), options(ind.options),
660 P(Polyhedron_Copy(ind.P)) {
661 map< indicator_term *, indicator_term * > old2new;
662 for (int i = 0; i < ind.term.size(); ++i) {
663 indicator_term *it = new indicator_term(*ind.term[i]);
664 old2new[ind.term[i]] = it;
665 term.push_back(it);
667 order.copy(ind.order, old2new);
668 set_domain(D);
670 ~indicator() {
671 for (int i = 0; i < term.size(); ++i)
672 delete term[i];
673 if (D)
674 delete D;
675 if (P)
676 Polyhedron_Free(P);
679 void set_domain(EDomain *D) {
680 if (this->D)
681 delete this->D;
682 this->D = D;
683 int nparam = ic.P->Dimension - ic.vertex.length();
684 Polyhedron *Q = Polyhedron_Project_Initial(D->D, nparam);
685 Q = DomainConstraintSimplify(Q, options->MaxRays);
686 if (!P || !PolyhedronIncludes(Q, P))
687 reduce_in_domain(Q);
688 if (P)
689 Polyhedron_Free(P);
690 P = Q;
693 void add(const indicator_term* it);
694 void remove(indicator_term* it);
695 void remove_initial_rational_vertices();
696 void expand_rational_vertex(indicator_term *initial);
698 void print(ostream& os, char **p);
699 void simplify();
700 void peel(int i, int j);
701 void combine(indicator_term *a, indicator_term *b);
702 void substitute(evalue *equation);
703 void reduce_in_domain(Polyhedron *D);
704 bool handle_equal_numerators(indicator_term *base);
706 max_term* create_max_term(indicator_term *it);
709 max_term* indicator::create_max_term(indicator_term *it)
711 int dim = it->den.NumCols();
712 int nparam = ic.P->Dimension - ic.vertex.length();
713 max_term *maximum = new max_term;
714 maximum->dim = nparam;
715 maximum->domain = new EDomain(D);
716 for (int j = 0; j < dim; ++j) {
717 evalue *E = new evalue;
718 value_init(E->d);
719 evalue_copy(E, it->vertex[j]);
720 if (evalue_frac2floor_in_domain(E, D->D))
721 reduce_evalue(E);
722 maximum->max.push_back(E);
724 return maximum;
727 /* compute a - b */
728 static evalue *ediff(const evalue *a, const evalue *b)
730 evalue mone;
731 value_init(mone.d);
732 evalue_set_si(&mone, -1, 1);
733 evalue *diff = new evalue;
734 value_init(diff->d);
735 evalue_copy(diff, b);
736 emul(&mone, diff);
737 eadd(a, diff);
738 reduce_evalue(diff);
739 free_evalue_refs(&mone);
740 return diff;
743 static order_sign evalue_sign(evalue *diff, EDomain *D, unsigned MaxRays)
745 order_sign sign = order_eq;
746 evalue mone;
747 value_init(mone.d);
748 evalue_set_si(&mone, -1, 1);
749 int len = 1 + D->D->Dimension + 1;
750 Vector *c = Vector_Alloc(len);
751 Matrix *T = Matrix_Alloc(2, len-1);
753 int fract = evalue2constraint(D, diff, c->p, len);
754 Vector_Copy(c->p+1, T->p[0], len-1);
755 value_assign(T->p[1][len-2], c->p[0]);
757 int min, max;
758 interval_minmax(D->D, T, &min, &max, MaxRays);
759 if (max < 0)
760 sign = order_lt;
761 else {
762 if (fract) {
763 emul(&mone, diff);
764 evalue2constraint(D, diff, c->p, len);
765 emul(&mone, diff);
766 Vector_Copy(c->p+1, T->p[0], len-1);
767 value_assign(T->p[1][len-2], c->p[0]);
769 int negmin, negmax;
770 interval_minmax(D->D, T, &negmin, &negmax, MaxRays);
771 min = -negmax;
773 if (min > 0)
774 sign = order_gt;
775 else if (max == 0 && min == 0)
776 sign = order_eq;
777 else if (min < 0 && max > 0)
778 sign = order_unknown;
779 else if (min < 0)
780 sign = order_le;
781 else
782 sign = order_ge;
785 Matrix_Free(T);
786 Vector_Free(c);
787 free_evalue_refs(&mone);
789 return sign;
792 order_sign partial_order::compare(indicator_term *a, indicator_term *b)
794 unsigned dim = a->den.NumCols();
795 order_sign sign = order_eq;
796 EDomain *D = ind->D;
797 unsigned MaxRays = ind->options->MaxRays;
798 if (MaxRays & POL_INTEGER && (a->sign == 0 || b->sign == 0))
799 MaxRays = 0;
801 for (int k = 0; k < dim; ++k) {
802 /* compute a->vertex[k] - b->vertex[k] */
803 evalue *diff = ediff(a->vertex[k], b->vertex[k]);
804 order_sign diff_sign = evalue_sign(diff, D, MaxRays);
806 if (diff_sign == order_lt) {
807 if (sign == order_eq || sign == order_le)
808 sign = order_lt;
809 else
810 sign = order_unknown;
811 free_evalue_refs(diff);
812 delete diff;
813 break;
815 if (diff_sign == order_gt) {
816 if (sign == order_eq || sign == order_ge)
817 sign = order_gt;
818 else
819 sign = order_unknown;
820 free_evalue_refs(diff);
821 delete diff;
822 break;
824 if (diff_sign == order_eq) {
825 if (D == ind->D && !EVALUE_IS_ZERO(*diff))
826 ind->substitute(diff);
827 free_evalue_refs(diff);
828 delete diff;
829 continue;
831 if ((diff_sign == order_unknown) ||
832 ((diff_sign == order_lt || diff_sign == order_le) && sign == order_ge) ||
833 ((diff_sign == order_gt || diff_sign == order_ge) && sign == order_le)) {
834 free_evalue_refs(diff);
835 delete diff;
836 sign = order_unknown;
837 break;
840 sign = diff_sign;
842 Matrix *M;
843 vector<EDomain_floor *> new_floors;
844 M = D->add_ge_constraint(diff, new_floors);
845 value_set_si(M->p[M->NbRows-1][0], 0);
846 Polyhedron *D2 = Constraints2Polyhedron(M, MaxRays);
847 EDomain *EDeq = new EDomain(D2, D, new_floors);
848 Polyhedron_Free(D2);
849 Matrix_Free(M);
850 for (int i = 0; i < new_floors.size(); ++i)
851 EDomain_floor::unref(new_floors[i]);
853 if (D != ind->D)
854 delete D;
855 D = EDeq;
857 free_evalue_refs(diff);
858 delete diff;
861 if (D != ind->D)
862 delete D;
864 return sign;
867 bool partial_order::compared(indicator_term* a, indicator_term* b)
869 map<indicator_term *, vector<indicator_term * > >::iterator j;
871 j = lt.find(a);
872 if (j != lt.end() && find(lt[a].begin(), lt[a].end(), b) != lt[a].end())
873 return true;
875 j = le.find(a);
876 if (j != le.end() && find(le[a].begin(), le[a].end(), b) != le[a].end())
877 return true;
879 return false;
882 void partial_order::add(indicator_term* it, std::set<indicator_term *> *filter)
884 if (eq.find(it) != eq.end() && eq[it].size() == 1)
885 return;
887 int it_pred = filter ? pred[it] : 0;
889 map<indicator_term *, int >::iterator i;
890 for (i = pred.begin(); i != pred.end(); ++i) {
891 if ((*i).second != 0)
892 continue;
893 if (eq.find((*i).first) != eq.end() && eq[(*i).first].size() == 1)
894 continue;
895 if (filter) {
896 if ((*i).first == it)
897 continue;
898 if (filter->find((*i).first) == filter->end())
899 continue;
900 if (compared((*i).first, it))
901 continue;
903 order_sign sign = compare(it, (*i).first);
904 if (sign == order_lt) {
905 lt[it].push_back((*i).first);
906 (*i).second++;
907 } else if (sign == order_le) {
908 le[it].push_back((*i).first);
909 (*i).second++;
910 } else if (sign == order_eq) {
911 pred[it] = it_pred;
912 set_equal(it, (*i).first);
913 return;
914 } else if (sign == order_gt) {
915 pending[(*i).first].push_back(it);
916 lt[(*i).first].push_back(it);
917 ++it_pred;
918 } else if (sign == order_ge) {
919 pending[(*i).first].push_back(it);
920 le[(*i).first].push_back(it);
921 ++it_pred;
924 pred[it] = it_pred;
927 void partial_order::remove(indicator_term* it)
929 std::set<indicator_term *> filter;
930 map<indicator_term *, vector<indicator_term * > >::iterator i;
932 assert(pred[it] == 0);
934 i = eq.find(it);
935 if (i != eq.end()) {
936 assert(eq[it].size() >= 1);
937 indicator_term *base;
938 if (eq[it].size() == 1) {
939 base = eq[it][0];
940 eq.erase(i);
942 vector<indicator_term * >::iterator j;
943 j = find(eq[base].begin(), eq[base].end(), it);
944 assert(j != eq[base].end());
945 eq[base].erase(j);
946 } else {
947 /* "it" may no longer be the smallest, since the order
948 * structure may have been copied from another one.
950 sort(eq[it].begin()+1, eq[it].end());
951 assert(eq[it][0] == it);
952 eq[it].erase(eq[it].begin());
953 base = eq[it][0];
954 eq[base] = eq[it];
955 eq.erase(i);
957 for (int j = 1; j < eq[base].size(); ++j)
958 eq[eq[base][j]][0] = base;
960 i = lt.find(it);
961 if (i != lt.end()) {
962 lt[base] = lt[it];
963 lt.erase(i);
966 i = le.find(it);
967 if (i != le.end()) {
968 le[base] = le[it];
969 le.erase(i);
972 i = pending.find(it);
973 if (i != pending.end()) {
974 pending[base] = pending[it];
975 pending.erase(i);
979 if (eq[base].size() == 1)
980 eq.erase(base);
982 map<indicator_term *, int >::iterator j;
983 j = pred.find(it);
984 pred.erase(j);
986 return;
989 i = lt.find(it);
990 if (i != lt.end()) {
991 for (int j = 0; j < lt[it].size(); ++j) {
992 filter.insert(lt[it][j]);
993 pred[lt[it][j]]--;
995 lt.erase(i);
998 i = le.find(it);
999 if (i != le.end()) {
1000 for (int j = 0; j < le[it].size(); ++j) {
1001 filter.insert(le[it][j]);
1002 pred[le[it][j]]--;
1004 le.erase(i);
1007 map<indicator_term *, int >::iterator j;
1008 j = pred.find(it);
1009 pred.erase(j);
1011 i = pending.find(it);
1012 if (i != pending.end()) {
1013 for (int j = 0; j < pending[it].size(); ++j) {
1014 filter.erase(pending[it][j]);
1015 add(pending[it][j], &filter);
1017 pending.erase(i);
1021 void partial_order::print(ostream& os, char **p)
1023 map<indicator_term *, vector<indicator_term * > >::iterator i;
1024 for (i = lt.begin(); i != lt.end(); ++i) {
1025 (*i).first->print(os, p);
1026 assert(pred.find((*i).first) != pred.end());
1027 os << "(" << pred[(*i).first] << ")";
1028 os << " < ";
1029 for (int j = 0; j < (*i).second.size(); ++j) {
1030 if (j)
1031 os << ", ";
1032 (*i).second[j]->print(os, p);
1033 assert(pred.find((*i).second[j]) != pred.end());
1034 os << "(" << pred[(*i).second[j]] << ")";
1036 os << endl;
1038 for (i = le.begin(); i != le.end(); ++i) {
1039 (*i).first->print(os, p);
1040 assert(pred.find((*i).first) != pred.end());
1041 os << "(" << pred[(*i).first] << ")";
1042 os << " <= ";
1043 for (int j = 0; j < (*i).second.size(); ++j) {
1044 if (j)
1045 os << ", ";
1046 (*i).second[j]->print(os, p);
1047 assert(pred.find((*i).second[j]) != pred.end());
1048 os << "(" << pred[(*i).second[j]] << ")";
1050 os << endl;
1052 for (i = eq.begin(); i != eq.end(); ++i) {
1053 if ((*i).second.size() <= 1)
1054 continue;
1055 (*i).first->print(os, p);
1056 assert(pred.find((*i).first) != pred.end());
1057 os << "(" << pred[(*i).first] << ")";
1058 for (int j = 1; j < (*i).second.size(); ++j) {
1059 if (j)
1060 os << " = ";
1061 (*i).second[j]->print(os, p);
1062 assert(pred.find((*i).second[j]) != pred.end());
1063 os << "(" << pred[(*i).second[j]] << ")";
1065 os << endl;
1067 for (i = pending.begin(); i != pending.end(); ++i) {
1068 os << "pending on ";
1069 (*i).first->print(os, p);
1070 assert(pred.find((*i).first) != pred.end());
1071 os << "(" << pred[(*i).first] << ")";
1072 os << ": ";
1073 for (int j = 0; j < (*i).second.size(); ++j) {
1074 if (j)
1075 os << ", ";
1076 (*i).second[j]->print(os, p);
1077 assert(pred.find((*i).second[j]) != pred.end());
1078 os << "(" << pred[(*i).second[j]] << ")";
1080 os << endl;
1084 void indicator::add(const indicator_term* it)
1086 indicator_term *nt = new indicator_term(*it);
1087 nt->reduce_in_domain(P ? P : D->D);
1088 term.push_back(nt);
1089 order.add(nt, NULL);
1090 assert(term.size() == order.pred.size());
1093 void indicator::remove(indicator_term* it)
1095 vector<indicator_term *>::iterator i;
1096 i = find(term.begin(), term.end(), it);
1097 assert(i!= term.end());
1098 order.remove(it);
1099 term.erase(i);
1100 assert(term.size() == order.pred.size());
1101 delete it;
1104 void indicator::expand_rational_vertex(indicator_term *initial)
1106 int pos = initial->pos;
1107 remove(initial);
1108 if (ic.terms[pos].size() == 0) {
1109 Param_Vertices *V;
1110 FORALL_PVertex_in_ParamPolyhedron(V, PD, ic.PP) // _i is internal counter
1111 if (_i == pos) {
1112 ic.decompose_at_vertex(V, pos, options->MaxRays);
1113 break;
1115 END_FORALL_PVertex_in_ParamPolyhedron;
1117 for (int j = 0; j < ic.terms[pos].size(); ++j)
1118 add(ic.terms[pos][j]);
1121 void indicator::remove_initial_rational_vertices()
1123 do {
1124 indicator_term *initial = NULL;
1125 map<indicator_term *, int >::iterator i;
1126 for (i = order.pred.begin(); i != order.pred.end(); ++i) {
1127 if ((*i).second != 0)
1128 continue;
1129 if ((*i).first->sign != 0)
1130 continue;
1131 if (order.eq.find((*i).first) != order.eq.end() &&
1132 order.eq[(*i).first].size() <= 1)
1133 continue;
1134 initial = (*i).first;
1135 break;
1137 if (!initial)
1138 break;
1139 expand_rational_vertex(initial);
1140 } while(1);
1143 void indicator::reduce_in_domain(Polyhedron *D)
1145 for (int i = 0; i < term.size(); ++i)
1146 term[i]->reduce_in_domain(D);
1149 void indicator::print(ostream& os, char **p)
1151 assert(term.size() == order.pred.size());
1152 for (int i = 0; i < term.size(); ++i) {
1153 term[i]->print(os, p);
1154 os << endl;
1156 order.print(os, p);
1159 /* Remove pairs of opposite terms */
1160 void indicator::simplify()
1162 for (int i = 0; i < term.size(); ++i) {
1163 for (int j = i+1; j < term.size(); ++j) {
1164 if (term[i]->sign + term[j]->sign != 0)
1165 continue;
1166 if (term[i]->den != term[j]->den)
1167 continue;
1168 int k;
1169 for (k = 0; k < term[i]->den.NumCols(); ++k)
1170 if (!eequal(term[i]->vertex[k], term[j]->vertex[k]))
1171 break;
1172 if (k < term[i]->den.NumCols())
1173 continue;
1174 delete term[i];
1175 delete term[j];
1176 term.erase(term.begin()+j);
1177 term.erase(term.begin()+i);
1178 --i;
1179 break;
1184 void indicator::peel(int i, int j)
1186 if (j < i) {
1187 int tmp = i;
1188 i = j;
1189 j = tmp;
1191 assert (i < j);
1192 int dim = term[i]->den.NumCols();
1194 mat_ZZ common;
1195 mat_ZZ rest_i;
1196 mat_ZZ rest_j;
1197 int n_common = 0, n_i = 0, n_j = 0;
1199 common.SetDims(min(term[i]->den.NumRows(), term[j]->den.NumRows()), dim);
1200 rest_i.SetDims(term[i]->den.NumRows(), dim);
1201 rest_j.SetDims(term[j]->den.NumRows(), dim);
1203 int k, l;
1204 for (k = 0, l = 0; k < term[i]->den.NumRows() && l < term[j]->den.NumRows(); ) {
1205 int s = lex_cmp(term[i]->den[k], term[j]->den[l]);
1206 if (s == 0) {
1207 common[n_common++] = term[i]->den[k];
1208 ++k;
1209 ++l;
1210 } else if (s < 0)
1211 rest_i[n_i++] = term[i]->den[k++];
1212 else
1213 rest_j[n_j++] = term[j]->den[l++];
1215 while (k < term[i]->den.NumRows())
1216 rest_i[n_i++] = term[i]->den[k++];
1217 while (l < term[j]->den.NumRows())
1218 rest_j[n_j++] = term[j]->den[l++];
1219 common.SetDims(n_common, dim);
1220 rest_i.SetDims(n_i, dim);
1221 rest_j.SetDims(n_j, dim);
1223 for (k = 0; k <= n_i; ++k) {
1224 indicator_term *it = new indicator_term(*term[i]);
1225 it->den.SetDims(n_common + k, dim);
1226 for (l = 0; l < n_common; ++l)
1227 it->den[l] = common[l];
1228 for (l = 0; l < k; ++l)
1229 it->den[n_common+l] = rest_i[l];
1230 lex_order_rows(it->den);
1231 if (k)
1232 for (l = 0; l < dim; ++l)
1233 evalue_add_constant(it->vertex[l], rest_i[k-1][l]);
1234 term.push_back(it);
1237 for (k = 0; k <= n_j; ++k) {
1238 indicator_term *it = new indicator_term(*term[j]);
1239 it->den.SetDims(n_common + k, dim);
1240 for (l = 0; l < n_common; ++l)
1241 it->den[l] = common[l];
1242 for (l = 0; l < k; ++l)
1243 it->den[n_common+l] = rest_j[l];
1244 lex_order_rows(it->den);
1245 if (k)
1246 for (l = 0; l < dim; ++l)
1247 evalue_add_constant(it->vertex[l], rest_j[k-1][l]);
1248 term.push_back(it);
1250 term.erase(term.begin()+j);
1251 term.erase(term.begin()+i);
1254 void indicator::combine(indicator_term *a, indicator_term *b)
1256 int dim = a->den.NumCols();
1258 mat_ZZ common;
1259 mat_ZZ rest_i;
1260 mat_ZZ rest_j;
1261 int n_common = 0, n_i = 0, n_j = 0;
1263 common.SetDims(min(a->den.NumRows(), b->den.NumRows()), dim);
1264 rest_i.SetDims(a->den.NumRows(), dim);
1265 rest_j.SetDims(b->den.NumRows(), dim);
1267 int k, l;
1268 for (k = 0, l = 0; k < a->den.NumRows() && l < b->den.NumRows(); ) {
1269 int s = lex_cmp(a->den[k], b->den[l]);
1270 if (s == 0) {
1271 common[n_common++] = a->den[k];
1272 ++k;
1273 ++l;
1274 } else if (s < 0)
1275 rest_i[n_i++] = a->den[k++];
1276 else
1277 rest_j[n_j++] = b->den[l++];
1279 while (k < a->den.NumRows())
1280 rest_i[n_i++] = a->den[k++];
1281 while (l < b->den.NumRows())
1282 rest_j[n_j++] = b->den[l++];
1283 common.SetDims(n_common, dim);
1284 rest_i.SetDims(n_i, dim);
1285 rest_j.SetDims(n_j, dim);
1287 assert(n_i < 30);
1288 assert(n_j < 30);
1290 for (k = (1 << n_i)-1; k >= 0; --k) {
1291 indicator_term *it = k ? new indicator_term(*b) : b;
1292 it->den.SetDims(n_common + n_i + n_j, dim);
1293 for (l = 0; l < n_common; ++l)
1294 it->den[l] = common[l];
1295 for (l = 0; l < n_i; ++l)
1296 it->den[n_common+l] = rest_i[l];
1297 for (l = 0; l < n_j; ++l)
1298 it->den[n_common+n_i+l] = rest_j[l];
1299 lex_order_rows(it->den);
1300 int change = 0;
1301 for (l = 0; l < n_i; ++l) {
1302 if (!(k & (1 <<l)))
1303 continue;
1304 change ^= 1;
1305 for (int m = 0; m < dim; ++m)
1306 evalue_add_constant(it->vertex[m], rest_i[l][m]);
1308 if (change)
1309 it->sign = -it->sign;
1310 if (it != b) {
1311 term.push_back(it);
1312 order.add(it, NULL);
1316 for (k = (1 << n_j)-1; k >= 0; --k) {
1317 indicator_term *it = k ? new indicator_term(*a) : a;
1318 it->den.SetDims(n_common + n_i + n_j, dim);
1319 for (l = 0; l < n_common; ++l)
1320 it->den[l] = common[l];
1321 for (l = 0; l < n_i; ++l)
1322 it->den[n_common+l] = rest_i[l];
1323 for (l = 0; l < n_j; ++l)
1324 it->den[n_common+n_i+l] = rest_j[l];
1325 lex_order_rows(it->den);
1326 int change = 0;
1327 for (l = 0; l < n_j; ++l) {
1328 if (!(k & (1 <<l)))
1329 continue;
1330 change ^= 1;
1331 for (int m = 0; m < dim; ++m)
1332 evalue_add_constant(it->vertex[m], rest_j[l][m]);
1334 if (change)
1335 it->sign = -it->sign;
1336 if (it != a) {
1337 term.push_back(it);
1338 order.add(it, NULL);
1343 bool indicator::handle_equal_numerators(indicator_term *base)
1345 for (int i = 0; i < order.eq[base].size(); ++i) {
1346 for (int j = i+1; j < order.eq[base].size(); ++j) {
1347 if (order.eq[base][i]->is_opposite(order.eq[base][j])) {
1348 remove(order.eq[base][j]);
1349 remove(i ? order.eq[base][i] : base);
1350 return true;
1354 for (int j = 1; j < order.eq[base].size(); ++j)
1355 if (order.eq[base][j]->sign != base->sign) {
1356 combine(base, order.eq[base][j]);
1357 return true;
1359 return false;
1362 void indicator::substitute(evalue *equation)
1364 evalue *fract = NULL;
1365 evalue *val = new evalue;
1366 value_init(val->d);
1367 evalue_copy(val, equation);
1369 evalue factor;
1370 value_init(factor.d);
1371 value_init(factor.x.n);
1373 evalue *e;
1374 for (e = val; value_zero_p(e->d) && e->x.p->type != fractional;
1375 e = &e->x.p->arr[type_offset(e->x.p)])
1378 if (value_zero_p(e->d) && e->x.p->type == fractional)
1379 fract = &e->x.p->arr[0];
1380 else {
1381 for (e = val; value_zero_p(e->d) && e->x.p->type != polynomial;
1382 e = &e->x.p->arr[type_offset(e->x.p)])
1384 assert(value_zero_p(e->d) && e->x.p->type == polynomial);
1387 int offset = type_offset(e->x.p);
1389 assert(value_notzero_p(e->x.p->arr[offset+1].d));
1390 assert(value_notzero_p(e->x.p->arr[offset+1].x.n));
1391 if (value_neg_p(e->x.p->arr[offset+1].x.n)) {
1392 value_oppose(factor.d, e->x.p->arr[offset+1].x.n);
1393 value_assign(factor.x.n, e->x.p->arr[offset+1].d);
1394 } else {
1395 value_assign(factor.d, e->x.p->arr[offset+1].x.n);
1396 value_oppose(factor.x.n, e->x.p->arr[offset+1].d);
1399 free_evalue_refs(&e->x.p->arr[offset+1]);
1400 enode *p = e->x.p;
1401 value_clear(e->d);
1402 *e = e->x.p->arr[offset];
1404 emul(&factor, val);
1406 if (fract) {
1407 for (int i = 0; i < term.size(); ++i)
1408 term[i]->substitute(fract, val);
1410 free_evalue_refs(&p->arr[0]);
1411 } else {
1412 for (int i = 0; i < term.size(); ++i)
1413 term[i]->substitute(p->pos, val);
1416 free_evalue_refs(&factor);
1417 free_evalue_refs(val);
1418 delete val;
1420 free(p);
1423 static void print_varlist(ostream& os, int n, char **names)
1425 int i;
1426 os << "[";
1427 for (i = 0; i < n; ++i) {
1428 if (i)
1429 os << ",";
1430 os << names[i];
1432 os << "]";
1435 void max_term::print(ostream& os, char **p, barvinok_options *options) const
1437 os << "{ ";
1438 print_varlist(os, dim, p);
1439 os << " -> ";
1440 os << "[";
1441 for (int i = 0; i < max.size(); ++i) {
1442 if (i)
1443 os << ",";
1444 evalue_print(os, max[i], p);
1446 os << "]";
1447 os << " : ";
1448 domain->print_constraints(os, p, options);
1449 os << " }" << endl;
1452 static void evalue_substitute(evalue *e, evalue **subs)
1454 evalue *v;
1456 if (value_notzero_p(e->d))
1457 return;
1459 enode *p = e->x.p;
1460 for (int i = 0; i < p->size; ++i)
1461 evalue_substitute(&p->arr[i], subs);
1463 if (p->type == polynomial)
1464 v = subs[p->pos-1];
1465 else {
1466 v = new evalue;
1467 value_init(v->d);
1468 value_set_si(v->d, 0);
1469 v->x.p = new_enode(p->type, 3, -1);
1470 value_clear(v->x.p->arr[0].d);
1471 v->x.p->arr[0] = p->arr[0];
1472 evalue_set_si(&v->x.p->arr[1], 0, 1);
1473 evalue_set_si(&v->x.p->arr[2], 1, 1);
1476 int offset = type_offset(p);
1478 for (int i = p->size-1; i >= offset+1; i--) {
1479 emul(v, &p->arr[i]);
1480 eadd(&p->arr[i], &p->arr[i-1]);
1481 free_evalue_refs(&(p->arr[i]));
1484 if (p->type != polynomial) {
1485 free_evalue_refs(v);
1486 delete v;
1489 value_clear(e->d);
1490 *e = p->arr[offset];
1491 free(p);
1494 /* "align" matrix to have nrows by inserting
1495 * the necessary number of rows and an equal number of columns at the end
1496 * right before the constant row/column
1498 static Matrix *align_matrix_initial(Matrix *M, int nrows)
1500 int i;
1501 int newrows = nrows - M->NbRows;
1502 Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
1503 for (i = 0; i < newrows; ++i)
1504 value_set_si(M2->p[M->NbRows-1+i][M->NbColumns-1+i], 1);
1505 for (i = 0; i < M->NbRows-1; ++i) {
1506 Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
1507 value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
1509 value_assign(M2->p[M2->NbRows-1][M2->NbColumns-1],
1510 M->p[M->NbRows-1][M->NbColumns-1]);
1511 return M2;
1514 /* T maps the compressed parameters to the original parameters,
1515 * while this max_term is based on the compressed parameters
1516 * and we want get the original parameters back.
1518 void max_term::substitute(Matrix *T, unsigned MaxRays)
1520 int nexist = domain->D->Dimension - (T->NbColumns-1);
1521 Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
1523 Polyhedron *D = DomainImage(domain->D, M, MaxRays);
1524 Polyhedron_Free(domain->D);
1525 domain->D = D;
1526 Matrix_Free(M);
1528 assert(T->NbRows == T->NbColumns);
1529 Matrix *T2 = Matrix_Copy(T);
1530 Matrix *inv = Matrix_Alloc(T->NbColumns, T->NbRows);
1531 int ok = Matrix_Inverse(T2, inv);
1532 Matrix_Free(T2);
1533 assert(ok);
1535 evalue denom;
1536 value_init(denom.d);
1537 value_init(denom.x.n);
1538 value_set_si(denom.x.n, 1);
1539 value_assign(denom.d, inv->p[inv->NbRows-1][inv->NbColumns-1]);
1541 vec_ZZ v;
1542 v.SetLength(inv->NbColumns);
1543 evalue* subs[inv->NbRows-1];
1544 for (int i = 0; i < inv->NbRows-1; ++i) {
1545 values2zz(inv->p[i], v, v.length());
1546 subs[i] = multi_monom(v);
1547 emul(&denom, subs[i]);
1549 free_evalue_refs(&denom);
1551 for (int i = 0; i < max.size(); ++i) {
1552 evalue_substitute(max[i], subs);
1553 reduce_evalue(max[i]);
1556 for (int i = 0; i < inv->NbRows-1; ++i) {
1557 free_evalue_refs(subs[i]);
1558 delete subs[i];
1560 Matrix_Free(inv);
1563 int Last_Non_Zero(Value *p, unsigned len)
1565 for (int i = len-1; i >= 0; --i)
1566 if (value_notzero_p(p[i]))
1567 return i;
1568 return -1;
1571 static void SwapColumns(Value **V, int n, int i, int j)
1573 for (int r = 0; r < n; ++r)
1574 value_swap(V[r][i], V[r][j]);
1577 static void SwapColumns(Polyhedron *P, int i, int j)
1579 SwapColumns(P->Constraint, P->NbConstraints, i, j);
1580 SwapColumns(P->Ray, P->NbRays, i, j);
1583 void compute_evalue(evalue *e, Value *val, Value *res)
1585 double d = compute_evalue(e, val);
1586 if (d > 0)
1587 d += .25;
1588 else
1589 d -= .25;
1590 value_set_double(*res, d);
1593 Vector *max_term::eval(Value *val, unsigned MaxRays) const
1595 if (!domain->contains(val, dim))
1596 return NULL;
1597 Vector *res = Vector_Alloc(max.size());
1598 for (int i = 0; i < max.size(); ++i) {
1599 compute_evalue(max[i], val, &res->p[i]);
1601 return res;
1604 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays);
1606 Vector *Polyhedron_not_empty(Polyhedron *P, unsigned MaxRays)
1608 Polyhedron *Porig = P;
1609 Vector *sample = NULL;
1611 POL_ENSURE_VERTICES(P);
1612 if (emptyQ2(P))
1613 return NULL;
1615 for (int i = 0; i < P->NbRays; ++i)
1616 if (value_one_p(P->Ray[i][1+P->Dimension])) {
1617 sample = Vector_Alloc(P->Dimension + 1);
1618 Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
1619 return sample;
1622 Matrix *T = remove_equalities(&P, 0, MaxRays);
1623 if (P)
1624 sample = Polyhedron_Sample(P, MaxRays);
1625 if (sample) {
1626 if (T) {
1627 Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
1628 Matrix_Vector_Product(T, sample->p, P_sample->p);
1629 Vector_Free(sample);
1630 sample = P_sample;
1633 if (T) {
1634 Polyhedron_Free(P);
1635 Matrix_Free(T);
1638 return sample;
1641 struct split {
1642 evalue *constraint;
1643 enum sign { le, ge, lge } sign;
1645 split (evalue *c, enum sign s) : constraint(c), sign(s) {}
1648 static void split_on(const split& sp, EDomain *D,
1649 EDomain **Dlt, EDomain **Deq, EDomain **Dgt,
1650 barvinok_options *options)
1652 Matrix *M, *M2;
1653 EDomain *ED[3];
1654 Polyhedron *D2;
1655 Value mone;
1656 value_init(mone);
1657 value_set_si(mone, -1);
1658 *Dlt = NULL;
1659 *Deq = NULL;
1660 *Dgt = NULL;
1661 vector<EDomain_floor *> new_floors;
1662 M = D->add_ge_constraint(sp.constraint, new_floors);
1663 if (sp.sign == split::lge || sp.sign == split::ge) {
1664 M2 = Matrix_Copy(M);
1665 value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1666 M2->p[M2->NbRows-1][M2->NbColumns-1]);
1667 D2 = Constraints2Polyhedron(M2, options->MaxRays);
1668 ED[2] = new EDomain(D2, D, new_floors);
1669 Polyhedron_Free(D2);
1670 Matrix_Free(M2);
1671 } else
1672 ED[2] = NULL;
1673 if (sp.sign == split::lge || sp.sign == split::le) {
1674 M2 = Matrix_Copy(M);
1675 Vector_Scale(M2->p[M2->NbRows-1]+1, M2->p[M2->NbRows-1]+1,
1676 mone, M2->NbColumns-1);
1677 value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
1678 M2->p[M2->NbRows-1][M2->NbColumns-1]);
1679 D2 = Constraints2Polyhedron(M2, options->MaxRays);
1680 ED[0] = new EDomain(D2, D, new_floors);
1681 Polyhedron_Free(D2);
1682 Matrix_Free(M2);
1683 } else
1684 ED[0] = NULL;
1686 assert(sp.sign == split::lge || sp.sign == split::ge || sp.sign == split::le);
1687 value_set_si(M->p[M->NbRows-1][0], 0);
1688 D2 = Constraints2Polyhedron(M, options->MaxRays);
1689 ED[1] = new EDomain(D2, D, new_floors);
1690 Polyhedron_Free(D2);
1691 Matrix_Free(M);
1693 Vector *sample = D->sample;
1694 if (sample && new_floors.size() > 0) {
1695 assert(sample->Size == D->D->Dimension+1);
1696 sample = Vector_Alloc(D->D->Dimension+new_floors.size()+1);
1697 Vector_Copy(D->sample->p, sample->p, D->D->Dimension);
1698 value_set_si(sample->p[D->D->Dimension+new_floors.size()], 1);
1699 for (int i = 0; i < new_floors.size(); ++i)
1700 new_floors[i]->eval(sample->p, &sample->p[D->D->Dimension+i]);
1703 for (int i = 0; i < new_floors.size(); ++i)
1704 EDomain_floor::unref(new_floors[i]);
1706 for (int i = 0; i < 3; ++i) {
1707 if (!ED[i])
1708 continue;
1709 if (sample && ED[i]->contains(sample->p, sample->Size-1)) {
1710 ED[i]->sample = Vector_Alloc(sample->Size);
1711 Vector_Copy(sample->p, ED[i]->sample->p, sample->Size);
1712 } else if (emptyQ2(ED[i]->D) ||
1713 (options->emptiness_check == 1 &&
1714 !(ED[i]->sample = Polyhedron_not_empty(ED[i]->D,
1715 options->MaxRays)))) {
1716 delete ED[i];
1717 ED[i] = NULL;
1720 *Dlt = ED[0];
1721 *Deq = ED[1];
1722 *Dgt = ED[2];
1723 value_clear(mone);
1724 if (sample != D->sample)
1725 Vector_Free(sample);
1728 ostream & operator<< (ostream & os, const vector<int> & v)
1730 os << "[";
1731 for (int i = 0; i < v.size(); ++i) {
1732 if (i)
1733 os << ", ";
1734 os << v[i];
1736 os << "]";
1737 return os;
1740 static bool isTranslation(Matrix *M)
1742 unsigned i, j;
1743 if (M->NbRows != M->NbColumns)
1744 return False;
1746 for (i = 0;i < M->NbRows; i ++)
1747 for (j = 0; j < M->NbColumns-1; j ++)
1748 if (i == j) {
1749 if(value_notone_p(M->p[i][j]))
1750 return False;
1751 } else {
1752 if(value_notzero_p(M->p[i][j]))
1753 return False;
1755 return value_one_p(M->p[M->NbRows-1][M->NbColumns-1]);
1758 static Matrix *compress_parameters(Polyhedron **P, Polyhedron **C,
1759 unsigned nparam, unsigned MaxRays)
1761 Matrix *M, *T, *CP;
1763 /* compress_parms doesn't like equalities that only involve parameters */
1764 for (int i = 0; i < (*P)->NbEq; ++i)
1765 assert(First_Non_Zero((*P)->Constraint[i]+1, (*P)->Dimension-nparam) != -1);
1767 M = Matrix_Alloc((*P)->NbEq, (*P)->Dimension+2);
1768 Vector_Copy((*P)->Constraint[0], M->p[0], (*P)->NbEq * ((*P)->Dimension+2));
1769 CP = compress_parms(M, nparam);
1770 Matrix_Free(M);
1772 if (isTranslation(CP)) {
1773 Matrix_Free(CP);
1774 return NULL;
1777 T = align_matrix(CP, (*P)->Dimension+1);
1778 *P = Polyhedron_Preimage(*P, T, MaxRays);
1779 Matrix_Free(T);
1781 *C = Polyhedron_Preimage(*C, CP, MaxRays);
1783 return CP;
1786 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
1788 /* Matrix "view" of equalities */
1789 Matrix M;
1790 M.NbRows = (*P)->NbEq;
1791 M.NbColumns = (*P)->Dimension+2;
1792 M.p_Init = (*P)->p_Init;
1793 M.p = (*P)->Constraint;
1795 Matrix *T = compress_variables(&M, nparam);
1797 if (!T) {
1798 *P = NULL;
1799 return NULL;
1801 if (isIdentity(T)) {
1802 Matrix_Free(T);
1803 T = NULL;
1804 } else
1805 *P = Polyhedron_Preimage(*P, T, MaxRays);
1807 return T;
1810 void construct_rational_vertices(Param_Polyhedron *PP, Matrix *T, unsigned dim,
1811 int nparam, vector<indicator_term *>& vertices)
1813 int i;
1814 Param_Vertices *PV;
1815 Value lcm, tmp;
1816 value_init(lcm);
1817 value_init(tmp);
1819 vec_ZZ v;
1820 v.SetLength(nparam+1);
1822 evalue factor;
1823 value_init(factor.d);
1824 value_init(factor.x.n);
1825 value_set_si(factor.x.n, 1);
1826 value_set_si(factor.d, 1);
1828 for (i = 0, PV = PP->V; PV; ++i, PV = PV->next) {
1829 indicator_term *term = new indicator_term(dim, i);
1830 vertices.push_back(term);
1831 Matrix *M = Matrix_Alloc(PV->Vertex->NbRows+nparam+1, nparam+1);
1832 value_set_si(lcm, 1);
1833 for (int j = 0; j < PV->Vertex->NbRows; ++j)
1834 value_lcm(lcm, PV->Vertex->p[j][nparam+1], &lcm);
1835 value_assign(M->p[M->NbRows-1][M->NbColumns-1], lcm);
1836 for (int j = 0; j < PV->Vertex->NbRows; ++j) {
1837 value_division(tmp, lcm, PV->Vertex->p[j][nparam+1]);
1838 Vector_Scale(PV->Vertex->p[j], M->p[j], tmp, nparam+1);
1840 for (int j = 0; j < nparam; ++j)
1841 value_assign(M->p[PV->Vertex->NbRows+j][j], lcm);
1842 if (T) {
1843 Matrix *M2 = Matrix_Alloc(T->NbRows, M->NbColumns);
1844 Matrix_Product(T, M, M2);
1845 Matrix_Free(M);
1846 M = M2;
1848 for (int j = 0; j < dim; ++j) {
1849 values2zz(M->p[j], v, nparam+1);
1850 term->vertex[j] = multi_monom(v);
1851 value_assign(factor.d, lcm);
1852 emul(&factor, term->vertex[j]);
1854 Matrix_Free(M);
1856 assert(i == PP->nbV);
1857 free_evalue_refs(&factor);
1858 value_clear(lcm);
1859 value_clear(tmp);
1862 /* An auxiliary class that keeps a reference to an evalue
1863 * and frees it when it goes out of scope.
1865 struct temp_evalue {
1866 evalue *E;
1867 temp_evalue() : E(NULL) {}
1868 temp_evalue(evalue *e) : E(e) {}
1869 operator evalue* () const { return E; }
1870 evalue *operator=(evalue *e) {
1871 if (E) {
1872 free_evalue_refs(E);
1873 delete E;
1875 E = e;
1876 return E;
1878 ~temp_evalue() {
1879 if (E) {
1880 free_evalue_refs(E);
1881 delete E;
1886 static vector<max_term*> lexmin(indicator& ind, unsigned nparam,
1887 vector<int> loc)
1889 vector<max_term*> maxima;
1890 map<indicator_term *, int >::iterator i;
1891 vector<int> best_score;
1892 vector<int> second_score;
1893 vector<int> neg_score;
1895 do {
1896 indicator_term *best = NULL, *second = NULL, *neg = NULL,
1897 *neg_eq = NULL, *neg_le = NULL;
1898 for (i = ind.order.pred.begin(); i != ind.order.pred.end(); ++i) {
1899 vector<int> score;
1900 if ((*i).second != 0)
1901 continue;
1902 indicator_term *term = (*i).first;
1903 if (term->sign == 0) {
1904 ind.expand_rational_vertex(term);
1905 break;
1908 if (ind.order.eq.find(term) != ind.order.eq.end()) {
1909 int j;
1910 if (ind.order.eq[term].size() <= 1)
1911 continue;
1912 for (j = 1; j < ind.order.eq[term].size(); ++j)
1913 if (ind.order.pred[ind.order.eq[term][j]] != 0)
1914 break;
1915 if (j < ind.order.eq[term].size())
1916 continue;
1917 score.push_back(ind.order.eq[term].size());
1918 } else
1919 score.push_back(0);
1920 if (ind.order.le.find(term) != ind.order.le.end())
1921 score.push_back(ind.order.le[term].size());
1922 else
1923 score.push_back(0);
1924 if (ind.order.lt.find(term) != ind.order.lt.end())
1925 score.push_back(-ind.order.lt[term].size());
1926 else
1927 score.push_back(0);
1929 if (term->sign > 0) {
1930 if (!best || score < best_score) {
1931 second = best;
1932 second_score = best_score;
1933 best = term;
1934 best_score = score;
1935 } else if (!second || score < second_score) {
1936 second = term;
1937 second_score = score;
1939 } else {
1940 if (!neg_eq && ind.order.eq.find(term) != ind.order.eq.end()) {
1941 for (int j = 1; j < ind.order.eq[term].size(); ++j)
1942 if (ind.order.eq[term][j]->sign != term->sign) {
1943 neg_eq = term;
1944 break;
1947 if (!neg_le && ind.order.le.find(term) != ind.order.le.end())
1948 neg_le = term;
1949 if (!neg || score < neg_score) {
1950 neg = term;
1951 neg_score = score;
1955 if (i != ind.order.pred.end())
1956 continue;
1958 if (!best && neg_eq) {
1959 assert(ind.order.eq[neg_eq].size() != 0);
1960 bool handled = ind.handle_equal_numerators(neg_eq);
1961 assert(handled);
1962 continue;
1965 if (!best && neg_le) {
1966 /* The smallest term is negative and <= some positive term */
1967 best = neg_le;
1968 neg = NULL;
1971 if (!best) {
1972 /* apparently there can be negative initial term on empty domains */
1973 if (ind.options->emptiness_check == 1)
1974 assert(!neg);
1975 break;
1978 if (!second && !neg) {
1979 indicator_term *rat = NULL;
1980 assert(best);
1981 if (ind.order.le[best].size() == 0) {
1982 if (ind.order.eq[best].size() != 0) {
1983 bool handled = ind.handle_equal_numerators(best);
1984 if (ind.options->emptiness_check == 1)
1985 assert(handled);
1986 /* If !handled then the leading coefficient is bigger than one;
1987 * must be an empty domain
1989 if (handled)
1990 continue;
1991 else
1992 break;
1994 maxima.push_back(ind.create_max_term(best));
1995 break;
1997 for (int j = 0; j < ind.order.le[best].size(); ++j) {
1998 if (ind.order.le[best][j]->sign == 0) {
1999 if (!rat && ind.order.pred[ind.order.le[best][j]] == 1)
2000 rat = ind.order.le[best][j];
2001 } else if (ind.order.le[best][j]->sign > 0) {
2002 second = ind.order.le[best][j];
2003 break;
2004 } else if (!neg)
2005 neg = ind.order.le[best][j];
2008 if (!second && !neg) {
2009 assert(rat);
2010 ind.order.unset_le(best, rat);
2011 ind.expand_rational_vertex(rat);
2012 continue;
2015 if (!second)
2016 second = neg;
2018 ind.order.unset_le(best, second);
2021 if (!second)
2022 second = neg;
2024 unsigned dim = best->den.NumCols();
2025 temp_evalue diff;
2026 order_sign sign;
2027 for (int k = 0; k < dim; ++k) {
2028 diff = ediff(best->vertex[k], second->vertex[k]);
2029 sign = evalue_sign(diff, ind.D, ind.options->MaxRays);
2031 /* neg can never be smaller than best, unless it may still cancel.
2032 * This can happen if positive terms have been determined to
2033 * be equal or less than or equal to some negative term.
2035 if (second == neg && !neg_eq && !neg_le) {
2036 if (sign == order_ge)
2037 sign = order_eq;
2038 if (sign == order_unknown)
2039 sign = order_le;
2042 if (sign != order_eq)
2043 break;
2044 if (!EVALUE_IS_ZERO(*diff))
2045 ind.substitute(diff);
2047 if (sign == order_eq) {
2048 ind.order.set_equal(best, second);
2049 continue;
2051 if (sign == order_lt) {
2052 ind.order.lt[best].push_back(second);
2053 ind.order.pred[second]++;
2054 continue;
2056 if (sign == order_gt) {
2057 ind.order.lt[second].push_back(best);
2058 ind.order.pred[best]++;
2059 continue;
2062 split sp(diff, sign == order_le ? split::le :
2063 sign == order_ge ? split::ge : split::lge);
2065 EDomain *Dlt, *Deq, *Dgt;
2066 split_on(sp, ind.D, &Dlt, &Deq, &Dgt, ind.options);
2067 if (ind.options->emptiness_check == 1)
2068 assert(Dlt || Deq || Dgt);
2069 else if (!(Dlt || Deq || Dgt))
2070 /* Must have been empty all along */
2071 break;
2073 if (Deq && (Dlt || Dgt)) {
2074 int locsize = loc.size();
2075 loc.push_back(0);
2076 indicator indeq(ind, Deq);
2077 Deq = NULL;
2078 indeq.substitute(diff);
2079 vector<max_term*> maxeq = lexmin(indeq, nparam, loc);
2080 maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
2081 loc.resize(locsize);
2083 if (Dgt && Dlt) {
2084 int locsize = loc.size();
2085 loc.push_back(1);
2086 indicator indgt(ind, Dgt);
2087 Dgt = NULL;
2088 /* we don't know the new location of these terms in indgt */
2090 indgt.order.lt[second].push_back(best);
2091 indgt.order.pred[best]++;
2093 vector<max_term*> maxgt = lexmin(indgt, nparam, loc);
2094 maxima.insert(maxima.end(), maxgt.begin(), maxgt.end());
2095 loc.resize(locsize);
2098 if (Deq) {
2099 loc.push_back(0);
2100 ind.substitute(diff);
2101 ind.set_domain(Deq);
2103 if (Dlt) {
2104 loc.push_back(-1);
2105 ind.order.lt[best].push_back(second);
2106 ind.order.pred[second]++;
2107 ind.set_domain(Dlt);
2109 if (Dgt) {
2110 loc.push_back(1);
2111 ind.order.lt[second].push_back(best);
2112 ind.order.pred[best]++;
2113 ind.set_domain(Dgt);
2115 } while(1);
2117 return maxima;
2120 static vector<max_term*> lexmin(Polyhedron *P, Polyhedron *C,
2121 barvinok_options *options)
2123 unsigned nparam = C->Dimension;
2124 Param_Polyhedron *PP = NULL;
2125 Polyhedron *CEq = NULL, *pVD;
2126 Matrix *CT = NULL;
2127 Matrix *T = NULL, *CP = NULL;
2128 Param_Domain *D, *next;
2129 Param_Vertices *V;
2130 Polyhedron *Porig = P;
2131 Polyhedron *Corig = C;
2132 vector<max_term*> all_max;
2133 Polyhedron *Q;
2134 unsigned P2PSD_MaxRays;
2136 if (emptyQ2(P))
2137 return all_max;
2139 POL_ENSURE_VERTICES(P);
2141 if (emptyQ2(P))
2142 return all_max;
2144 assert(P->NbBid == 0);
2146 if (P->NbEq > 0) {
2147 if (nparam > 0)
2148 CP = compress_parameters(&P, &C, nparam, options->MaxRays);
2149 Q = P;
2150 T = remove_equalities(&P, nparam, options->MaxRays);
2151 if (P != Q && Q != Porig)
2152 Polyhedron_Free(Q);
2153 if (!P) {
2154 if (C != Corig)
2155 Polyhedron_Free(C);
2156 return all_max;
2160 if (options->MaxRays & POL_NO_DUAL)
2161 P2PSD_MaxRays = 0;
2162 else
2163 P2PSD_MaxRays = options->MaxRays;
2165 Q = P;
2166 PP = Polyhedron2Param_SimplifiedDomain(&P, C, P2PSD_MaxRays, &CEq, &CT);
2167 if (P != Q && Q != Porig)
2168 Polyhedron_Free(Q);
2170 if (CT) {
2171 if (isIdentity(CT)) {
2172 Matrix_Free(CT);
2173 CT = NULL;
2174 } else
2175 nparam = CT->NbRows - 1;
2177 assert(!CT);
2179 unsigned dim = P->Dimension - nparam;
2181 int nd;
2182 for (nd = 0, D=PP->D; D; ++nd, D=D->next);
2183 Polyhedron **fVD = new Polyhedron*[nd];
2185 indicator_constructor ic(P, dim, PP, T);
2187 vector<indicator_term *> all_vertices;
2188 construct_rational_vertices(PP, T, T ? T->NbRows-nparam-1 : dim,
2189 nparam, all_vertices);
2191 for (nd = 0, D=PP->D; D; D=next) {
2192 next = D->next;
2194 Polyhedron *rVD = reduce_domain(D->Domain, CT, CEq,
2195 fVD, nd, options->MaxRays);
2196 if (!rVD)
2197 continue;
2199 pVD = CT ? DomainImage(rVD,CT,options->MaxRays) : rVD;
2201 EDomain *epVD = new EDomain(pVD);
2202 indicator ind(ic, D, epVD, options);
2204 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
2205 ind.add(all_vertices[_i]);
2206 END_FORALL_PVertex_in_ParamPolyhedron;
2208 ind.remove_initial_rational_vertices();
2210 vector<int> loc;
2211 vector<max_term*> maxima = lexmin(ind, nparam, loc);
2212 if (CP)
2213 for (int j = 0; j < maxima.size(); ++j)
2214 maxima[j]->substitute(CP, options->MaxRays);
2215 all_max.insert(all_max.end(), maxima.begin(), maxima.end());
2217 ++nd;
2218 if (rVD != pVD)
2219 Domain_Free(pVD);
2220 Domain_Free(rVD);
2222 for (int i = 0; i < all_vertices.size(); ++i)
2223 delete all_vertices[i];
2224 if (CP)
2225 Matrix_Free(CP);
2226 if (T)
2227 Matrix_Free(T);
2228 Param_Polyhedron_Free(PP);
2229 if (CEq)
2230 Polyhedron_Free(CEq);
2231 for (--nd; nd >= 0; --nd) {
2232 Domain_Free(fVD[nd]);
2234 delete [] fVD;
2235 if (C != Corig)
2236 Polyhedron_Free(C);
2237 if (P != Porig)
2238 Polyhedron_Free(P);
2240 return all_max;
2243 static void verify_results(Polyhedron *A, Polyhedron *C,
2244 vector<max_term*>& maxima, int m, int M,
2245 int print_all, unsigned MaxRays);
2247 int main(int argc, char **argv)
2249 Polyhedron *A, *C;
2250 Matrix *MA;
2251 int bignum;
2252 char **iter_names, **param_names;
2253 int c, ind = 0;
2254 int range = 0;
2255 int verify = 0;
2256 int print_all = 0;
2257 int m = INT_MAX, M = INT_MIN, r;
2258 int print_solution = 1;
2259 struct barvinok_options *options;
2261 options = barvinok_options_new_with_defaults();
2263 while ((c = getopt_long(argc, argv, "TAm:M:r:V", lexmin_options, &ind)) != -1) {
2264 switch (c) {
2265 case NO_EMPTINESS_CHECK:
2266 options->emptiness_check = 0;
2267 break;
2268 case 'T':
2269 verify = 1;
2270 break;
2271 case 'A':
2272 print_all = 1;
2273 break;
2274 case 'm':
2275 m = atoi(optarg);
2276 verify = 1;
2277 break;
2278 case 'M':
2279 M = atoi(optarg);
2280 verify = 1;
2281 break;
2282 case 'r':
2283 M = atoi(optarg);
2284 m = -M;
2285 verify = 1;
2286 break;
2287 case 'V':
2288 printf(barvinok_version());
2289 exit(0);
2290 break;
2294 MA = Matrix_Read();
2295 C = Constraints2Polyhedron(MA, options->MaxRays);
2296 Matrix_Free(MA);
2297 fscanf(stdin, " %d", &bignum);
2298 assert(bignum == -1);
2299 MA = Matrix_Read();
2300 A = Constraints2Polyhedron(MA, options->MaxRays);
2301 Matrix_Free(MA);
2303 if (A->Dimension >= VBIGDIM)
2304 r = VSRANGE;
2305 else if (A->Dimension >= BIGDIM)
2306 r = SRANGE;
2307 else
2308 r = RANGE;
2309 if (M == INT_MIN)
2310 M = r;
2311 if (m == INT_MAX)
2312 m = -r;
2314 if (verify && m > M) {
2315 fprintf(stderr,"Nothing to do: min > max !\n");
2316 return(0);
2318 if (verify)
2319 print_solution = 0;
2321 iter_names = util_generate_names(A->Dimension - C->Dimension, "i");
2322 param_names = util_generate_names(C->Dimension, "p");
2323 if (print_solution) {
2324 Polyhedron_Print(stdout, P_VALUE_FMT, A);
2325 Polyhedron_Print(stdout, P_VALUE_FMT, C);
2327 vector<max_term*> maxima = lexmin(A, C, options);
2328 if (print_solution)
2329 for (int i = 0; i < maxima.size(); ++i)
2330 maxima[i]->print(cout, param_names, options);
2332 if (verify)
2333 verify_results(A, C, maxima, m, M, print_all, options->MaxRays);
2335 for (int i = 0; i < maxima.size(); ++i)
2336 delete maxima[i];
2338 util_free_names(A->Dimension - C->Dimension, iter_names);
2339 util_free_names(C->Dimension, param_names);
2340 Polyhedron_Free(A);
2341 Polyhedron_Free(C);
2343 free(options);
2345 return 0;
2348 static bool lexmin(int pos, Polyhedron *P, Value *context)
2350 Value LB, UB, k;
2351 int lu_flags;
2352 bool found = false;
2354 if (emptyQ(P))
2355 return false;
2357 value_init(LB); value_init(UB); value_init(k);
2358 value_set_si(LB,0);
2359 value_set_si(UB,0);
2360 lu_flags = lower_upper_bounds(pos,P,context,&LB,&UB);
2361 assert(!(lu_flags & LB_INFINITY));
2363 value_set_si(context[pos],0);
2364 if (!lu_flags && value_lt(UB,LB)) {
2365 value_clear(LB); value_clear(UB); value_clear(k);
2366 return false;
2368 if (!P->next) {
2369 value_assign(context[pos], LB);
2370 value_clear(LB); value_clear(UB); value_clear(k);
2371 return true;
2373 for (value_assign(k,LB); lu_flags || value_le(k,UB); value_increment(k,k)) {
2374 value_assign(context[pos],k);
2375 if ((found = lexmin(pos+1, P->next, context)))
2376 break;
2378 if (!found)
2379 value_set_si(context[pos],0);
2380 value_clear(LB); value_clear(UB); value_clear(k);
2381 return found;
2384 static void print_list(FILE *out, Value *z, char* brackets, int len)
2386 fprintf(out, "%c", brackets[0]);
2387 value_print(out, VALUE_FMT,z[0]);
2388 for (int k = 1; k < len; ++k) {
2389 fprintf(out, ", ");
2390 value_print(out, VALUE_FMT,z[k]);
2392 fprintf(out, "%c", brackets[1]);
2395 static int check_poly(Polyhedron *S, Polyhedron *CS, vector<max_term*>& maxima,
2396 int nparam, int pos, Value *z, int print_all, int st,
2397 unsigned MaxRays)
2399 if (pos == nparam) {
2400 int k;
2401 bool found = lexmin(1, S, z);
2403 if (print_all) {
2404 printf("lexmin");
2405 print_list(stdout, z+S->Dimension-nparam+1, "()", nparam);
2406 printf(" = ");
2407 if (found)
2408 print_list(stdout, z+1, "[]", S->Dimension-nparam);
2409 printf(" ");
2412 Vector *min = NULL;
2413 for (int i = 0; i < maxima.size(); ++i)
2414 if ((min = maxima[i]->eval(z+S->Dimension-nparam+1, MaxRays)))
2415 break;
2417 int ok = !(found ^ !!min);
2418 if (found && min)
2419 for (int i = 0; i < S->Dimension-nparam; ++i)
2420 if (value_ne(z[1+i], min->p[i])) {
2421 ok = 0;
2422 break;
2424 if (!ok) {
2425 printf("\n");
2426 fflush(stdout);
2427 fprintf(stderr, "Error !\n");
2428 fprintf(stderr, "lexmin");
2429 print_list(stderr, z+S->Dimension-nparam+1, "()", nparam);
2430 fprintf(stderr, " should be ");
2431 if (found)
2432 print_list(stderr, z+1, "[]", S->Dimension-nparam);
2433 fprintf(stderr, " while digging gives ");
2434 if (min)
2435 print_list(stderr, min->p, "[]", S->Dimension-nparam);
2436 fprintf(stderr, ".\n");
2437 return 0;
2438 } else if (print_all)
2439 printf("OK.\n");
2440 if (min)
2441 Vector_Free(min);
2443 for (k = 1; k <= S->Dimension-nparam; ++k)
2444 value_set_si(z[k], 0);
2445 } else {
2446 Value tmp;
2447 Value LB, UB;
2448 value_init(tmp);
2449 value_init(LB);
2450 value_init(UB);
2451 int ok =
2452 !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
2453 for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
2454 if (!print_all) {
2455 int k = VALUE_TO_INT(tmp);
2456 if (!pos && !(k%st)) {
2457 printf("o");
2458 fflush(stdout);
2461 value_assign(z[pos+S->Dimension-nparam+1],tmp);
2462 if (!check_poly(S, CS->next, maxima, nparam, pos+1, z, print_all, st,
2463 MaxRays)) {
2464 value_clear(tmp);
2465 value_clear(LB);
2466 value_clear(UB);
2467 return 0;
2470 value_set_si(z[pos+S->Dimension-nparam+1],0);
2471 value_clear(tmp);
2472 value_clear(LB);
2473 value_clear(UB);
2475 return 1;
2478 void verify_results(Polyhedron *A, Polyhedron *C, vector<max_term*>& maxima,
2479 int m, int M, int print_all, unsigned MaxRays)
2481 Polyhedron *CC, *CC2, *CS, *S;
2482 unsigned nparam = C->Dimension;
2483 Value *p;
2484 int i;
2485 int st;
2487 CC = Polyhedron_Project(A, nparam);
2488 CC2 = DomainIntersection(C, CC, MaxRays);
2489 Domain_Free(CC);
2490 CC = CC2;
2492 /* Intersect context with range */
2493 if (nparam > 0) {
2494 Matrix *MM;
2495 Polyhedron *U;
2497 MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
2498 for (int i = 0; i < C->Dimension; ++i) {
2499 value_set_si(MM->p[2*i][0], 1);
2500 value_set_si(MM->p[2*i][1+i], 1);
2501 value_set_si(MM->p[2*i][1+C->Dimension], -m);
2502 value_set_si(MM->p[2*i+1][0], 1);
2503 value_set_si(MM->p[2*i+1][1+i], -1);
2504 value_set_si(MM->p[2*i+1][1+C->Dimension], M);
2506 CC2 = AddConstraints(MM->p[0], 2*CC->Dimension, CC, MaxRays);
2507 U = Universe_Polyhedron(0);
2508 CS = Polyhedron_Scan(CC2, U, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2509 Polyhedron_Free(U);
2510 Polyhedron_Free(CC2);
2511 Matrix_Free(MM);
2512 } else
2513 CS = NULL;
2515 p = ALLOCN(Value, A->Dimension+2);
2516 for (i=0; i <= A->Dimension; i++) {
2517 value_init(p[i]);
2518 value_set_si(p[i],0);
2520 value_init(p[i]);
2521 value_set_si(p[i], 1);
2523 S = Polyhedron_Scan(A, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2525 if (!print_all && C->Dimension > 0) {
2526 if (M-m > 80)
2527 st = 1+(M-m)/80;
2528 else
2529 st = 1;
2530 for (int i = m; i <= M; i += st)
2531 printf(".");
2532 printf( "\r" );
2533 fflush(stdout);
2536 if (S) {
2537 if (!(CS && emptyQ2(CS)))
2538 check_poly(S, CS, maxima, nparam, 0, p, print_all, st, MaxRays);
2539 Domain_Free(S);
2542 if (!print_all)
2543 printf("\n");
2545 for (i=0; i <= (A->Dimension+1); i++)
2546 value_clear(p[i]);
2547 free(p);
2548 if (CS)
2549 Domain_Free(CS);
2550 Polyhedron_Free(CC);