keep track of number of Bernoulli sums
[barvinok.git] / lexmin.cc
blob0ac3f08c14bed5a7c2080048d5b0e7e8b1b6d861
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <map>
5 #include <set>
6 #define partition STL_PARTITION
7 #include <algorithm>
8 #undef partition
9 #include <gmp.h>
10 #include <NTL/vec_ZZ.h>
11 #include <NTL/mat_ZZ.h>
12 #include <barvinok/barvinok.h>
13 #include <barvinok/evalue.h>
14 #include <barvinok/options.h>
15 #include <barvinok/util.h>
16 #include "argp.h"
17 #include "progname.h"
18 #include "conversion.h"
19 #include "decomposer.h"
20 #include "lattice_point.h"
21 #include "reduce_domain.h"
22 #include "mat_util.h"
23 #include "combine.h"
24 #include "edomain.h"
25 #include "evalue_util.h"
26 #include "remove_equalities.h"
27 #include "polysign.h"
28 #include "verify.h"
29 #include "lexmin.h"
30 #include "param_util.h"
32 #undef CS /* for Solaris 10 */
34 #ifdef NTL_STD_CXX
35 using namespace NTL;
36 #endif
38 using std::vector;
39 using std::map;
40 using std::cerr;
41 using std::cout;
42 using std::endl;
43 using std::ostream;
45 #define EMPTINESS_CHECK (BV_OPT_LAST+1)
46 #define NO_REDUCTION (BV_OPT_LAST+2)
48 struct argp_option argp_options[] = {
49 { "emptiness-check", EMPTINESS_CHECK, "[none|count]", 0 },
50 { "no-reduction", NO_REDUCTION, 0, 0 },
51 { 0 }
54 static error_t parse_opt(int key, char *arg, struct argp_state *state)
56 struct lexmin_options *options = (struct lexmin_options *)(state->input);
57 struct barvinok_options *bv_options = options->verify.barvinok;
59 switch (key) {
60 case ARGP_KEY_INIT:
61 state->child_inputs[0] = options->verify.barvinok;
62 state->child_inputs[1] = &options->verify;
63 options->emptiness_check = BV_LEXMIN_EMPTINESS_CHECK_SAMPLE;
64 options->reduce = 1;
65 break;
66 case EMPTINESS_CHECK:
67 if (!strcmp(arg, "none"))
68 options->emptiness_check = BV_LEXMIN_EMPTINESS_CHECK_NONE;
69 else if (!strcmp(arg, "count")) {
70 options->emptiness_check = BV_LEXMIN_EMPTINESS_CHECK_COUNT;
71 bv_options->count_sample_infinite = 0;
73 break;
74 case NO_REDUCTION:
75 options->reduce = 0;
76 break;
77 default:
78 return ARGP_ERR_UNKNOWN;
80 return 0;
83 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
85 static int type_offset(enode *p)
87 return p->type == fractional ? 1 :
88 p->type == flooring ? 1 : 0;
91 void compute_evalue(evalue *e, Value *val, Value *res)
93 double d = compute_evalue(e, val);
94 if (d > 0)
95 d += .25;
96 else
97 d -= .25;
98 value_set_double(*res, d);
101 struct indicator_term {
102 int sign;
103 int pos; /* number of rational vertex */
104 int n; /* number of cone associated to given rational vertex */
105 mat_ZZ den;
106 evalue **vertex;
108 indicator_term(unsigned dim, int pos) {
109 den.SetDims(0, dim);
110 vertex = new evalue* [dim];
111 this->pos = pos;
112 n = -1;
113 sign = 0;
115 indicator_term(unsigned dim, int pos, int n) {
116 den.SetDims(dim, dim);
117 vertex = new evalue* [dim];
118 this->pos = pos;
119 this->n = n;
121 indicator_term(const indicator_term& src) {
122 sign = src.sign;
123 pos = src.pos;
124 n = src.n;
125 den = src.den;
126 unsigned dim = den.NumCols();
127 vertex = new evalue* [dim];
128 for (int i = 0; i < dim; ++i) {
129 vertex[i] = new evalue();
130 value_init(vertex[i]->d);
131 evalue_copy(vertex[i], src.vertex[i]);
134 void swap(indicator_term *other) {
135 int tmp;
136 tmp = sign; sign = other->sign; other->sign = tmp;
137 tmp = pos; pos = other->pos; other->pos = tmp;
138 tmp = n; n = other->n; other->n = tmp;
139 mat_ZZ tmp_den = den; den = other->den; other->den = tmp_den;
140 unsigned dim = den.NumCols();
141 for (int i = 0; i < dim; ++i) {
142 evalue *tmp = vertex[i];
143 vertex[i] = other->vertex[i];
144 other->vertex[i] = tmp;
147 ~indicator_term() {
148 unsigned dim = den.NumCols();
149 for (int i = 0; i < dim; ++i) {
150 free_evalue_refs(vertex[i]);
151 delete vertex[i];
153 delete [] vertex;
155 void print(ostream& os, char **p) const;
156 void substitute(Matrix *T);
157 void normalize();
158 void substitute(evalue *fract, evalue *val);
159 void substitute(int pos, evalue *val);
160 void reduce_in_domain(Polyhedron *D);
161 bool is_opposite(const indicator_term *neg) const;
162 vec_ZZ eval(Value *val) const {
163 vec_ZZ v;
164 unsigned dim = den.NumCols();
165 v.SetLength(dim);
166 Value tmp;
167 value_init(tmp);
168 for (int i = 0; i < dim; ++i) {
169 compute_evalue(vertex[i], val, &tmp);
170 value2zz(tmp, v[i]);
172 value_clear(tmp);
173 return v;
177 static int evalue_rational_cmp(const evalue *e1, const evalue *e2)
179 int r;
180 Value m;
181 Value m2;
182 value_init(m);
183 value_init(m2);
185 assert(value_notzero_p(e1->d));
186 assert(value_notzero_p(e2->d));
187 value_multiply(m, e1->x.n, e2->d);
188 value_multiply(m2, e2->x.n, e1->d);
189 if (value_lt(m, m2))
190 r = -1;
191 else if (value_gt(m, m2))
192 r = 1;
193 else
194 r = 0;
195 value_clear(m);
196 value_clear(m2);
198 return r;
201 static int evalue_cmp(const evalue *e1, const evalue *e2)
203 if (value_notzero_p(e1->d)) {
204 if (value_zero_p(e2->d))
205 return -1;
206 return evalue_rational_cmp(e1, e2);
208 if (value_notzero_p(e2->d))
209 return 1;
210 if (e1->x.p->type != e2->x.p->type)
211 return e1->x.p->type - e2->x.p->type;
212 if (e1->x.p->size != e2->x.p->size)
213 return e1->x.p->size - e2->x.p->size;
214 if (e1->x.p->pos != e2->x.p->pos)
215 return e1->x.p->pos - e2->x.p->pos;
216 assert(e1->x.p->type == polynomial ||
217 e1->x.p->type == fractional ||
218 e1->x.p->type == flooring);
219 for (int i = 0; i < e1->x.p->size; ++i) {
220 int s = evalue_cmp(&e1->x.p->arr[i], &e2->x.p->arr[i]);
221 if (s)
222 return s;
224 return 0;
227 void evalue_length(evalue *e, int len[2])
229 len[0] = 0;
230 len[1] = 0;
232 while (value_zero_p(e->d)) {
233 assert(e->x.p->type == polynomial ||
234 e->x.p->type == fractional ||
235 e->x.p->type == flooring);
236 if (e->x.p->type == polynomial)
237 len[1]++;
238 else
239 len[0]++;
240 int offset = type_offset(e->x.p);
241 assert(e->x.p->size == offset+2);
242 e = &e->x.p->arr[offset];
246 static bool it_smaller(const indicator_term* it1, const indicator_term* it2)
248 if (it1 == it2)
249 return false;
250 int len1[2], len2[2];
251 unsigned dim = it1->den.NumCols();
252 for (int i = 0; i < dim; ++i) {
253 evalue_length(it1->vertex[i], len1);
254 evalue_length(it2->vertex[i], len2);
255 if (len1[0] != len2[0])
256 return len1[0] < len2[0];
257 if (len1[1] != len2[1])
258 return len1[1] < len2[1];
260 if (it1->pos != it2->pos)
261 return it1->pos < it2->pos;
262 if (it1->n != it2->n)
263 return it1->n < it2->n;
264 int s = lex_cmp(it1->den, it2->den);
265 if (s)
266 return s < 0;
267 for (int i = 0; i < dim; ++i) {
268 s = evalue_cmp(it1->vertex[i], it2->vertex[i]);
269 if (s)
270 return s < 0;
272 assert(it1->sign != 0);
273 assert(it2->sign != 0);
274 if (it1->sign != it2->sign)
275 return it1->sign > 0;
276 return it1 < it2;
279 struct smaller_it {
280 static const int requires_resort;
281 bool operator()(const indicator_term* it1, const indicator_term* it2) const {
282 return it_smaller(it1, it2);
285 const int smaller_it::requires_resort = 1;
287 struct smaller_it_p {
288 static const int requires_resort;
289 bool operator()(const indicator_term* it1, const indicator_term* it2) const {
290 return it1 < it2;
293 const int smaller_it_p::requires_resort = 0;
295 /* Returns true if this and neg are opposite using the knowledge
296 * that they have the same numerator.
297 * In particular, we check that the signs are different and that
298 * the denominator is the same.
300 bool indicator_term::is_opposite(const indicator_term *neg) const
302 if (sign + neg->sign != 0)
303 return false;
304 if (den != neg->den)
305 return false;
306 return true;
309 void indicator_term::reduce_in_domain(Polyhedron *D)
311 for (int k = 0; k < den.NumCols(); ++k) {
312 reduce_evalue_in_domain(vertex[k], D);
313 if (evalue_range_reduction_in_domain(vertex[k], D))
314 reduce_evalue(vertex[k]);
318 void indicator_term::print(ostream& os, char **p) const
320 unsigned dim = den.NumCols();
321 unsigned factors = den.NumRows();
322 if (sign == 0)
323 os << " s ";
324 else if (sign > 0)
325 os << " + ";
326 else
327 os << " - ";
328 os << '[';
329 for (int i = 0; i < dim; ++i) {
330 if (i)
331 os << ',';
332 evalue_print(os, vertex[i], p);
334 os << ']';
335 for (int i = 0; i < factors; ++i) {
336 os << " + t" << i << "*[";
337 for (int j = 0; j < dim; ++j) {
338 if (j)
339 os << ',';
340 os << den[i][j];
342 os << "]";
344 os << " ((" << pos << ", " << n << ", " << (void*)this << "))";
347 /* Perform the substitution specified by T on the variables.
348 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
349 * The substitution is performed as in gen_fun::substitute
351 void indicator_term::substitute(Matrix *T)
353 unsigned dim = den.NumCols();
354 unsigned nparam = T->NbColumns - dim - 1;
355 unsigned newdim = T->NbRows - nparam - 1;
356 evalue **newvertex;
357 mat_ZZ trans;
358 matrix2zz(T, trans, newdim, dim);
359 trans = transpose(trans);
360 den *= trans;
361 newvertex = new evalue* [newdim];
363 vec_ZZ v;
364 v.SetLength(nparam+1);
366 evalue factor;
367 value_init(factor.d);
368 value_set_si(factor.d, 1);
369 value_init(factor.x.n);
370 for (int i = 0; i < newdim; ++i) {
371 values2zz(T->p[i]+dim, v, nparam+1);
372 newvertex[i] = multi_monom(v);
374 for (int j = 0; j < dim; ++j) {
375 if (value_zero_p(T->p[i][j]))
376 continue;
377 evalue term;
378 value_init(term.d);
379 evalue_copy(&term, vertex[j]);
380 value_assign(factor.x.n, T->p[i][j]);
381 emul(&factor, &term);
382 eadd(&term, newvertex[i]);
383 free_evalue_refs(&term);
386 free_evalue_refs(&factor);
387 for (int i = 0; i < dim; ++i) {
388 free_evalue_refs(vertex[i]);
389 delete vertex[i];
391 delete [] vertex;
392 vertex = newvertex;
395 static void evalue_add_constant(evalue *e, ZZ v)
397 Value tmp;
398 value_init(tmp);
400 /* go down to constant term */
401 while (value_zero_p(e->d))
402 e = &e->x.p->arr[type_offset(e->x.p)];
403 /* and add v */
404 zz2value(v, tmp);
405 value_multiply(tmp, tmp, e->d);
406 value_addto(e->x.n, e->x.n, tmp);
408 value_clear(tmp);
411 /* Make all powers in denominator lexico-positive */
412 void indicator_term::normalize()
414 vec_ZZ extra_vertex;
415 extra_vertex.SetLength(den.NumCols());
416 for (int r = 0; r < den.NumRows(); ++r) {
417 for (int k = 0; k < den.NumCols(); ++k) {
418 if (den[r][k] == 0)
419 continue;
420 if (den[r][k] > 0)
421 break;
422 sign = -sign;
423 den[r] = -den[r];
424 extra_vertex += den[r];
425 break;
428 for (int k = 0; k < extra_vertex.length(); ++k)
429 if (extra_vertex[k] != 0)
430 evalue_add_constant(vertex[k], extra_vertex[k]);
433 static void substitute(evalue *e, evalue *fract, evalue *val)
435 evalue *t;
437 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
438 if (t->x.p->type == fractional && eequal(&t->x.p->arr[0], fract))
439 break;
441 if (value_notzero_p(t->d))
442 return;
444 free_evalue_refs(&t->x.p->arr[0]);
445 evalue *term = &t->x.p->arr[2];
446 enode *p = t->x.p;
447 value_clear(t->d);
448 *t = t->x.p->arr[1];
450 emul(val, term);
451 eadd(term, e);
452 free_evalue_refs(term);
453 free(p);
455 reduce_evalue(e);
458 void indicator_term::substitute(evalue *fract, evalue *val)
460 unsigned dim = den.NumCols();
461 for (int i = 0; i < dim; ++i) {
462 ::substitute(vertex[i], fract, val);
466 static void substitute(evalue *e, int pos, evalue *val)
468 evalue *t;
470 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
471 if (t->x.p->type == polynomial && t->x.p->pos == pos)
472 break;
474 if (value_notzero_p(t->d))
475 return;
477 evalue *term = &t->x.p->arr[1];
478 enode *p = t->x.p;
479 value_clear(t->d);
480 *t = t->x.p->arr[0];
482 emul(val, term);
483 eadd(term, e);
484 free_evalue_refs(term);
485 free(p);
487 reduce_evalue(e);
490 void indicator_term::substitute(int pos, evalue *val)
492 unsigned dim = den.NumCols();
493 for (int i = 0; i < dim; ++i) {
494 ::substitute(vertex[i], pos, val);
498 struct indicator_constructor : public signed_cone_consumer,
499 public vertex_decomposer {
500 vec_ZZ vertex;
501 vector<indicator_term*> *terms;
502 Matrix *T; /* Transformation to original space */
503 int nbV;
504 int pos;
505 int n;
507 indicator_constructor(Polyhedron *P, unsigned dim, Param_Polyhedron *PP,
508 Matrix *T) :
509 vertex_decomposer(PP, *this), T(T), nbV(PP->nbV) {
510 vertex.SetLength(dim);
511 terms = new vector<indicator_term*>[PP->nbV];
513 ~indicator_constructor() {
514 for (int i = 0; i < nbV; ++i)
515 for (int j = 0; j < terms[i].size(); ++j)
516 delete terms[i][j];
517 delete [] terms;
519 void substitute(Matrix *T);
520 void normalize();
521 void print(ostream& os, char **p);
523 virtual void handle(const signed_cone& sc, barvinok_options *options);
524 void decompose_at_vertex(Param_Vertices *V, int _i,
525 barvinok_options *options) {
526 pos = _i;
527 n = 0;
528 vertex_decomposer::decompose_at_vertex(V, _i, options);
532 void indicator_constructor::handle(const signed_cone& sc, barvinok_options *options)
534 assert(sc.det == 1);
535 unsigned dim = vertex.length();
537 assert(sc.rays.NumRows() == dim);
539 indicator_term *term = new indicator_term(dim, pos, n++);
540 term->sign = sc.sign;
541 terms[vert].push_back(term);
543 lattice_point(V, sc.rays, vertex, term->vertex, options);
545 term->den = sc.rays;
546 for (int r = 0; r < dim; ++r) {
547 for (int j = 0; j < dim; ++j) {
548 if (term->den[r][j] == 0)
549 continue;
550 if (term->den[r][j] > 0)
551 break;
552 term->sign = -term->sign;
553 term->den[r] = -term->den[r];
554 vertex += term->den[r];
555 break;
559 for (int i = 0; i < dim; ++i) {
560 if (!term->vertex[i]) {
561 term->vertex[i] = new evalue();
562 value_init(term->vertex[i]->d);
563 value_init(term->vertex[i]->x.n);
564 zz2value(vertex[i], term->vertex[i]->x.n);
565 value_set_si(term->vertex[i]->d, 1);
566 continue;
568 if (vertex[i] == 0)
569 continue;
570 evalue_add_constant(term->vertex[i], vertex[i]);
573 if (T) {
574 term->substitute(T);
575 term->normalize();
578 lex_order_rows(term->den);
581 void indicator_constructor::print(ostream& os, char **p)
583 for (int i = 0; i < PP->nbV; ++i)
584 for (int j = 0; j < terms[i].size(); ++j) {
585 os << "i: " << i << ", j: " << j << endl;
586 terms[i][j]->print(os, p);
587 os << endl;
591 void indicator_constructor::normalize()
593 for (int i = 0; i < PP->nbV; ++i)
594 for (int j = 0; j < terms[i].size(); ++j) {
595 vec_ZZ vertex;
596 vertex.SetLength(terms[i][j]->den.NumCols());
597 for (int r = 0; r < terms[i][j]->den.NumRows(); ++r) {
598 for (int k = 0; k < terms[i][j]->den.NumCols(); ++k) {
599 if (terms[i][j]->den[r][k] == 0)
600 continue;
601 if (terms[i][j]->den[r][k] > 0)
602 break;
603 terms[i][j]->sign = -terms[i][j]->sign;
604 terms[i][j]->den[r] = -terms[i][j]->den[r];
605 vertex += terms[i][j]->den[r];
606 break;
609 lex_order_rows(terms[i][j]->den);
610 for (int k = 0; k < vertex.length(); ++k)
611 if (vertex[k] != 0)
612 evalue_add_constant(terms[i][j]->vertex[k], vertex[k]);
616 struct order_cache_el {
617 vector<evalue *> e;
618 order_cache_el copy() const {
619 order_cache_el n;
620 for (int i = 0; i < e.size(); ++i) {
621 evalue *c = new evalue;
622 value_init(c->d);
623 evalue_copy(c, e[i]);
624 n.e.push_back(c);
626 return n;
628 void free() {
629 for (int i = 0; i < e.size(); ++i) {
630 free_evalue_refs(e[i]);
631 delete e[i];
634 void negate() {
635 evalue mone;
636 value_init(mone.d);
637 evalue_set_si(&mone, -1, 1);
638 for (int i = 0; i < e.size(); ++i)
639 emul(&mone, e[i]);
640 free_evalue_refs(&mone);
642 void print(ostream& os, char **p);
645 void order_cache_el::print(ostream& os, char **p)
647 os << "[";
648 for (int i = 0; i < e.size(); ++i) {
649 if (i)
650 os << ",";
651 evalue_print(os, e[i], p);
653 os << "]";
656 struct order_cache {
657 vector<order_cache_el> lt;
658 vector<order_cache_el> le;
659 vector<order_cache_el> unknown;
661 void clear_transients() {
662 for (int i = 0; i < le.size(); ++i)
663 le[i].free();
664 for (int i = 0; i < unknown.size(); ++i)
665 unknown[i].free();
666 le.resize(0);
667 unknown.resize(0);
669 ~order_cache() {
670 clear_transients();
671 for (int i = 0; i < lt.size(); ++i)
672 lt[i].free();
673 lt.resize(0);
675 void add(order_cache_el& cache_el, order_sign sign);
676 order_sign check_lt(vector<order_cache_el>* list,
677 const indicator_term *a, const indicator_term *b,
678 order_cache_el& cache_el);
679 order_sign check_lt(const indicator_term *a, const indicator_term *b,
680 order_cache_el& cache_el);
681 order_sign check_direct(const indicator_term *a, const indicator_term *b,
682 order_cache_el& cache_el);
683 order_sign check(const indicator_term *a, const indicator_term *b,
684 order_cache_el& cache_el);
685 void copy(const order_cache& cache);
686 void print(ostream& os, char **p);
689 void order_cache::copy(const order_cache& cache)
691 for (int i = 0; i < cache.lt.size(); ++i) {
692 order_cache_el n = cache.lt[i].copy();
693 add(n, order_lt);
697 void order_cache::add(order_cache_el& cache_el, order_sign sign)
699 if (sign == order_lt) {
700 lt.push_back(cache_el);
701 } else if (sign == order_gt) {
702 cache_el.negate();
703 lt.push_back(cache_el);
704 } else if (sign == order_le) {
705 le.push_back(cache_el);
706 } else if (sign == order_ge) {
707 cache_el.negate();
708 le.push_back(cache_el);
709 } else if (sign == order_unknown) {
710 unknown.push_back(cache_el);
711 } else {
712 assert(sign == order_eq);
713 cache_el.free();
715 return;
718 /* compute a - b */
719 static evalue *ediff(const evalue *a, const evalue *b)
721 evalue mone;
722 value_init(mone.d);
723 evalue_set_si(&mone, -1, 1);
724 evalue *diff = new evalue;
725 value_init(diff->d);
726 evalue_copy(diff, b);
727 emul(&mone, diff);
728 eadd(a, diff);
729 reduce_evalue(diff);
730 free_evalue_refs(&mone);
731 return diff;
734 static bool evalue_first_difference(const evalue *e1, const evalue *e2,
735 const evalue **d1, const evalue **d2)
737 *d1 = e1;
738 *d2 = e2;
740 if (value_ne(e1->d, e2->d))
741 return true;
743 if (value_notzero_p(e1->d)) {
744 if (value_eq(e1->x.n, e2->x.n))
745 return false;
746 return true;
748 if (e1->x.p->type != e2->x.p->type)
749 return true;
750 if (e1->x.p->size != e2->x.p->size)
751 return true;
752 if (e1->x.p->pos != e2->x.p->pos)
753 return true;
755 assert(e1->x.p->type == polynomial ||
756 e1->x.p->type == fractional ||
757 e1->x.p->type == flooring);
758 int offset = type_offset(e1->x.p);
759 assert(e1->x.p->size == offset+2);
760 for (int i = 0; i < e1->x.p->size; ++i)
761 if (i != type_offset(e1->x.p) &&
762 !eequal(&e1->x.p->arr[i], &e2->x.p->arr[i]))
763 return true;
765 return evalue_first_difference(&e1->x.p->arr[offset],
766 &e2->x.p->arr[offset], d1, d2);
769 static order_sign evalue_diff_constant_sign(const evalue *e1, const evalue *e2)
771 if (!evalue_first_difference(e1, e2, &e1, &e2))
772 return order_eq;
773 if (value_zero_p(e1->d) || value_zero_p(e2->d))
774 return order_undefined;
775 int s = evalue_rational_cmp(e1, e2);
776 if (s < 0)
777 return order_lt;
778 else if (s > 0)
779 return order_gt;
780 else
781 return order_eq;
784 order_sign order_cache::check_lt(vector<order_cache_el>* list,
785 const indicator_term *a, const indicator_term *b,
786 order_cache_el& cache_el)
788 order_sign sign = order_undefined;
789 for (int i = 0; i < list->size(); ++i) {
790 int j;
791 for (j = cache_el.e.size(); j < (*list)[i].e.size(); ++j)
792 cache_el.e.push_back(ediff(a->vertex[j], b->vertex[j]));
793 for (j = 0; j < (*list)[i].e.size(); ++j) {
794 order_sign diff_sign;
795 diff_sign = evalue_diff_constant_sign((*list)[i].e[j], cache_el.e[j]);
796 if (diff_sign == order_gt) {
797 sign = order_lt;
798 break;
799 } else if (diff_sign == order_lt)
800 break;
801 else if (diff_sign == order_undefined)
802 break;
803 else
804 assert(diff_sign == order_eq);
806 if (j == (*list)[i].e.size())
807 sign = list == &lt ? order_lt : order_le;
808 if (sign != order_undefined)
809 break;
811 return sign;
814 order_sign order_cache::check_direct(const indicator_term *a,
815 const indicator_term *b,
816 order_cache_el& cache_el)
818 order_sign sign = check_lt(&lt, a, b, cache_el);
819 if (sign != order_undefined)
820 return sign;
821 sign = check_lt(&le, a, b, cache_el);
822 if (sign != order_undefined)
823 return sign;
825 for (int i = 0; i < unknown.size(); ++i) {
826 int j;
827 for (j = cache_el.e.size(); j < unknown[i].e.size(); ++j)
828 cache_el.e.push_back(ediff(a->vertex[j], b->vertex[j]));
829 for (j = 0; j < unknown[i].e.size(); ++j) {
830 if (!eequal(unknown[i].e[j], cache_el.e[j]))
831 break;
833 if (j == unknown[i].e.size()) {
834 sign = order_unknown;
835 break;
838 return sign;
841 order_sign order_cache::check(const indicator_term *a, const indicator_term *b,
842 order_cache_el& cache_el)
844 order_sign sign = check_direct(a, b, cache_el);
845 if (sign != order_undefined)
846 return sign;
847 int size = cache_el.e.size();
848 cache_el.negate();
849 sign = check_direct(a, b, cache_el);
850 cache_el.negate();
851 assert(cache_el.e.size() == size);
852 if (sign == order_undefined)
853 return sign;
854 if (sign == order_lt)
855 sign = order_gt;
856 else if (sign == order_le)
857 sign = order_ge;
858 else
859 assert(sign == order_unknown);
860 return sign;
863 struct indicator;
865 struct partial_order {
866 indicator *ind;
868 std::set<const indicator_term *, smaller_it > head;
869 map<const indicator_term *, int, smaller_it > pred;
870 map<const indicator_term *, vector<const indicator_term * >, smaller_it > lt;
871 map<const indicator_term *, vector<const indicator_term * >, smaller_it > le;
872 map<const indicator_term *, vector<const indicator_term * >, smaller_it > eq;
874 map<const indicator_term *, vector<const indicator_term * >, smaller_it > pending;
876 order_cache cache;
878 partial_order(indicator *ind) : ind(ind) {}
879 void copy(const partial_order& order,
880 map< const indicator_term *, indicator_term * > old2new);
881 void resort() {
882 map<const indicator_term *, vector<const indicator_term * > >::iterator i;
883 map<const indicator_term *, int >::iterator j;
884 std::set<const indicator_term *>::iterator k;
886 if (head.key_comp().requires_resort) {
887 typeof(head) new_head;
888 for (k = head.begin(); k != head.end(); ++k)
889 new_head.insert(*k);
890 head.swap(new_head);
891 new_head.clear();
894 if (pred.key_comp().requires_resort) {
895 typeof(pred) new_pred;
896 for (j = pred.begin(); j != pred.end(); ++j)
897 new_pred[(*j).first] = (*j).second;
898 pred.swap(new_pred);
899 new_pred.clear();
902 if (lt.key_comp().requires_resort) {
903 typeof(lt) m;
904 for (i = lt.begin(); i != lt.end(); ++i)
905 m[(*i).first] = (*i).second;
906 lt.swap(m);
907 m.clear();
910 if (le.key_comp().requires_resort) {
911 typeof(le) m;
912 for (i = le.begin(); i != le.end(); ++i)
913 m[(*i).first] = (*i).second;
914 le.swap(m);
915 m.clear();
918 if (eq.key_comp().requires_resort) {
919 typeof(eq) m;
920 for (i = eq.begin(); i != eq.end(); ++i)
921 m[(*i).first] = (*i).second;
922 eq.swap(m);
923 m.clear();
926 if (pending.key_comp().requires_resort) {
927 typeof(pending) m;
928 for (i = pending.begin(); i != pending.end(); ++i)
929 m[(*i).first] = (*i).second;
930 pending.swap(m);
931 m.clear();
935 order_sign compare(const indicator_term *a, const indicator_term *b);
936 void set_equal(const indicator_term *a, const indicator_term *b);
937 void unset_le(const indicator_term *a, const indicator_term *b);
938 void dec_pred(const indicator_term *it) {
939 if (--pred[it] == 0) {
940 pred.erase(it);
941 head.insert(it);
944 void inc_pred(const indicator_term *it) {
945 if (head.find(it) != head.end())
946 head.erase(it);
947 pred[it]++;
950 bool compared(const indicator_term* a, const indicator_term* b);
951 void add(const indicator_term* it, std::set<const indicator_term *> *filter);
952 void remove(const indicator_term* it);
954 void print(ostream& os, char **p);
956 /* replace references to orig to references to replacement */
957 void replace(const indicator_term* orig, indicator_term* replacement);
958 void sanity_check() const;
961 /* We actually replace the contents of orig by that of replacement,
962 * but we have to be careful since replacing the content changes
963 * the order of orig in the maps.
965 void partial_order::replace(const indicator_term* orig, indicator_term* replacement)
967 std::set<const indicator_term *>::iterator k;
968 k = head.find(orig);
969 bool is_head = k != head.end();
970 int orig_pred;
971 if (is_head) {
972 head.erase(orig);
973 } else {
974 orig_pred = pred[orig];
975 pred.erase(orig);
977 vector<const indicator_term * > orig_lt;
978 vector<const indicator_term * > orig_le;
979 vector<const indicator_term * > orig_eq;
980 vector<const indicator_term * > orig_pending;
981 map<const indicator_term *, vector<const indicator_term * > >::iterator i;
982 bool in_lt = ((i = lt.find(orig)) != lt.end());
983 if (in_lt) {
984 orig_lt = (*i).second;
985 lt.erase(orig);
987 bool in_le = ((i = le.find(orig)) != le.end());
988 if (in_le) {
989 orig_le = (*i).second;
990 le.erase(orig);
992 bool in_eq = ((i = eq.find(orig)) != eq.end());
993 if (in_eq) {
994 orig_eq = (*i).second;
995 eq.erase(orig);
997 bool in_pending = ((i = pending.find(orig)) != pending.end());
998 if (in_pending) {
999 orig_pending = (*i).second;
1000 pending.erase(orig);
1002 indicator_term *old = const_cast<indicator_term *>(orig);
1003 old->swap(replacement);
1004 if (is_head)
1005 head.insert(old);
1006 else
1007 pred[old] = orig_pred;
1008 if (in_lt)
1009 lt[old] = orig_lt;
1010 if (in_le)
1011 le[old] = orig_le;
1012 if (in_eq)
1013 eq[old] = orig_eq;
1014 if (in_pending)
1015 pending[old] = orig_pending;
1018 void partial_order::unset_le(const indicator_term *a, const indicator_term *b)
1020 vector<const indicator_term *>::iterator i;
1021 i = find(le[a].begin(), le[a].end(), b);
1022 le[a].erase(i);
1023 if (le[a].size() == 0)
1024 le.erase(a);
1025 dec_pred(b);
1026 i = find(pending[a].begin(), pending[a].end(), b);
1027 if (i != pending[a].end())
1028 pending[a].erase(i);
1031 void partial_order::set_equal(const indicator_term *a, const indicator_term *b)
1033 if (eq[a].size() == 0)
1034 eq[a].push_back(a);
1035 if (eq[b].size() == 0)
1036 eq[b].push_back(b);
1037 a = eq[a][0];
1038 b = eq[b][0];
1039 assert(a != b);
1040 if (pred.key_comp()(b, a)) {
1041 const indicator_term *c = a;
1042 a = b;
1043 b = c;
1046 const indicator_term *base = a;
1048 map<const indicator_term *, vector<const indicator_term * > >::iterator i;
1050 for (int j = 0; j < eq[b].size(); ++j) {
1051 eq[base].push_back(eq[b][j]);
1052 eq[eq[b][j]][0] = base;
1054 eq[b].resize(1);
1056 i = lt.find(b);
1057 if (i != lt.end()) {
1058 for (int j = 0; j < lt[b].size(); ++j) {
1059 if (find(eq[base].begin(), eq[base].end(), lt[b][j]) != eq[base].end())
1060 dec_pred(lt[b][j]);
1061 else if (find(lt[base].begin(), lt[base].end(), lt[b][j])
1062 != lt[base].end())
1063 dec_pred(lt[b][j]);
1064 else
1065 lt[base].push_back(lt[b][j]);
1067 lt.erase(b);
1070 i = le.find(b);
1071 if (i != le.end()) {
1072 for (int j = 0; j < le[b].size(); ++j) {
1073 if (find(eq[base].begin(), eq[base].end(), le[b][j]) != eq[base].end())
1074 dec_pred(le[b][j]);
1075 else if (find(le[base].begin(), le[base].end(), le[b][j])
1076 != le[base].end())
1077 dec_pred(le[b][j]);
1078 else
1079 le[base].push_back(le[b][j]);
1081 le.erase(b);
1084 i = pending.find(base);
1085 if (i != pending.end()) {
1086 vector<const indicator_term * > old = pending[base];
1087 pending[base].clear();
1088 for (int j = 0; j < old.size(); ++j) {
1089 if (find(eq[base].begin(), eq[base].end(), old[j]) == eq[base].end())
1090 pending[base].push_back(old[j]);
1094 i = pending.find(b);
1095 if (i != pending.end()) {
1096 for (int j = 0; j < pending[b].size(); ++j) {
1097 if (find(eq[base].begin(), eq[base].end(), pending[b][j]) == eq[base].end())
1098 pending[base].push_back(pending[b][j]);
1100 pending.erase(b);
1104 void partial_order::copy(const partial_order& order,
1105 map< const indicator_term *, indicator_term * > old2new)
1107 cache.copy(order.cache);
1109 map<const indicator_term *, vector<const indicator_term * > >::const_iterator i;
1110 map<const indicator_term *, int >::const_iterator j;
1111 std::set<const indicator_term *>::const_iterator k;
1113 for (k = order.head.begin(); k != order.head.end(); ++k)
1114 head.insert(old2new[*k]);
1116 for (j = order.pred.begin(); j != order.pred.end(); ++j)
1117 pred[old2new[(*j).first]] = (*j).second;
1119 for (i = order.lt.begin(); i != order.lt.end(); ++i) {
1120 for (int j = 0; j < (*i).second.size(); ++j)
1121 lt[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1123 for (i = order.le.begin(); i != order.le.end(); ++i) {
1124 for (int j = 0; j < (*i).second.size(); ++j)
1125 le[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1127 for (i = order.eq.begin(); i != order.eq.end(); ++i) {
1128 for (int j = 0; j < (*i).second.size(); ++j)
1129 eq[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1131 for (i = order.pending.begin(); i != order.pending.end(); ++i) {
1132 for (int j = 0; j < (*i).second.size(); ++j)
1133 pending[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1137 struct max_term {
1138 EDomain *domain;
1139 vector<evalue *> max;
1141 void print(ostream& os, char **p, barvinok_options *options) const;
1142 void substitute(Matrix *T, barvinok_options *options);
1143 Vector *eval(Value *val, unsigned MaxRays) const;
1145 ~max_term() {
1146 for (int i = 0; i < max.size(); ++i) {
1147 free_evalue_refs(max[i]);
1148 delete max[i];
1150 delete domain;
1155 * Project on first dim dimensions
1157 Polyhedron* Polyhedron_Project_Initial(Polyhedron *P, int dim)
1159 int i;
1160 Matrix *T;
1161 Polyhedron *I;
1163 if (P->Dimension == dim)
1164 return Polyhedron_Copy(P);
1166 T = Matrix_Alloc(dim+1, P->Dimension+1);
1167 for (i = 0; i < dim; ++i)
1168 value_set_si(T->p[i][i], 1);
1169 value_set_si(T->p[dim][P->Dimension], 1);
1170 I = Polyhedron_Image(P, T, P->NbConstraints);
1171 Matrix_Free(T);
1172 return I;
1175 struct indicator {
1176 vector<indicator_term*> term;
1177 indicator_constructor& ic;
1178 partial_order order;
1179 EDomain *D;
1180 Polyhedron *P;
1181 Param_Domain *PD;
1182 lexmin_options *options;
1183 vector<evalue *> substitutions;
1185 indicator(indicator_constructor& ic, Param_Domain *PD, EDomain *D,
1186 lexmin_options *options) :
1187 ic(ic), PD(PD), D(D), order(this), options(options), P(NULL) {}
1188 indicator(const indicator& ind, EDomain *D) :
1189 ic(ind.ic), PD(ind.PD), D(NULL), order(this), options(ind.options),
1190 P(Polyhedron_Copy(ind.P)) {
1191 map< const indicator_term *, indicator_term * > old2new;
1192 for (int i = 0; i < ind.term.size(); ++i) {
1193 indicator_term *it = new indicator_term(*ind.term[i]);
1194 old2new[ind.term[i]] = it;
1195 term.push_back(it);
1197 order.copy(ind.order, old2new);
1198 set_domain(D);
1200 ~indicator() {
1201 for (int i = 0; i < term.size(); ++i)
1202 delete term[i];
1203 if (D)
1204 delete D;
1205 if (P)
1206 Polyhedron_Free(P);
1209 void set_domain(EDomain *D) {
1210 order.cache.clear_transients();
1211 if (this->D)
1212 delete this->D;
1213 this->D = D;
1214 int nparam = ic.PP->Constraints->NbColumns-2 - ic.vertex.length();
1215 if (options->reduce) {
1216 Polyhedron *Q = Polyhedron_Project_Initial(D->D, nparam);
1217 Q = DomainConstraintSimplify(Q, options->verify.barvinok->MaxRays);
1218 if (!P || !PolyhedronIncludes(Q, P))
1219 reduce_in_domain(Q);
1220 if (P)
1221 Polyhedron_Free(P);
1222 P = Q;
1223 order.resort();
1227 void add(const indicator_term* it);
1228 void remove(const indicator_term* it);
1229 void remove_initial_rational_vertices();
1230 void expand_rational_vertex(const indicator_term *initial);
1232 void print(ostream& os, char **p);
1233 void simplify();
1234 void peel(int i, int j);
1235 void combine(const indicator_term *a, const indicator_term *b);
1236 void add_substitution(evalue *equation);
1237 void perform_pending_substitutions();
1238 void reduce_in_domain(Polyhedron *D);
1239 bool handle_equal_numerators(const indicator_term *base);
1241 max_term* create_max_term(const indicator_term *it);
1242 private:
1243 void substitute(evalue *equation);
1246 void partial_order::sanity_check() const
1248 map<const indicator_term *, vector<const indicator_term * > >::const_iterator i;
1249 map<const indicator_term *, vector<const indicator_term * > >::const_iterator prev;
1250 map<const indicator_term *, vector<const indicator_term * > >::const_iterator l;
1251 map<const indicator_term *, int >::const_iterator k, prev_k;
1253 for (k = pred.begin(); k != pred.end(); prev_k = k, ++k)
1254 if (k != pred.begin())
1255 assert(pred.key_comp()((*prev_k).first, (*k).first));
1256 for (i = lt.begin(); i != lt.end(); prev = i, ++i) {
1257 vec_ZZ i_v;
1258 if (ind->D->sample)
1259 i_v = (*i).first->eval(ind->D->sample->p);
1260 if (i != lt.begin())
1261 assert(lt.key_comp()((*prev).first, (*i).first));
1262 l = eq.find((*i).first);
1263 if (l != eq.end())
1264 assert((*l).second.size() > 1);
1265 assert(head.find((*i).first) != head.end() ||
1266 pred.find((*i).first) != pred.end());
1267 for (int j = 0; j < (*i).second.size(); ++j) {
1268 k = pred.find((*i).second[j]);
1269 assert(k != pred.end());
1270 assert((*k).second != 0);
1271 if ((*i).first->sign != 0 &&
1272 (*i).second[j]->sign != 0 && ind->D->sample) {
1273 vec_ZZ j_v = (*i).second[j]->eval(ind->D->sample->p);
1274 assert(lex_cmp(i_v, j_v) < 0);
1278 for (i = le.begin(); i != le.end(); ++i) {
1279 assert((*i).second.size() > 0);
1280 assert(head.find((*i).first) != head.end() ||
1281 pred.find((*i).first) != pred.end());
1282 for (int j = 0; j < (*i).second.size(); ++j) {
1283 k = pred.find((*i).second[j]);
1284 assert(k != pred.end());
1285 assert((*k).second != 0);
1288 for (i = eq.begin(); i != eq.end(); ++i) {
1289 assert(head.find((*i).first) != head.end() ||
1290 pred.find((*i).first) != pred.end());
1291 assert((*i).second.size() >= 1);
1293 for (i = pending.begin(); i != pending.end(); ++i) {
1294 assert(head.find((*i).first) != head.end() ||
1295 pred.find((*i).first) != pred.end());
1296 for (int j = 0; j < (*i).second.size(); ++j)
1297 assert(head.find((*i).second[j]) != head.end() ||
1298 pred.find((*i).second[j]) != pred.end());
1302 max_term* indicator::create_max_term(const indicator_term *it)
1304 int dim = it->den.NumCols();
1305 int nparam = ic.PP->Constraints->NbColumns-2 - ic.vertex.length();
1306 max_term *maximum = new max_term;
1307 maximum->domain = new EDomain(D);
1308 for (int j = 0; j < dim; ++j) {
1309 evalue *E = new evalue;
1310 value_init(E->d);
1311 evalue_copy(E, it->vertex[j]);
1312 if (evalue_frac2floor_in_domain(E, D->D))
1313 reduce_evalue(E);
1314 maximum->max.push_back(E);
1316 return maximum;
1319 static order_sign evalue_sign(evalue *diff, EDomain *D, barvinok_options *options)
1321 order_sign sign = order_eq;
1322 evalue mone;
1323 value_init(mone.d);
1324 evalue_set_si(&mone, -1, 1);
1325 int len = 1 + D->D->Dimension + 1;
1326 Vector *c = Vector_Alloc(len);
1327 Matrix *T = Matrix_Alloc(2, len-1);
1329 int fract = evalue2constraint(D, diff, c->p, len);
1330 Vector_Copy(c->p+1, T->p[0], len-1);
1331 value_assign(T->p[1][len-2], c->p[0]);
1333 order_sign upper_sign = polyhedron_affine_sign(D->D, T, options);
1334 if (upper_sign == order_lt || !fract)
1335 sign = upper_sign;
1336 else {
1337 emul(&mone, diff);
1338 evalue2constraint(D, diff, c->p, len);
1339 emul(&mone, diff);
1340 Vector_Copy(c->p+1, T->p[0], len-1);
1341 value_assign(T->p[1][len-2], c->p[0]);
1343 order_sign neg_lower_sign = polyhedron_affine_sign(D->D, T, options);
1345 if (neg_lower_sign == order_lt)
1346 sign = order_gt;
1347 else if (neg_lower_sign == order_eq || neg_lower_sign == order_le) {
1348 if (upper_sign == order_eq || upper_sign == order_le)
1349 sign = order_eq;
1350 else
1351 sign = order_ge;
1352 } else {
1353 if (upper_sign == order_lt || upper_sign == order_le ||
1354 upper_sign == order_eq)
1355 sign = order_le;
1356 else
1357 sign = order_unknown;
1361 Matrix_Free(T);
1362 Vector_Free(c);
1363 free_evalue_refs(&mone);
1365 return sign;
1368 /* An auxiliary class that keeps a reference to an evalue
1369 * and frees it when it goes out of scope.
1371 struct temp_evalue {
1372 evalue *E;
1373 temp_evalue() : E(NULL) {}
1374 temp_evalue(evalue *e) : E(e) {}
1375 operator evalue* () const { return E; }
1376 evalue *operator=(evalue *e) {
1377 if (E) {
1378 free_evalue_refs(E);
1379 delete E;
1381 E = e;
1382 return E;
1384 ~temp_evalue() {
1385 if (E) {
1386 free_evalue_refs(E);
1387 delete E;
1392 static void substitute(vector<indicator_term*>& term, evalue *equation)
1394 evalue *fract = NULL;
1395 evalue *val = new evalue;
1396 value_init(val->d);
1397 evalue_copy(val, equation);
1399 evalue factor;
1400 value_init(factor.d);
1401 value_init(factor.x.n);
1403 evalue *e;
1404 for (e = val; value_zero_p(e->d) && e->x.p->type != fractional;
1405 e = &e->x.p->arr[type_offset(e->x.p)])
1408 if (value_zero_p(e->d) && e->x.p->type == fractional)
1409 fract = &e->x.p->arr[0];
1410 else {
1411 for (e = val; value_zero_p(e->d) && e->x.p->type != polynomial;
1412 e = &e->x.p->arr[type_offset(e->x.p)])
1414 assert(value_zero_p(e->d) && e->x.p->type == polynomial);
1417 int offset = type_offset(e->x.p);
1419 assert(value_notzero_p(e->x.p->arr[offset+1].d));
1420 assert(value_notzero_p(e->x.p->arr[offset+1].x.n));
1421 if (value_neg_p(e->x.p->arr[offset+1].x.n)) {
1422 value_oppose(factor.d, e->x.p->arr[offset+1].x.n);
1423 value_assign(factor.x.n, e->x.p->arr[offset+1].d);
1424 } else {
1425 value_assign(factor.d, e->x.p->arr[offset+1].x.n);
1426 value_oppose(factor.x.n, e->x.p->arr[offset+1].d);
1429 free_evalue_refs(&e->x.p->arr[offset+1]);
1430 enode *p = e->x.p;
1431 value_clear(e->d);
1432 *e = e->x.p->arr[offset];
1434 emul(&factor, val);
1436 if (fract) {
1437 for (int i = 0; i < term.size(); ++i)
1438 term[i]->substitute(fract, val);
1440 free_evalue_refs(&p->arr[0]);
1441 } else {
1442 for (int i = 0; i < term.size(); ++i)
1443 term[i]->substitute(p->pos, val);
1446 free_evalue_refs(&factor);
1447 free_evalue_refs(val);
1448 delete val;
1450 free(p);
1453 order_sign partial_order::compare(const indicator_term *a, const indicator_term *b)
1455 unsigned dim = a->den.NumCols();
1456 order_sign sign = order_eq;
1457 EDomain *D = ind->D;
1458 unsigned MaxRays = ind->options->verify.barvinok->MaxRays;
1459 bool rational = a->sign == 0 || b->sign == 0;
1461 order_sign cached_sign = order_eq;
1462 for (int k = 0; k < dim; ++k) {
1463 cached_sign = evalue_diff_constant_sign(a->vertex[k], b->vertex[k]);
1464 if (cached_sign != order_eq)
1465 break;
1467 if (cached_sign != order_undefined)
1468 return cached_sign;
1470 order_cache_el cache_el;
1471 cached_sign = order_undefined;
1472 if (!rational)
1473 cached_sign = cache.check(a, b, cache_el);
1474 if (cached_sign != order_undefined) {
1475 cache_el.free();
1476 return cached_sign;
1479 if (rational && POL_ISSET(MaxRays, POL_INTEGER)) {
1480 ind->options->verify.barvinok->MaxRays &= ~POL_INTEGER;
1481 if (ind->options->verify.barvinok->MaxRays)
1482 ind->options->verify.barvinok->MaxRays |= POL_HIGH_BIT;
1485 sign = order_eq;
1487 vector<indicator_term *> term;
1489 for (int k = 0; k < dim; ++k) {
1490 /* compute a->vertex[k] - b->vertex[k] */
1491 evalue *diff;
1492 if (cache_el.e.size() <= k) {
1493 diff = ediff(a->vertex[k], b->vertex[k]);
1494 cache_el.e.push_back(diff);
1495 } else
1496 diff = cache_el.e[k];
1497 temp_evalue tdiff;
1498 if (term.size())
1499 tdiff = diff = ediff(term[0]->vertex[k], term[1]->vertex[k]);
1500 order_sign diff_sign;
1501 if (!D)
1502 diff_sign = order_undefined;
1503 else if (eequal(a->vertex[k], b->vertex[k]))
1504 diff_sign = order_eq;
1505 else
1506 diff_sign = evalue_sign(diff, D, ind->options->verify.barvinok);
1508 if (diff_sign == order_undefined) {
1509 assert(sign == order_le || sign == order_ge);
1510 if (sign == order_le)
1511 sign = order_lt;
1512 else
1513 sign = order_gt;
1514 break;
1516 if (diff_sign == order_lt) {
1517 if (sign == order_eq || sign == order_le)
1518 sign = order_lt;
1519 else
1520 sign = order_unknown;
1521 break;
1523 if (diff_sign == order_gt) {
1524 if (sign == order_eq || sign == order_ge)
1525 sign = order_gt;
1526 else
1527 sign = order_unknown;
1528 break;
1530 if (diff_sign == order_eq) {
1531 if (D == ind->D && term.size() == 0 && !rational &&
1532 !EVALUE_IS_ZERO(*diff))
1533 ind->add_substitution(diff);
1534 continue;
1536 if ((diff_sign == order_unknown) ||
1537 ((diff_sign == order_lt || diff_sign == order_le) && sign == order_ge) ||
1538 ((diff_sign == order_gt || diff_sign == order_ge) && sign == order_le)) {
1539 sign = order_unknown;
1540 break;
1543 sign = diff_sign;
1545 if (!term.size()) {
1546 term.push_back(new indicator_term(*a));
1547 term.push_back(new indicator_term(*b));
1549 substitute(term, diff);
1552 if (!rational)
1553 cache.add(cache_el, sign);
1554 else
1555 cache_el.free();
1557 if (D && D != ind->D)
1558 delete D;
1560 if (term.size()) {
1561 delete term[0];
1562 delete term[1];
1565 ind->options->verify.barvinok->MaxRays = MaxRays;
1566 return sign;
1569 bool partial_order::compared(const indicator_term* a, const indicator_term* b)
1571 map<const indicator_term *, vector<const indicator_term * > >::iterator j;
1573 j = lt.find(a);
1574 if (j != lt.end() && find(lt[a].begin(), lt[a].end(), b) != lt[a].end())
1575 return true;
1577 j = le.find(a);
1578 if (j != le.end() && find(le[a].begin(), le[a].end(), b) != le[a].end())
1579 return true;
1581 return false;
1584 void partial_order::add(const indicator_term* it,
1585 std::set<const indicator_term *> *filter)
1587 if (eq.find(it) != eq.end() && eq[it].size() == 1)
1588 return;
1590 typeof(head) head_copy(head);
1592 if (!filter)
1593 head.insert(it);
1595 std::set<const indicator_term *>::iterator i;
1596 for (i = head_copy.begin(); i != head_copy.end(); ++i) {
1597 if (*i == it)
1598 continue;
1599 if (eq.find(*i) != eq.end() && eq[*i].size() == 1)
1600 continue;
1601 if (filter) {
1602 if (filter->find(*i) == filter->end())
1603 continue;
1604 if (compared(*i, it))
1605 continue;
1607 order_sign sign = compare(it, *i);
1608 if (sign == order_lt) {
1609 lt[it].push_back(*i);
1610 inc_pred(*i);
1611 } else if (sign == order_le) {
1612 le[it].push_back(*i);
1613 inc_pred(*i);
1614 } else if (sign == order_eq) {
1615 set_equal(it, *i);
1616 return;
1617 } else if (sign == order_gt) {
1618 pending[*i].push_back(it);
1619 lt[*i].push_back(it);
1620 inc_pred(it);
1621 } else if (sign == order_ge) {
1622 pending[*i].push_back(it);
1623 le[*i].push_back(it);
1624 inc_pred(it);
1629 void partial_order::remove(const indicator_term* it)
1631 std::set<const indicator_term *> filter;
1632 map<const indicator_term *, vector<const indicator_term * > >::iterator i;
1634 assert(head.find(it) != head.end());
1636 i = eq.find(it);
1637 if (i != eq.end()) {
1638 assert(eq[it].size() >= 1);
1639 const indicator_term *base;
1640 if (eq[it].size() == 1) {
1641 base = eq[it][0];
1642 eq.erase(it);
1644 vector<const indicator_term * >::iterator j;
1645 j = find(eq[base].begin(), eq[base].end(), it);
1646 assert(j != eq[base].end());
1647 eq[base].erase(j);
1648 } else {
1649 /* "it" may no longer be the smallest, since the order
1650 * structure may have been copied from another one.
1652 sort(eq[it].begin()+1, eq[it].end(), pred.key_comp());
1653 assert(eq[it][0] == it);
1654 eq[it].erase(eq[it].begin());
1655 base = eq[it][0];
1656 eq[base] = eq[it];
1657 eq.erase(it);
1659 for (int j = 1; j < eq[base].size(); ++j)
1660 eq[eq[base][j]][0] = base;
1662 i = lt.find(it);
1663 if (i != lt.end()) {
1664 lt[base] = lt[it];
1665 lt.erase(it);
1668 i = le.find(it);
1669 if (i != le.end()) {
1670 le[base] = le[it];
1671 le.erase(it);
1674 i = pending.find(it);
1675 if (i != pending.end()) {
1676 pending[base] = pending[it];
1677 pending.erase(it);
1681 if (eq[base].size() == 1)
1682 eq.erase(base);
1684 head.erase(it);
1686 return;
1689 i = lt.find(it);
1690 if (i != lt.end()) {
1691 for (int j = 0; j < lt[it].size(); ++j) {
1692 filter.insert(lt[it][j]);
1693 dec_pred(lt[it][j]);
1695 lt.erase(it);
1698 i = le.find(it);
1699 if (i != le.end()) {
1700 for (int j = 0; j < le[it].size(); ++j) {
1701 filter.insert(le[it][j]);
1702 dec_pred(le[it][j]);
1704 le.erase(it);
1707 head.erase(it);
1709 i = pending.find(it);
1710 if (i != pending.end()) {
1711 vector<const indicator_term *> it_pending = pending[it];
1712 pending.erase(it);
1713 for (int j = 0; j < it_pending.size(); ++j) {
1714 filter.erase(it_pending[j]);
1715 add(it_pending[j], &filter);
1720 void partial_order::print(ostream& os, char **p)
1722 map<const indicator_term *, vector<const indicator_term * > >::iterator i;
1723 map<const indicator_term *, int >::iterator j;
1724 std::set<const indicator_term *>::iterator k;
1725 for (k = head.begin(); k != head.end(); ++k) {
1726 (*k)->print(os, p);
1727 os << endl;
1729 for (j = pred.begin(); j != pred.end(); ++j) {
1730 (*j).first->print(os, p);
1731 os << ": " << (*j).second << endl;
1733 for (i = lt.begin(); i != lt.end(); ++i) {
1734 (*i).first->print(os, p);
1735 assert(head.find((*i).first) != head.end() ||
1736 pred.find((*i).first) != pred.end());
1737 if (pred.find((*i).first) != pred.end())
1738 os << "(" << pred[(*i).first] << ")";
1739 os << " < ";
1740 for (int j = 0; j < (*i).second.size(); ++j) {
1741 if (j)
1742 os << ", ";
1743 (*i).second[j]->print(os, p);
1744 assert(pred.find((*i).second[j]) != pred.end());
1745 os << "(" << pred[(*i).second[j]] << ")";
1747 os << endl;
1749 for (i = le.begin(); i != le.end(); ++i) {
1750 (*i).first->print(os, p);
1751 assert(head.find((*i).first) != head.end() ||
1752 pred.find((*i).first) != pred.end());
1753 if (pred.find((*i).first) != pred.end())
1754 os << "(" << pred[(*i).first] << ")";
1755 os << " <= ";
1756 for (int j = 0; j < (*i).second.size(); ++j) {
1757 if (j)
1758 os << ", ";
1759 (*i).second[j]->print(os, p);
1760 assert(pred.find((*i).second[j]) != pred.end());
1761 os << "(" << pred[(*i).second[j]] << ")";
1763 os << endl;
1765 for (i = eq.begin(); i != eq.end(); ++i) {
1766 if ((*i).second.size() <= 1)
1767 continue;
1768 (*i).first->print(os, p);
1769 assert(head.find((*i).first) != head.end() ||
1770 pred.find((*i).first) != pred.end());
1771 if (pred.find((*i).first) != pred.end())
1772 os << "(" << pred[(*i).first] << ")";
1773 for (int j = 1; j < (*i).second.size(); ++j) {
1774 if (j)
1775 os << " = ";
1776 (*i).second[j]->print(os, p);
1777 assert(head.find((*i).second[j]) != head.end() ||
1778 pred.find((*i).second[j]) != pred.end());
1779 if (pred.find((*i).second[j]) != pred.end())
1780 os << "(" << pred[(*i).second[j]] << ")";
1782 os << endl;
1784 for (i = pending.begin(); i != pending.end(); ++i) {
1785 os << "pending on ";
1786 (*i).first->print(os, p);
1787 assert(head.find((*i).first) != head.end() ||
1788 pred.find((*i).first) != pred.end());
1789 if (pred.find((*i).first) != pred.end())
1790 os << "(" << pred[(*i).first] << ")";
1791 os << ": ";
1792 for (int j = 0; j < (*i).second.size(); ++j) {
1793 if (j)
1794 os << ", ";
1795 (*i).second[j]->print(os, p);
1796 assert(pred.find((*i).second[j]) != pred.end());
1797 os << "(" << pred[(*i).second[j]] << ")";
1799 os << endl;
1803 void indicator::add(const indicator_term* it)
1805 indicator_term *nt = new indicator_term(*it);
1806 if (options->reduce)
1807 nt->reduce_in_domain(P ? P : D->D);
1808 term.push_back(nt);
1809 order.add(nt, NULL);
1810 assert(term.size() == order.head.size() + order.pred.size());
1813 void indicator::remove(const indicator_term* it)
1815 vector<indicator_term *>::iterator i;
1816 i = find(term.begin(), term.end(), it);
1817 assert(i!= term.end());
1818 order.remove(it);
1819 term.erase(i);
1820 assert(term.size() == order.head.size() + order.pred.size());
1821 delete it;
1824 void indicator::expand_rational_vertex(const indicator_term *initial)
1826 int pos = initial->pos;
1827 remove(initial);
1828 if (ic.terms[pos].size() == 0) {
1829 Param_Vertices *V;
1830 FORALL_PVertex_in_ParamPolyhedron(V, PD, ic.PP) // _i is internal counter
1831 if (_i == pos) {
1832 ic.decompose_at_vertex(V, pos, options->verify.barvinok);
1833 break;
1835 END_FORALL_PVertex_in_ParamPolyhedron;
1837 for (int j = 0; j < ic.terms[pos].size(); ++j)
1838 add(ic.terms[pos][j]);
1841 void indicator::remove_initial_rational_vertices()
1843 do {
1844 const indicator_term *initial = NULL;
1845 std::set<const indicator_term *>::iterator i;
1846 for (i = order.head.begin(); i != order.head.end(); ++i) {
1847 if ((*i)->sign != 0)
1848 continue;
1849 if (order.eq.find(*i) != order.eq.end() && order.eq[*i].size() <= 1)
1850 continue;
1851 initial = *i;
1852 break;
1854 if (!initial)
1855 break;
1856 expand_rational_vertex(initial);
1857 } while(1);
1860 void indicator::reduce_in_domain(Polyhedron *D)
1862 for (int i = 0; i < term.size(); ++i)
1863 term[i]->reduce_in_domain(D);
1866 void indicator::print(ostream& os, char **p)
1868 assert(term.size() == order.head.size() + order.pred.size());
1869 for (int i = 0; i < term.size(); ++i) {
1870 term[i]->print(os, p);
1871 if (D->sample) {
1872 os << ": " << term[i]->eval(D->sample->p);
1874 os << endl;
1876 order.print(os, p);
1879 /* Remove pairs of opposite terms */
1880 void indicator::simplify()
1882 for (int i = 0; i < term.size(); ++i) {
1883 for (int j = i+1; j < term.size(); ++j) {
1884 if (term[i]->sign + term[j]->sign != 0)
1885 continue;
1886 if (term[i]->den != term[j]->den)
1887 continue;
1888 int k;
1889 for (k = 0; k < term[i]->den.NumCols(); ++k)
1890 if (!eequal(term[i]->vertex[k], term[j]->vertex[k]))
1891 break;
1892 if (k < term[i]->den.NumCols())
1893 continue;
1894 delete term[i];
1895 delete term[j];
1896 term.erase(term.begin()+j);
1897 term.erase(term.begin()+i);
1898 --i;
1899 break;
1904 void indicator::peel(int i, int j)
1906 if (j < i) {
1907 int tmp = i;
1908 i = j;
1909 j = tmp;
1911 assert (i < j);
1912 int dim = term[i]->den.NumCols();
1914 mat_ZZ common;
1915 mat_ZZ rest_i;
1916 mat_ZZ rest_j;
1917 int n_common = 0, n_i = 0, n_j = 0;
1919 common.SetDims(min(term[i]->den.NumRows(), term[j]->den.NumRows()), dim);
1920 rest_i.SetDims(term[i]->den.NumRows(), dim);
1921 rest_j.SetDims(term[j]->den.NumRows(), dim);
1923 int k, l;
1924 for (k = 0, l = 0; k < term[i]->den.NumRows() && l < term[j]->den.NumRows(); ) {
1925 int s = lex_cmp(term[i]->den[k], term[j]->den[l]);
1926 if (s == 0) {
1927 common[n_common++] = term[i]->den[k];
1928 ++k;
1929 ++l;
1930 } else if (s < 0)
1931 rest_i[n_i++] = term[i]->den[k++];
1932 else
1933 rest_j[n_j++] = term[j]->den[l++];
1935 while (k < term[i]->den.NumRows())
1936 rest_i[n_i++] = term[i]->den[k++];
1937 while (l < term[j]->den.NumRows())
1938 rest_j[n_j++] = term[j]->den[l++];
1939 common.SetDims(n_common, dim);
1940 rest_i.SetDims(n_i, dim);
1941 rest_j.SetDims(n_j, dim);
1943 for (k = 0; k <= n_i; ++k) {
1944 indicator_term *it = new indicator_term(*term[i]);
1945 it->den.SetDims(n_common + k, dim);
1946 for (l = 0; l < n_common; ++l)
1947 it->den[l] = common[l];
1948 for (l = 0; l < k; ++l)
1949 it->den[n_common+l] = rest_i[l];
1950 lex_order_rows(it->den);
1951 if (k)
1952 for (l = 0; l < dim; ++l)
1953 evalue_add_constant(it->vertex[l], rest_i[k-1][l]);
1954 term.push_back(it);
1957 for (k = 0; k <= n_j; ++k) {
1958 indicator_term *it = new indicator_term(*term[j]);
1959 it->den.SetDims(n_common + k, dim);
1960 for (l = 0; l < n_common; ++l)
1961 it->den[l] = common[l];
1962 for (l = 0; l < k; ++l)
1963 it->den[n_common+l] = rest_j[l];
1964 lex_order_rows(it->den);
1965 if (k)
1966 for (l = 0; l < dim; ++l)
1967 evalue_add_constant(it->vertex[l], rest_j[k-1][l]);
1968 term.push_back(it);
1970 term.erase(term.begin()+j);
1971 term.erase(term.begin()+i);
1974 void indicator::combine(const indicator_term *a, const indicator_term *b)
1976 int dim = a->den.NumCols();
1978 mat_ZZ common;
1979 mat_ZZ rest_i; /* factors in a, but not in b */
1980 mat_ZZ rest_j; /* factors in b, but not in a */
1981 int n_common = 0, n_i = 0, n_j = 0;
1983 common.SetDims(min(a->den.NumRows(), b->den.NumRows()), dim);
1984 rest_i.SetDims(a->den.NumRows(), dim);
1985 rest_j.SetDims(b->den.NumRows(), dim);
1987 int k, l;
1988 for (k = 0, l = 0; k < a->den.NumRows() && l < b->den.NumRows(); ) {
1989 int s = lex_cmp(a->den[k], b->den[l]);
1990 if (s == 0) {
1991 common[n_common++] = a->den[k];
1992 ++k;
1993 ++l;
1994 } else if (s < 0)
1995 rest_i[n_i++] = a->den[k++];
1996 else
1997 rest_j[n_j++] = b->den[l++];
1999 while (k < a->den.NumRows())
2000 rest_i[n_i++] = a->den[k++];
2001 while (l < b->den.NumRows())
2002 rest_j[n_j++] = b->den[l++];
2003 common.SetDims(n_common, dim);
2004 rest_i.SetDims(n_i, dim);
2005 rest_j.SetDims(n_j, dim);
2007 assert(order.eq[a].size() > 1);
2008 indicator_term *prev;
2010 prev = NULL;
2011 for (int k = n_i-1; k >= 0; --k) {
2012 indicator_term *it = new indicator_term(*b);
2013 it->den.SetDims(n_common + n_j + n_i-k, dim);
2014 for (int l = k; l < n_i; ++l)
2015 it->den[n_common+n_j+l-k] = rest_i[l];
2016 lex_order_rows(it->den);
2017 for (int m = 0; m < dim; ++m)
2018 evalue_add_constant(it->vertex[m], rest_i[k][m]);
2019 it->sign = -it->sign;
2020 if (prev) {
2021 order.pending[it].push_back(prev);
2022 order.lt[it].push_back(prev);
2023 order.inc_pred(prev);
2025 term.push_back(it);
2026 order.head.insert(it);
2027 prev = it;
2029 if (n_i) {
2030 indicator_term *it = new indicator_term(*b);
2031 it->den.SetDims(n_common + n_i + n_j, dim);
2032 for (l = 0; l < n_i; ++l)
2033 it->den[n_common+n_j+l] = rest_i[l];
2034 lex_order_rows(it->den);
2035 assert(prev);
2036 order.pending[a].push_back(prev);
2037 order.lt[a].push_back(prev);
2038 order.inc_pred(prev);
2039 order.replace(b, it);
2040 delete it;
2043 prev = NULL;
2044 for (int k = n_j-1; k >= 0; --k) {
2045 indicator_term *it = new indicator_term(*a);
2046 it->den.SetDims(n_common + n_i + n_j-k, dim);
2047 for (int l = k; l < n_j; ++l)
2048 it->den[n_common+n_i+l-k] = rest_j[l];
2049 lex_order_rows(it->den);
2050 for (int m = 0; m < dim; ++m)
2051 evalue_add_constant(it->vertex[m], rest_j[k][m]);
2052 it->sign = -it->sign;
2053 if (prev) {
2054 order.pending[it].push_back(prev);
2055 order.lt[it].push_back(prev);
2056 order.inc_pred(prev);
2058 term.push_back(it);
2059 order.head.insert(it);
2060 prev = it;
2062 if (n_j) {
2063 indicator_term *it = new indicator_term(*a);
2064 it->den.SetDims(n_common + n_i + n_j, dim);
2065 for (l = 0; l < n_j; ++l)
2066 it->den[n_common+n_i+l] = rest_j[l];
2067 lex_order_rows(it->den);
2068 assert(prev);
2069 order.pending[a].push_back(prev);
2070 order.lt[a].push_back(prev);
2071 order.inc_pred(prev);
2072 order.replace(a, it);
2073 delete it;
2076 assert(term.size() == order.head.size() + order.pred.size());
2079 bool indicator::handle_equal_numerators(const indicator_term *base)
2081 for (int i = 0; i < order.eq[base].size(); ++i) {
2082 for (int j = i+1; j < order.eq[base].size(); ++j) {
2083 if (order.eq[base][i]->is_opposite(order.eq[base][j])) {
2084 remove(order.eq[base][j]);
2085 remove(i ? order.eq[base][i] : base);
2086 return true;
2090 for (int j = 1; j < order.eq[base].size(); ++j)
2091 if (order.eq[base][j]->sign != base->sign) {
2092 combine(base, order.eq[base][j]);
2093 return true;
2095 return false;
2098 void indicator::substitute(evalue *equation)
2100 ::substitute(term, equation);
2103 void indicator::add_substitution(evalue *equation)
2105 for (int i = 0; i < substitutions.size(); ++i)
2106 if (eequal(substitutions[i], equation))
2107 return;
2108 evalue *copy = new evalue();
2109 value_init(copy->d);
2110 evalue_copy(copy, equation);
2111 substitutions.push_back(copy);
2114 void indicator::perform_pending_substitutions()
2116 if (substitutions.size() == 0)
2117 return;
2119 for (int i = 0; i < substitutions.size(); ++i) {
2120 substitute(substitutions[i]);
2121 free_evalue_refs(substitutions[i]);
2122 delete substitutions[i];
2124 substitutions.clear();
2125 order.resort();
2128 static void print_varlist(ostream& os, int n, char **names)
2130 int i;
2131 os << "[";
2132 for (i = 0; i < n; ++i) {
2133 if (i)
2134 os << ",";
2135 os << names[i];
2137 os << "]";
2140 void max_term::print(ostream& os, char **p, barvinok_options *options) const
2142 os << "{ ";
2143 print_varlist(os, domain->dimension(), p);
2144 os << " -> ";
2145 os << "[";
2146 for (int i = 0; i < max.size(); ++i) {
2147 if (i)
2148 os << ",";
2149 evalue_print(os, max[i], p);
2151 os << "]";
2152 os << " : ";
2153 domain->print_constraints(os, p, options);
2154 os << " }" << endl;
2157 /* T maps the compressed parameters to the original parameters,
2158 * while this max_term is based on the compressed parameters
2159 * and we want get the original parameters back.
2161 void max_term::substitute(Matrix *T, barvinok_options *options)
2163 assert(domain->dimension() == T->NbColumns-1);
2164 int nexist = domain->D->Dimension - (T->NbColumns-1);
2165 Matrix *Eq;
2166 Matrix *inv = left_inverse(T, &Eq);
2168 evalue denom;
2169 value_init(denom.d);
2170 value_init(denom.x.n);
2171 value_set_si(denom.x.n, 1);
2172 value_assign(denom.d, inv->p[inv->NbRows-1][inv->NbColumns-1]);
2174 vec_ZZ v;
2175 v.SetLength(inv->NbColumns);
2176 evalue* subs[inv->NbRows-1];
2177 for (int i = 0; i < inv->NbRows-1; ++i) {
2178 values2zz(inv->p[i], v, v.length());
2179 subs[i] = multi_monom(v);
2180 emul(&denom, subs[i]);
2182 free_evalue_refs(&denom);
2184 domain->substitute(subs, inv, Eq, options->MaxRays);
2185 Matrix_Free(Eq);
2187 for (int i = 0; i < max.size(); ++i) {
2188 evalue_substitute(max[i], subs);
2189 reduce_evalue(max[i]);
2192 for (int i = 0; i < inv->NbRows-1; ++i) {
2193 free_evalue_refs(subs[i]);
2194 delete subs[i];
2196 Matrix_Free(inv);
2199 int Last_Non_Zero(Value *p, unsigned len)
2201 for (int i = len-1; i >= 0; --i)
2202 if (value_notzero_p(p[i]))
2203 return i;
2204 return -1;
2207 Vector *max_term::eval(Value *val, unsigned MaxRays) const
2209 if (!domain->contains(val, domain->dimension()))
2210 return NULL;
2211 Vector *res = Vector_Alloc(max.size());
2212 for (int i = 0; i < max.size(); ++i) {
2213 compute_evalue(max[i], val, &res->p[i]);
2215 return res;
2218 struct split {
2219 evalue *constraint;
2220 enum sign { le, ge, lge } sign;
2222 split (evalue *c, enum sign s) : constraint(c), sign(s) {}
2225 static void split_on(const split& sp, EDomain *D,
2226 EDomain **Dlt, EDomain **Deq, EDomain **Dgt,
2227 lexmin_options *options)
2229 EDomain *ED[3];
2230 *Dlt = NULL;
2231 *Deq = NULL;
2232 *Dgt = NULL;
2233 ge_constraint *ge = D->compute_ge_constraint(sp.constraint);
2234 if (sp.sign == split::lge || sp.sign == split::ge)
2235 ED[2] = EDomain::new_from_ge_constraint(ge, 1, options->verify.barvinok);
2236 else
2237 ED[2] = NULL;
2238 if (sp.sign == split::lge || sp.sign == split::le)
2239 ED[0] = EDomain::new_from_ge_constraint(ge, -1, options->verify.barvinok);
2240 else
2241 ED[0] = NULL;
2243 assert(sp.sign == split::lge || sp.sign == split::ge || sp.sign == split::le);
2244 ED[1] = EDomain::new_from_ge_constraint(ge, 0, options->verify.barvinok);
2246 delete ge;
2248 for (int i = 0; i < 3; ++i) {
2249 if (!ED[i])
2250 continue;
2251 if (D->sample && ED[i]->contains(D->sample->p, D->sample->Size-1)) {
2252 ED[i]->sample = Vector_Alloc(D->sample->Size);
2253 Vector_Copy(D->sample->p, ED[i]->sample->p, D->sample->Size);
2254 } else if (emptyQ2(ED[i]->D) ||
2255 (options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2256 !(ED[i]->not_empty(options)))) {
2257 delete ED[i];
2258 ED[i] = NULL;
2261 *Dlt = ED[0];
2262 *Deq = ED[1];
2263 *Dgt = ED[2];
2266 ostream & operator<< (ostream & os, const vector<int> & v)
2268 os << "[";
2269 for (int i = 0; i < v.size(); ++i) {
2270 if (i)
2271 os << ", ";
2272 os << v[i];
2274 os << "]";
2275 return os;
2278 static bool isTranslation(Matrix *M)
2280 unsigned i, j;
2281 if (M->NbRows != M->NbColumns)
2282 return False;
2284 for (i = 0;i < M->NbRows; i ++)
2285 for (j = 0; j < M->NbColumns-1; j ++)
2286 if (i == j) {
2287 if(value_notone_p(M->p[i][j]))
2288 return False;
2289 } else {
2290 if(value_notzero_p(M->p[i][j]))
2291 return False;
2293 return value_one_p(M->p[M->NbRows-1][M->NbColumns-1]);
2296 static Matrix *compress_parameters(Polyhedron **P, Polyhedron **C,
2297 unsigned nparam, unsigned MaxRays)
2299 Matrix *M, *T, *CP;
2301 /* compress_parms doesn't like equalities that only involve parameters */
2302 for (int i = 0; i < (*P)->NbEq; ++i)
2303 assert(First_Non_Zero((*P)->Constraint[i]+1, (*P)->Dimension-nparam) != -1);
2305 M = Matrix_Alloc((*P)->NbEq, (*P)->Dimension+2);
2306 Vector_Copy((*P)->Constraint[0], M->p[0], (*P)->NbEq * ((*P)->Dimension+2));
2307 CP = compress_parms(M, nparam);
2308 Matrix_Free(M);
2310 if (isTranslation(CP)) {
2311 Matrix_Free(CP);
2312 return NULL;
2315 T = align_matrix(CP, (*P)->Dimension+1);
2316 *P = Polyhedron_Preimage(*P, T, MaxRays);
2317 Matrix_Free(T);
2319 *C = Polyhedron_Preimage(*C, CP, MaxRays);
2321 return CP;
2324 void construct_rational_vertices(Param_Polyhedron *PP, Matrix *T, unsigned dim,
2325 int nparam, vector<indicator_term *>& vertices)
2327 int i;
2328 Param_Vertices *PV;
2329 Value lcm, tmp;
2330 value_init(lcm);
2331 value_init(tmp);
2333 vec_ZZ v;
2334 v.SetLength(nparam+1);
2336 evalue factor;
2337 value_init(factor.d);
2338 value_init(factor.x.n);
2339 value_set_si(factor.x.n, 1);
2340 value_set_si(factor.d, 1);
2342 for (i = 0, PV = PP->V; PV; ++i, PV = PV->next) {
2343 indicator_term *term = new indicator_term(dim, i);
2344 vertices.push_back(term);
2345 Matrix *M = Matrix_Alloc(PV->Vertex->NbRows+nparam+1, nparam+1);
2346 value_set_si(lcm, 1);
2347 for (int j = 0; j < PV->Vertex->NbRows; ++j)
2348 value_lcm(lcm, lcm, PV->Vertex->p[j][nparam+1]);
2349 value_assign(M->p[M->NbRows-1][M->NbColumns-1], lcm);
2350 for (int j = 0; j < PV->Vertex->NbRows; ++j) {
2351 value_division(tmp, lcm, PV->Vertex->p[j][nparam+1]);
2352 Vector_Scale(PV->Vertex->p[j], M->p[j], tmp, nparam+1);
2354 for (int j = 0; j < nparam; ++j)
2355 value_assign(M->p[PV->Vertex->NbRows+j][j], lcm);
2356 if (T) {
2357 Matrix *M2 = Matrix_Alloc(T->NbRows, M->NbColumns);
2358 Matrix_Product(T, M, M2);
2359 Matrix_Free(M);
2360 M = M2;
2362 for (int j = 0; j < dim; ++j) {
2363 values2zz(M->p[j], v, nparam+1);
2364 term->vertex[j] = multi_monom(v);
2365 value_assign(factor.d, lcm);
2366 emul(&factor, term->vertex[j]);
2368 Matrix_Free(M);
2370 assert(i == PP->nbV);
2371 free_evalue_refs(&factor);
2372 value_clear(lcm);
2373 value_clear(tmp);
2376 static vector<max_term*> lexmin(indicator& ind, unsigned nparam,
2377 vector<int> loc)
2379 vector<max_term*> maxima;
2380 std::set<const indicator_term *>::iterator i;
2381 vector<int> best_score;
2382 vector<int> second_score;
2383 vector<int> neg_score;
2385 do {
2386 ind.perform_pending_substitutions();
2387 const indicator_term *best = NULL, *second = NULL, *neg = NULL,
2388 *neg_eq = NULL, *neg_le = NULL;
2389 for (i = ind.order.head.begin(); i != ind.order.head.end(); ++i) {
2390 vector<int> score;
2391 const indicator_term *term = *i;
2392 if (term->sign == 0) {
2393 ind.expand_rational_vertex(term);
2394 break;
2397 if (ind.order.eq.find(term) != ind.order.eq.end()) {
2398 int j;
2399 if (ind.order.eq[term].size() <= 1)
2400 continue;
2401 for (j = 1; j < ind.order.eq[term].size(); ++j)
2402 if (ind.order.pred.find(ind.order.eq[term][j]) !=
2403 ind.order.pred.end())
2404 break;
2405 if (j < ind.order.eq[term].size())
2406 continue;
2407 score.push_back(ind.order.eq[term].size());
2408 } else
2409 score.push_back(0);
2410 if (ind.order.le.find(term) != ind.order.le.end())
2411 score.push_back(ind.order.le[term].size());
2412 else
2413 score.push_back(0);
2414 if (ind.order.lt.find(term) != ind.order.lt.end())
2415 score.push_back(-ind.order.lt[term].size());
2416 else
2417 score.push_back(0);
2419 if (term->sign > 0) {
2420 if (!best || score < best_score) {
2421 second = best;
2422 second_score = best_score;
2423 best = term;
2424 best_score = score;
2425 } else if (!second || score < second_score) {
2426 second = term;
2427 second_score = score;
2429 } else {
2430 if (!neg_eq && ind.order.eq.find(term) != ind.order.eq.end()) {
2431 for (int j = 1; j < ind.order.eq[term].size(); ++j)
2432 if (ind.order.eq[term][j]->sign != term->sign) {
2433 neg_eq = term;
2434 break;
2437 if (!neg_le && ind.order.le.find(term) != ind.order.le.end())
2438 neg_le = term;
2439 if (!neg || score < neg_score) {
2440 neg = term;
2441 neg_score = score;
2445 if (i != ind.order.head.end())
2446 continue;
2448 if (!best && neg_eq) {
2449 assert(ind.order.eq[neg_eq].size() != 0);
2450 bool handled = ind.handle_equal_numerators(neg_eq);
2451 assert(handled);
2452 continue;
2455 if (!best && neg_le) {
2456 /* The smallest term is negative and <= some positive term */
2457 best = neg_le;
2458 neg = NULL;
2461 if (!best) {
2462 /* apparently there can be negative initial term on empty domains */
2463 if (ind.options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2464 ind.options->verify.barvinok->lp_solver == BV_LP_POLYLIB)
2465 assert(!neg);
2466 break;
2469 if (!second && !neg) {
2470 const indicator_term *rat = NULL;
2471 assert(best);
2472 if (ind.order.le.find(best) == ind.order.le.end()) {
2473 if (ind.order.eq.find(best) != ind.order.eq.end()) {
2474 bool handled = ind.handle_equal_numerators(best);
2475 if (ind.options->emptiness_check !=
2476 BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2477 ind.options->verify.barvinok->lp_solver == BV_LP_POLYLIB)
2478 assert(handled);
2479 /* If !handled then the leading coefficient is bigger than one;
2480 * must be an empty domain
2482 if (handled)
2483 continue;
2484 else
2485 break;
2487 maxima.push_back(ind.create_max_term(best));
2488 break;
2490 for (int j = 0; j < ind.order.le[best].size(); ++j) {
2491 if (ind.order.le[best][j]->sign == 0) {
2492 if (!rat && ind.order.pred[ind.order.le[best][j]] == 1)
2493 rat = ind.order.le[best][j];
2494 } else if (ind.order.le[best][j]->sign > 0) {
2495 second = ind.order.le[best][j];
2496 break;
2497 } else if (!neg)
2498 neg = ind.order.le[best][j];
2501 if (!second && !neg) {
2502 assert(rat);
2503 ind.order.unset_le(best, rat);
2504 ind.expand_rational_vertex(rat);
2505 continue;
2508 if (!second)
2509 second = neg;
2511 ind.order.unset_le(best, second);
2514 if (!second)
2515 second = neg;
2517 unsigned dim = best->den.NumCols();
2518 temp_evalue diff;
2519 order_sign sign;
2520 for (int k = 0; k < dim; ++k) {
2521 diff = ediff(best->vertex[k], second->vertex[k]);
2522 sign = evalue_sign(diff, ind.D, ind.options->verify.barvinok);
2524 /* neg can never be smaller than best, unless it may still cancel.
2525 * This can happen if positive terms have been determined to
2526 * be equal or less than or equal to some negative term.
2528 if (second == neg && !neg_eq && !neg_le) {
2529 if (sign == order_ge)
2530 sign = order_eq;
2531 if (sign == order_unknown)
2532 sign = order_le;
2535 if (sign != order_eq)
2536 break;
2537 if (!EVALUE_IS_ZERO(*diff)) {
2538 ind.add_substitution(diff);
2539 ind.perform_pending_substitutions();
2542 if (sign == order_eq) {
2543 ind.order.set_equal(best, second);
2544 continue;
2546 if (sign == order_lt) {
2547 ind.order.lt[best].push_back(second);
2548 ind.order.inc_pred(second);
2549 continue;
2551 if (sign == order_gt) {
2552 ind.order.lt[second].push_back(best);
2553 ind.order.inc_pred(best);
2554 continue;
2557 split sp(diff, sign == order_le ? split::le :
2558 sign == order_ge ? split::ge : split::lge);
2560 EDomain *Dlt, *Deq, *Dgt;
2561 split_on(sp, ind.D, &Dlt, &Deq, &Dgt, ind.options);
2562 if (ind.options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE)
2563 assert(Dlt || Deq || Dgt);
2564 else if (!(Dlt || Deq || Dgt))
2565 /* Must have been empty all along */
2566 break;
2568 if (Deq && (Dlt || Dgt)) {
2569 int locsize = loc.size();
2570 loc.push_back(0);
2571 indicator indeq(ind, Deq);
2572 Deq = NULL;
2573 indeq.add_substitution(diff);
2574 indeq.perform_pending_substitutions();
2575 vector<max_term*> maxeq = lexmin(indeq, nparam, loc);
2576 maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
2577 loc.resize(locsize);
2579 if (Dgt && Dlt) {
2580 int locsize = loc.size();
2581 loc.push_back(1);
2582 indicator indgt(ind, Dgt);
2583 Dgt = NULL;
2584 /* we don't know the new location of these terms in indgt */
2586 indgt.order.lt[second].push_back(best);
2587 indgt.order.inc_pred(best);
2589 vector<max_term*> maxgt = lexmin(indgt, nparam, loc);
2590 maxima.insert(maxima.end(), maxgt.begin(), maxgt.end());
2591 loc.resize(locsize);
2594 if (Deq) {
2595 loc.push_back(0);
2596 ind.set_domain(Deq);
2597 ind.add_substitution(diff);
2598 ind.perform_pending_substitutions();
2600 if (Dlt) {
2601 loc.push_back(-1);
2602 ind.set_domain(Dlt);
2603 ind.order.lt[best].push_back(second);
2604 ind.order.inc_pred(second);
2606 if (Dgt) {
2607 loc.push_back(1);
2608 ind.set_domain(Dgt);
2609 ind.order.lt[second].push_back(best);
2610 ind.order.inc_pred(best);
2612 } while(1);
2614 return maxima;
2617 static void lexmin_base(Polyhedron *P, Polyhedron *C,
2618 Matrix *CP, Matrix *T,
2619 vector<max_term*>& all_max,
2620 lexmin_options *options)
2622 unsigned nparam = C->Dimension;
2623 Param_Polyhedron *PP = NULL;
2625 PP = Polyhedron2Param_Polyhedron(P, C, options->verify.barvinok);
2627 unsigned dim = P->Dimension - nparam;
2629 int nd = -1;
2631 indicator_constructor ic(P, dim, PP, T);
2633 vector<indicator_term *> all_vertices;
2634 construct_rational_vertices(PP, T, T ? T->NbRows-nparam-1 : dim,
2635 nparam, all_vertices);
2637 Polyhedron *TC = true_context(P, C, options->verify.barvinok->MaxRays);
2638 FORALL_REDUCED_DOMAIN(PP, TC, nd, options->verify.barvinok, i, D, rVD)
2639 Param_Vertices *V;
2641 EDomain *epVD = new EDomain(rVD);
2642 indicator ind(ic, D, epVD, options);
2644 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
2645 ind.add(all_vertices[_i]);
2646 END_FORALL_PVertex_in_ParamPolyhedron;
2648 ind.remove_initial_rational_vertices();
2650 vector<int> loc;
2651 vector<max_term*> maxima = lexmin(ind, nparam, loc);
2652 if (CP)
2653 for (int j = 0; j < maxima.size(); ++j)
2654 maxima[j]->substitute(CP, options->verify.barvinok);
2655 all_max.insert(all_max.end(), maxima.begin(), maxima.end());
2657 Domain_Free(rVD);
2658 END_FORALL_REDUCED_DOMAIN
2659 Polyhedron_Free(TC);
2660 for (int i = 0; i < all_vertices.size(); ++i)
2661 delete all_vertices[i];
2662 Param_Polyhedron_Free(PP);
2665 static vector<max_term*> lexmin(Polyhedron *P, Polyhedron *C,
2666 lexmin_options *options)
2668 unsigned nparam = C->Dimension;
2669 Matrix *T = NULL, *CP = NULL;
2670 Polyhedron *Porig = P;
2671 Polyhedron *Corig = C;
2672 vector<max_term*> all_max;
2674 if (emptyQ2(P))
2675 return all_max;
2677 POL_ENSURE_VERTICES(P);
2679 if (emptyQ2(P))
2680 return all_max;
2682 assert(P->NbBid == 0);
2684 if (P->NbEq > 0)
2685 remove_all_equalities(&P, &C, &CP, &T, nparam,
2686 options->verify.barvinok->MaxRays);
2687 if (!emptyQ2(P))
2688 lexmin_base(P, C, CP, T, all_max, options);
2689 done:
2690 if (CP)
2691 Matrix_Free(CP);
2692 if (T)
2693 Matrix_Free(T);
2694 if (C != Corig)
2695 Polyhedron_Free(C);
2696 if (P != Porig)
2697 Polyhedron_Free(P);
2699 return all_max;
2702 static void verify_results(Polyhedron *A, Polyhedron *C,
2703 vector<max_term*>& maxima,
2704 struct verify_options *options);
2706 int main(int argc, char **argv)
2708 Polyhedron *A, *C;
2709 Matrix *MA;
2710 int bignum;
2711 char **iter_names, **param_names;
2712 int print_solution = 1;
2713 struct lexmin_options options;
2714 static struct argp_child argp_children[] = {
2715 { &barvinok_argp, 0, 0, 0 },
2716 { &verify_argp, 0, "verification", 1 },
2717 { 0 }
2719 static struct argp argp = { argp_options, parse_opt, 0, 0, argp_children };
2720 struct barvinok_options *bv_options;
2722 bv_options = barvinok_options_new_with_defaults();
2723 bv_options->lookup_table = 0;
2725 options.verify.barvinok = bv_options;
2726 set_program_name(argv[0]);
2727 argp_parse(&argp, argc, argv, 0, 0, &options);
2729 MA = Matrix_Read();
2730 C = Constraints2Polyhedron(MA, bv_options->MaxRays);
2731 Matrix_Free(MA);
2732 fscanf(stdin, " %d", &bignum);
2733 assert(bignum == -1);
2734 MA = Matrix_Read();
2735 A = Constraints2Polyhedron(MA, bv_options->MaxRays);
2736 Matrix_Free(MA);
2738 verify_options_set_range(&options.verify, A->Dimension);
2740 if (options.verify.verify)
2741 print_solution = 0;
2743 iter_names = util_generate_names(A->Dimension - C->Dimension, "i");
2744 param_names = util_generate_names(C->Dimension, "p");
2745 if (print_solution) {
2746 Polyhedron_Print(stdout, P_VALUE_FMT, A);
2747 Polyhedron_Print(stdout, P_VALUE_FMT, C);
2749 vector<max_term*> maxima = lexmin(A, C, &options);
2750 if (print_solution)
2751 for (int i = 0; i < maxima.size(); ++i)
2752 maxima[i]->print(cout, param_names, options.verify.barvinok);
2754 if (options.verify.verify)
2755 verify_results(A, C, maxima, &options.verify);
2757 for (int i = 0; i < maxima.size(); ++i)
2758 delete maxima[i];
2760 util_free_names(A->Dimension - C->Dimension, iter_names);
2761 util_free_names(C->Dimension, param_names);
2762 Polyhedron_Free(A);
2763 Polyhedron_Free(C);
2765 barvinok_options_free(bv_options);
2767 return 0;
2770 static bool lexmin(int pos, Polyhedron *P, Value *context)
2772 Value LB, UB, k;
2773 int lu_flags;
2774 bool found = false;
2776 if (emptyQ(P))
2777 return false;
2779 value_init(LB); value_init(UB); value_init(k);
2780 value_set_si(LB,0);
2781 value_set_si(UB,0);
2782 lu_flags = lower_upper_bounds(pos,P,context,&LB,&UB);
2783 assert(!(lu_flags & LB_INFINITY));
2785 value_set_si(context[pos],0);
2786 if (!lu_flags && value_lt(UB,LB)) {
2787 value_clear(LB); value_clear(UB); value_clear(k);
2788 return false;
2790 if (!P->next) {
2791 value_assign(context[pos], LB);
2792 value_clear(LB); value_clear(UB); value_clear(k);
2793 return true;
2795 for (value_assign(k,LB); lu_flags || value_le(k,UB); value_increment(k,k)) {
2796 value_assign(context[pos],k);
2797 if ((found = lexmin(pos+1, P->next, context)))
2798 break;
2800 if (!found)
2801 value_set_si(context[pos],0);
2802 value_clear(LB); value_clear(UB); value_clear(k);
2803 return found;
2806 static void print_list(FILE *out, Value *z, const char* brackets, int len)
2808 fprintf(out, "%c", brackets[0]);
2809 value_print(out, VALUE_FMT,z[0]);
2810 for (int k = 1; k < len; ++k) {
2811 fprintf(out, ", ");
2812 value_print(out, VALUE_FMT,z[k]);
2814 fprintf(out, "%c", brackets[1]);
2817 static int check_poly_lexmin(const struct check_poly_data *data,
2818 int nparam, Value *z,
2819 const struct verify_options *options);
2821 struct check_poly_lexmin_data : public check_poly_data {
2822 Polyhedron *S;
2823 vector<max_term*>& maxima;
2825 check_poly_lexmin_data(Polyhedron *S, Value *z,
2826 vector<max_term*>& maxima) : S(S), maxima(maxima) {
2827 this->z = z;
2828 this->check = check_poly_lexmin;
2832 static int check_poly_lexmin(const struct check_poly_data *data,
2833 int nparam, Value *z,
2834 const struct verify_options *options)
2836 const check_poly_lexmin_data *lexmin_data;
2837 lexmin_data = static_cast<const check_poly_lexmin_data *>(data);
2838 Polyhedron *S = lexmin_data->S;
2839 vector<max_term*>& maxima = lexmin_data->maxima;
2840 int k;
2841 bool found = lexmin(1, S, lexmin_data->z);
2843 if (options->print_all) {
2844 printf("lexmin");
2845 print_list(stdout, z, "()", nparam);
2846 printf(" = ");
2847 if (found)
2848 print_list(stdout, lexmin_data->z+1, "[]", S->Dimension-nparam);
2849 printf(" ");
2852 Vector *min = NULL;
2853 for (int i = 0; i < maxima.size(); ++i)
2854 if ((min = maxima[i]->eval(z, options->barvinok->MaxRays)))
2855 break;
2857 int ok = !(found ^ !!min);
2858 if (found && min)
2859 for (int i = 0; i < S->Dimension-nparam; ++i)
2860 if (value_ne(lexmin_data->z[1+i], min->p[i])) {
2861 ok = 0;
2862 break;
2864 if (!ok) {
2865 printf("\n");
2866 fflush(stdout);
2867 fprintf(stderr, "Error !\n");
2868 fprintf(stderr, "lexmin");
2869 print_list(stderr, z, "()", nparam);
2870 fprintf(stderr, " should be ");
2871 if (found)
2872 print_list(stderr, lexmin_data->z+1, "[]", S->Dimension-nparam);
2873 fprintf(stderr, " while digging gives ");
2874 if (min)
2875 print_list(stderr, min->p, "[]", S->Dimension-nparam);
2876 fprintf(stderr, ".\n");
2877 return 0;
2878 } else if (options->print_all)
2879 printf("OK.\n");
2880 if (min)
2881 Vector_Free(min);
2883 for (k = 1; k <= S->Dimension-nparam; ++k)
2884 value_set_si(lexmin_data->z[k], 0);
2887 void verify_results(Polyhedron *A, Polyhedron *C, vector<max_term*>& maxima,
2888 struct verify_options *options)
2890 Polyhedron *CS, *S;
2891 unsigned nparam = C->Dimension;
2892 unsigned MaxRays = options->barvinok->MaxRays;
2893 Vector *p;
2894 int i;
2895 int st;
2897 CS = check_poly_context_scan(A, &C, nparam, options);
2899 p = Vector_Alloc(A->Dimension+2);
2900 value_set_si(p->p[A->Dimension+1], 1);
2902 S = Polyhedron_Scan(A, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2904 check_poly_init(C, options);
2906 if (S) {
2907 if (!(CS && emptyQ2(CS))) {
2908 check_poly_lexmin_data data(S, p->p, maxima);
2909 check_poly(CS, &data, nparam, 0, p->p+S->Dimension-nparam+1, options);
2911 Domain_Free(S);
2914 if (!options->print_all)
2915 printf("\n");
2917 Vector_Free(p);
2918 if (CS) {
2919 Domain_Free(CS);
2920 Polyhedron_Free(C);