?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
counter.h: directly include required headers
[barvinok.git] / lexmin.cc
blobe735003dc6161530bf48e7cb312dc24e93064719
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 <isl_set_polylib.h>
13 #include <barvinok/barvinok.h>
14 #include <barvinok/evalue.h>
15 #include <barvinok/options.h>
16 #include <barvinok/util.h>
17 #include "conversion.h"
18 #include "decomposer.h"
19 #include "lattice_point.h"
20 #include "reduce_domain.h"
21 #include "mat_util.h"
22 #include "edomain.h"
23 #include "evalue_util.h"
24 #include "remove_equalities.h"
25 #include "polysign.h"
26 #include "verify.h"
27 #include "lexmin.h"
28 #include "param_util.h"
29
30 #undef CS /* for Solaris 10 */
31
32 using namespace NTL;
33
34 using std::vector;
35 using std::map;
36 using std::cerr;
37 using std::cout;
38 using std::endl;
39 using std::ostream;
40
41 #define ALLOC(type) (type*)malloc(sizeof(type))
42 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
43
44 static int type_offset(enode *p)
45 {
46 return p->type == fractional ? 1 :
47 p->type == flooring ? 1 : 0;
48 }
49
50 void compute_evalue(evalue *e, Value *val, Value *res)
51 {
52 double d = compute_evalue(e, val);
53 if (d > 0)
54 d += .25;
55 else
56 d -= .25;
57 value_set_double(*res, d);
58 }
59
60 struct indicator_term {
61 int sign;
62 int pos; /* number of rational vertex */
63 int n; /* number of cone associated to given rational vertex */
64 mat_ZZ den;
65 evalue **vertex;
66
67 indicator_term(unsigned dim, int pos) {
68 den.SetDims(0, dim);
69 vertex = new evalue* [dim];
70 this->pos = pos;
71 n = -1;
72 sign = 0;
73 }
74 indicator_term(unsigned dim, int pos, int n) {
75 den.SetDims(dim, dim);
76 vertex = new evalue* [dim];
77 this->pos = pos;
78 this->n = n;
79 }
80 indicator_term(const indicator_term& src) {
81 sign = src.sign;
82 pos = src.pos;
83 n = src.n;
84 den = src.den;
85 unsigned dim = den.NumCols();
86 vertex = new evalue* [dim];
87 for (int i = 0; i < dim; ++i) {
88 vertex[i] = ALLOC(evalue);
89 value_init(vertex[i]->d);
90 evalue_copy(vertex[i], src.vertex[i]);
91 }
92 }
93 void swap(indicator_term *other) {
94 int tmp;
95 tmp = sign; sign = other->sign; other->sign = tmp;
96 tmp = pos; pos = other->pos; other->pos = tmp;
97 tmp = n; n = other->n; other->n = tmp;
98 mat_ZZ tmp_den = den; den = other->den; other->den = tmp_den;
99 unsigned dim = den.NumCols();
100 for (int i = 0; i < dim; ++i) {
101 evalue *tmp = vertex[i];
102 vertex[i] = other->vertex[i];
103 other->vertex[i] = tmp;
104 }
105 }
106 ~indicator_term() {
107 unsigned dim = den.NumCols();
108 for (int i = 0; i < dim; ++i)
109 evalue_free(vertex[i]);
110 delete [] vertex;
111 }
112 void print(ostream& os, char **p) const;
113 void substitute(Matrix *T);
114 void normalize();
115 void substitute(evalue *fract, evalue *val);
116 void substitute(int pos, evalue *val);
117 void reduce_in_domain(Polyhedron *D);
118 bool is_opposite(const indicator_term *neg) const;
119 vec_ZZ eval(Value *val) const {
120 vec_ZZ v;
121 unsigned dim = den.NumCols();
122 v.SetLength(dim);
123 Value tmp;
124 value_init(tmp);
125 for (int i = 0; i < dim; ++i) {
126 compute_evalue(vertex[i], val, &tmp);
127 value2zz(tmp, v[i]);
128 }
129 value_clear(tmp);
130 return v;
131 }
132 };
133
134 static int evalue_rational_cmp(const evalue *e1, const evalue *e2)
135 {
136 int r;
137 Value m;
138 Value m2;
139 value_init(m);
140 value_init(m2);
141
142 assert(value_notzero_p(e1->d));
143 assert(value_notzero_p(e2->d));
144 value_multiply(m, e1->x.n, e2->d);
145 value_multiply(m2, e2->x.n, e1->d);
146 if (value_lt(m, m2))
147 r = -1;
148 else if (value_gt(m, m2))
149 r = 1;
150 else
151 r = 0;
152 value_clear(m);
153 value_clear(m2);
154
155 return r;
156 }
157
158 static int evalue_cmp(const evalue *e1, const evalue *e2)
159 {
160 if (value_notzero_p(e1->d)) {
161 if (value_zero_p(e2->d))
162 return -1;
163 return evalue_rational_cmp(e1, e2);
164 }
165 if (value_notzero_p(e2->d))
166 return 1;
167 if (e1->x.p->type != e2->x.p->type)
168 return e1->x.p->type - e2->x.p->type;
169 if (e1->x.p->size != e2->x.p->size)
170 return e1->x.p->size - e2->x.p->size;
171 if (e1->x.p->pos != e2->x.p->pos)
172 return e1->x.p->pos - e2->x.p->pos;
173 assert(e1->x.p->type == polynomial ||
174 e1->x.p->type == fractional ||
175 e1->x.p->type == flooring);
176 for (int i = 0; i < e1->x.p->size; ++i) {
177 int s = evalue_cmp(&e1->x.p->arr[i], &e2->x.p->arr[i]);
178 if (s)
179 return s;
180 }
181 return 0;
182 }
183
184 void evalue_length(evalue *e, int len[2])
185 {
186 len[0] = 0;
187 len[1] = 0;
188
189 while (value_zero_p(e->d)) {
190 assert(e->x.p->type == polynomial ||
191 e->x.p->type == fractional ||
192 e->x.p->type == flooring);
193 if (e->x.p->type == polynomial)
194 len[1]++;
195 else
196 len[0]++;
197 int offset = type_offset(e->x.p);
198 assert(e->x.p->size == offset+2);
199 e = &e->x.p->arr[offset];
200 }
201 }
202
203 static bool it_smaller(const indicator_term* it1, const indicator_term* it2)
204 {
205 if (it1 == it2)
206 return false;
207 int len1[2], len2[2];
208 unsigned dim = it1->den.NumCols();
209 for (int i = 0; i < dim; ++i) {
210 evalue_length(it1->vertex[i], len1);
211 evalue_length(it2->vertex[i], len2);
212 if (len1[0] != len2[0])
213 return len1[0] < len2[0];
214 if (len1[1] != len2[1])
215 return len1[1] < len2[1];
216 }
217 if (it1->pos != it2->pos)
218 return it1->pos < it2->pos;
219 if (it1->n != it2->n)
220 return it1->n < it2->n;
221 int s = lex_cmp(it1->den, it2->den);
222 if (s)
223 return s < 0;
224 for (int i = 0; i < dim; ++i) {
225 s = evalue_cmp(it1->vertex[i], it2->vertex[i]);
226 if (s)
227 return s < 0;
228 }
229 assert(it1->sign != 0);
230 assert(it2->sign != 0);
231 if (it1->sign != it2->sign)
232 return it1->sign > 0;
233 return it1 < it2;
234 }
235
236 struct smaller_it {
237 static const int requires_resort;
238 bool operator()(const indicator_term* it1, const indicator_term* it2) const {
239 return it_smaller(it1, it2);
240 }
241 };
242 const int smaller_it::requires_resort = 1;
243
244 struct smaller_it_p {
245 static const int requires_resort;
246 bool operator()(const indicator_term* it1, const indicator_term* it2) const {
247 return it1 < it2;
248 }
249 };
250 const int smaller_it_p::requires_resort = 0;
251
252 /* Returns true if this and neg are opposite using the knowledge
253 * that they have the same numerator.
254 * In particular, we check that the signs are different and that
255 * the denominator is the same.
256 */
257 bool indicator_term::is_opposite(const indicator_term *neg) const
258 {
259 if (sign + neg->sign != 0)
260 return false;
261 if (den != neg->den)
262 return false;
263 return true;
264 }
265
266 void indicator_term::reduce_in_domain(Polyhedron *D)
267 {
268 for (int k = 0; k < den.NumCols(); ++k) {
269 reduce_evalue_in_domain(vertex[k], D);
270 if (evalue_range_reduction_in_domain(vertex[k], D))
271 reduce_evalue(vertex[k]);
272 }
273 }
274
275 void indicator_term::print(ostream& os, char **p) const
276 {
277 unsigned dim = den.NumCols();
278 unsigned factors = den.NumRows();
279 if (sign == 0)
280 os << " s ";
281 else if (sign > 0)
282 os << " + ";
283 else
284 os << " - ";
285 os << '[';
286 for (int i = 0; i < dim; ++i) {
287 if (i)
288 os << ',';
289 evalue_print(os, vertex[i], p);
290 }
291 os << ']';
292 for (int i = 0; i < factors; ++i) {
293 os << " + t" << i << "*[";
294 for (int j = 0; j < dim; ++j) {
295 if (j)
296 os << ',';
297 os << den[i][j];
298 }
299 os << "]";
300 }
301 os << " ((" << pos << ", " << n << ", " << (void*)this << "))";
302 }
303
304 /* Perform the substitution specified by T on the variables.
305 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
306 * The substitution is performed as in gen_fun::substitute
307 */
308 void indicator_term::substitute(Matrix *T)
309 {
310 unsigned dim = den.NumCols();
311 unsigned nparam = T->NbColumns - dim - 1;
312 unsigned newdim = T->NbRows - nparam - 1;
313 evalue **newvertex;
314 mat_ZZ trans;
315 matrix2zz(T, trans, newdim, dim);
316 trans = transpose(trans);
317 den *= trans;
318 newvertex = new evalue* [newdim];
319
320 vec_ZZ v;
321 v.SetLength(nparam+1);
322
323 evalue factor;
324 value_init(factor.d);
325 value_set_si(factor.d, 1);
326 value_init(factor.x.n);
327 for (int i = 0; i < newdim; ++i) {
328 values2zz(T->p[i]+dim, v, nparam+1);
329 newvertex[i] = multi_monom(v);
330
331 for (int j = 0; j < dim; ++j) {
332 if (value_zero_p(T->p[i][j]))
333 continue;
334 evalue term;
335 value_init(term.d);
336 evalue_copy(&term, vertex[j]);
337 value_assign(factor.x.n, T->p[i][j]);
338 emul(&factor, &term);
339 eadd(&term, newvertex[i]);
340 free_evalue_refs(&term);
341 }
342 }
343 free_evalue_refs(&factor);
344 for (int i = 0; i < dim; ++i)
345 evalue_free(vertex[i]);
346 delete [] vertex;
347 vertex = newvertex;
348 }
349
350 static void evalue_add_constant(evalue *e, ZZ v)
351 {
352 Value tmp;
353 value_init(tmp);
354
355 /* go down to constant term */
356 while (value_zero_p(e->d))
357 e = &e->x.p->arr[type_offset(e->x.p)];
358 /* and add v */
359 zz2value(v, tmp);
360 value_multiply(tmp, tmp, e->d);
361 value_addto(e->x.n, e->x.n, tmp);
362
363 value_clear(tmp);
364 }
365
366 /* Make all powers in denominator lexico-positive */
367 void indicator_term::normalize()
368 {
369 vec_ZZ extra_vertex;
370 extra_vertex.SetLength(den.NumCols());
371 for (int r = 0; r < den.NumRows(); ++r) {
372 for (int k = 0; k < den.NumCols(); ++k) {
373 if (den[r][k] == 0)
374 continue;
375 if (den[r][k] > 0)
376 break;
377 sign = -sign;
378 den[r] = -den[r];
379 extra_vertex += den[r];
380 break;
381 }
382 }
383 for (int k = 0; k < extra_vertex.length(); ++k)
384 if (extra_vertex[k] != 0)
385 evalue_add_constant(vertex[k], extra_vertex[k]);
386 }
387
388 static void substitute(evalue *e, evalue *fract, evalue *val)
389 {
390 evalue *t;
391
392 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
393 if (t->x.p->type == fractional && eequal(&t->x.p->arr[0], fract))
394 break;
395 }
396 if (value_notzero_p(t->d))
397 return;
398
399 free_evalue_refs(&t->x.p->arr[0]);
400 evalue *term = &t->x.p->arr[2];
401 enode *p = t->x.p;
402 value_clear(t->d);
403 *t = t->x.p->arr[1];
404
405 emul(val, term);
406 eadd(term, e);
407 free_evalue_refs(term);
408 free(p);
409
410 reduce_evalue(e);
411 }
412
413 void indicator_term::substitute(evalue *fract, evalue *val)
414 {
415 unsigned dim = den.NumCols();
416 for (int i = 0; i < dim; ++i) {
417 ::substitute(vertex[i], fract, val);
418 }
419 }
420
421 static void substitute(evalue *e, int pos, evalue *val)
422 {
423 evalue *t;
424
425 for (t = e; value_zero_p(t->d); t = &t->x.p->arr[type_offset(t->x.p)]) {
426 if (t->x.p->type == polynomial && t->x.p->pos == pos)
427 break;
428 }
429 if (value_notzero_p(t->d))
430 return;
431
432 evalue *term = &t->x.p->arr[1];
433 enode *p = t->x.p;
434 value_clear(t->d);
435 *t = t->x.p->arr[0];
436
437 emul(val, term);
438 eadd(term, e);
439 free_evalue_refs(term);
440 free(p);
441
442 reduce_evalue(e);
443 }
444
445 void indicator_term::substitute(int pos, evalue *val)
446 {
447 unsigned dim = den.NumCols();
448 for (int i = 0; i < dim; ++i) {
449 ::substitute(vertex[i], pos, val);
450 }
451 }
452
453 struct indicator_constructor : public signed_cone_consumer,
454 public vertex_decomposer {
455 vec_ZZ vertex;
456 vector<indicator_term*> *terms;
457 Matrix *T; /* Transformation to original space */
458 int nbV;
459 int pos;
460 int n;
461
462 indicator_constructor(unsigned dim, Param_Polyhedron *PP, Matrix *T) :
463 vertex_decomposer(PP, *this), T(T), nbV(PP->nbV) {
464 vertex.SetLength(dim);
465 terms = new vector<indicator_term*>[PP->nbV];
466 }
467 ~indicator_constructor() {
468 for (int i = 0; i < nbV; ++i)
469 for (int j = 0; j < terms[i].size(); ++j)
470 delete terms[i][j];
471 delete [] terms;
472 }
473 void print(ostream& os, char **p);
474
475 virtual void handle(const signed_cone& sc, barvinok_options *options);
476 void decompose_at_vertex(Param_Vertices *V, int _i,
477 barvinok_options *options) {
478 pos = _i;
479 n = 0;
480 vertex_decomposer::decompose_at_vertex(V, _i, options);
481 }
482 };
483
484 void indicator_constructor::handle(const signed_cone& sc, barvinok_options *options)
485 {
486 assert(sc.det == 1);
487 unsigned dim = vertex.length();
488
489 assert(sc.rays.NumRows() == dim);
490
491 indicator_term *term = new indicator_term(dim, pos, n++);
492 term->sign = sc.sign;
493 terms[vert].push_back(term);
494
495 lattice_point(V, sc.rays, vertex, term->vertex, options);
496
497 term->den = sc.rays;
498 for (int r = 0; r < dim; ++r) {
499 for (int j = 0; j < dim; ++j) {
500 if (term->den[r][j] == 0)
501 continue;
502 if (term->den[r][j] > 0)
503 break;
504 term->sign = -term->sign;
505 term->den[r] = -term->den[r];
506 vertex += term->den[r];
507 break;
508 }
509 }
510
511 for (int i = 0; i < dim; ++i) {
512 if (!term->vertex[i]) {
513 term->vertex[i] = ALLOC(evalue);
514 value_init(term->vertex[i]->d);
515 value_init(term->vertex[i]->x.n);
516 zz2value(vertex[i], term->vertex[i]->x.n);
517 value_set_si(term->vertex[i]->d, 1);
518 continue;
519 }
520 if (vertex[i] == 0)
521 continue;
522 evalue_add_constant(term->vertex[i], vertex[i]);
523 }
524
525 if (T) {
526 term->substitute(T);
527 term->normalize();
528 }
529
530 lex_order_rows(term->den);
531 }
532
533 void indicator_constructor::print(ostream& os, char **p)
534 {
535 for (int i = 0; i < PP->nbV; ++i)
536 for (int j = 0; j < terms[i].size(); ++j) {
537 os << "i: " << i << ", j: " << j << endl;
538 terms[i][j]->print(os, p);
539 os << endl;
540 }
541 }
542
543 struct order_cache_el {
544 vector<evalue *> e;
545 order_cache_el copy() const {
546 order_cache_el n;
547 for (int i = 0; i < e.size(); ++i) {
548 evalue *c = new evalue;
549 value_init(c->d);
550 evalue_copy(c, e[i]);
551 n.e.push_back(c);
552 }
553 return n;
554 }
555 void free() {
556 for (int i = 0; i < e.size(); ++i) {
557 free_evalue_refs(e[i]);
558 delete e[i];
559 }
560 }
561 void negate() {
562 evalue mone;
563 value_init(mone.d);
564 evalue_set_si(&mone, -1, 1);
565 for (int i = 0; i < e.size(); ++i)
566 emul(&mone, e[i]);
567 free_evalue_refs(&mone);
568 }
569 void print(ostream& os, char **p);
570 };
571
572 void order_cache_el::print(ostream& os, char **p)
573 {
574 os << "[";
575 for (int i = 0; i < e.size(); ++i) {
576 if (i)
577 os << ",";
578 evalue_print(os, e[i], p);
579 }
580 os << "]";
581 }
582
583 struct order_cache {
584 vector<order_cache_el> lt;
585 vector<order_cache_el> le;
586 vector<order_cache_el> unknown;
587
588 void clear_transients() {
589 for (int i = 0; i < le.size(); ++i)
590 le[i].free();
591 for (int i = 0; i < unknown.size(); ++i)
592 unknown[i].free();
593 le.clear();
594 unknown.clear();
595 }
596 ~order_cache() {
597 clear_transients();
598 for (int i = 0; i < lt.size(); ++i)
599 lt[i].free();
600 lt.clear();
601 }
602 void add(order_cache_el& cache_el, order_sign sign);
603 order_sign check_lt(vector<order_cache_el>* list,
604 const indicator_term *a, const indicator_term *b,
605 order_cache_el& cache_el);
606 order_sign check_lt(const indicator_term *a, const indicator_term *b,
607 order_cache_el& cache_el);
608 order_sign check_direct(const indicator_term *a, const indicator_term *b,
609 order_cache_el& cache_el);
610 order_sign check(const indicator_term *a, const indicator_term *b,
611 order_cache_el& cache_el);
612 void copy(const order_cache& cache);
613 void print(ostream& os, char **p);
614 };
615
616 void order_cache::copy(const order_cache& cache)
617 {
618 for (int i = 0; i < cache.lt.size(); ++i) {
619 order_cache_el n = cache.lt[i].copy();
620 add(n, order_lt);
621 }
622 }
623
624 void order_cache::add(order_cache_el& cache_el, order_sign sign)
625 {
626 if (sign == order_lt) {
627 lt.push_back(cache_el);
628 } else if (sign == order_gt) {
629 cache_el.negate();
630 lt.push_back(cache_el);
631 } else if (sign == order_le) {
632 le.push_back(cache_el);
633 } else if (sign == order_ge) {
634 cache_el.negate();
635 le.push_back(cache_el);
636 } else if (sign == order_unknown) {
637 unknown.push_back(cache_el);
638 } else {
639 assert(sign == order_eq);
640 cache_el.free();
641 }
642 return;
643 }
644
645 /* compute a - b */
646 static evalue *ediff(const evalue *a, const evalue *b)
647 {
648 evalue mone;
649 value_init(mone.d);
650 evalue_set_si(&mone, -1, 1);
651 evalue *diff = new evalue;
652 value_init(diff->d);
653 evalue_copy(diff, b);
654 emul(&mone, diff);
655 eadd(a, diff);
656 reduce_evalue(diff);
657 free_evalue_refs(&mone);
658 return diff;
659 }
660
661 static bool evalue_first_difference(const evalue *e1, const evalue *e2,
662 const evalue **d1, const evalue **d2)
663 {
664 *d1 = e1;
665 *d2 = e2;
666
667 if (value_ne(e1->d, e2->d))
668 return true;
669
670 if (value_notzero_p(e1->d)) {
671 if (value_eq(e1->x.n, e2->x.n))
672 return false;
673 return true;
674 }
675 if (e1->x.p->type != e2->x.p->type)
676 return true;
677 if (e1->x.p->size != e2->x.p->size)
678 return true;
679 if (e1->x.p->pos != e2->x.p->pos)
680 return true;
681
682 assert(e1->x.p->type == polynomial ||
683 e1->x.p->type == fractional ||
684 e1->x.p->type == flooring);
685 int offset = type_offset(e1->x.p);
686 assert(e1->x.p->size == offset+2);
687 for (int i = 0; i < e1->x.p->size; ++i)
688 if (i != type_offset(e1->x.p) &&
689 !eequal(&e1->x.p->arr[i], &e2->x.p->arr[i]))
690 return true;
691
692 return evalue_first_difference(&e1->x.p->arr[offset],
693 &e2->x.p->arr[offset], d1, d2);
694 }
695
696 static order_sign evalue_diff_constant_sign(const evalue *e1, const evalue *e2)
697 {
698 if (!evalue_first_difference(e1, e2, &e1, &e2))
699 return order_eq;
700 if (value_zero_p(e1->d) || value_zero_p(e2->d))
701 return order_undefined;
702 int s = evalue_rational_cmp(e1, e2);
703 if (s < 0)
704 return order_lt;
705 else if (s > 0)
706 return order_gt;
707 else
708 return order_eq;
709 }
710
711 order_sign order_cache::check_lt(vector<order_cache_el>* list,
712 const indicator_term *a, const indicator_term *b,
713 order_cache_el& cache_el)
714 {
715 order_sign sign = order_undefined;
716 for (int i = 0; i < list->size(); ++i) {
717 int j;
718 for (j = cache_el.e.size(); j < (*list)[i].e.size(); ++j)
719 cache_el.e.push_back(ediff(a->vertex[j], b->vertex[j]));
720 for (j = 0; j < (*list)[i].e.size(); ++j) {
721 order_sign diff_sign;
722 diff_sign = evalue_diff_constant_sign((*list)[i].e[j], cache_el.e[j]);
723 if (diff_sign == order_gt) {
724 sign = order_lt;
725 break;
726 } else if (diff_sign == order_lt)
727 break;
728 else if (diff_sign == order_undefined)
729 break;
730 else
731 assert(diff_sign == order_eq);
732 }
733 if (j == (*list)[i].e.size())
734 sign = list == &lt ? order_lt : order_le;
735 if (sign != order_undefined)
736 break;
737 }
738 return sign;
739 }
740
741 order_sign order_cache::check_direct(const indicator_term *a,
742 const indicator_term *b,
743 order_cache_el& cache_el)
744 {
745 order_sign sign = check_lt(&lt, a, b, cache_el);
746 if (sign != order_undefined)
747 return sign;
748 sign = check_lt(&le, a, b, cache_el);
749 if (sign != order_undefined)
750 return sign;
751
752 for (int i = 0; i < unknown.size(); ++i) {
753 int j;
754 for (j = cache_el.e.size(); j < unknown[i].e.size(); ++j)
755 cache_el.e.push_back(ediff(a->vertex[j], b->vertex[j]));
756 for (j = 0; j < unknown[i].e.size(); ++j) {
757 if (!eequal(unknown[i].e[j], cache_el.e[j]))
758 break;
759 }
760 if (j == unknown[i].e.size()) {
761 sign = order_unknown;
762 break;
763 }
764 }
765 return sign;
766 }
767
768 order_sign order_cache::check(const indicator_term *a, const indicator_term *b,
769 order_cache_el& cache_el)
770 {
771 order_sign sign = check_direct(a, b, cache_el);
772 if (sign != order_undefined)
773 return sign;
774 int size = cache_el.e.size();
775 cache_el.negate();
776 sign = check_direct(a, b, cache_el);
777 cache_el.negate();
778 assert(cache_el.e.size() == size);
779 if (sign == order_undefined)
780 return sign;
781 if (sign == order_lt)
782 sign = order_gt;
783 else if (sign == order_le)
784 sign = order_ge;
785 else
786 assert(sign == order_unknown);
787 return sign;
788 }
789
790 struct indicator;
791
792 struct partial_order {
793 indicator *ind;
794
795 typedef std::set<const indicator_term *, smaller_it > head_type;
796 head_type head;
797 typedef map<const indicator_term *, int, smaller_it > pred_type;
798 pred_type pred;
799 typedef map<const indicator_term *, vector<const indicator_term * >, smaller_it > order_type;
800 order_type lt;
801 order_type le;
802 order_type eq;
803
804 order_type pending;
805
806 order_cache cache;
807
808 partial_order(indicator *ind) : ind(ind) {}
809 void copy(const partial_order& order,
810 map< const indicator_term *, indicator_term * > old2new);
811 void resort() {
812 order_type::iterator i;
813 pred_type::iterator j;
814 head_type::iterator k;
815
816 if (head.key_comp().requires_resort) {
817 head_type new_head;
818 for (k = head.begin(); k != head.end(); ++k)
819 new_head.insert(*k);
820 head.swap(new_head);
821 new_head.clear();
822 }
823
824 if (pred.key_comp().requires_resort) {
825 pred_type new_pred;
826 for (j = pred.begin(); j != pred.end(); ++j)
827 new_pred[(*j).first] = (*j).second;
828 pred.swap(new_pred);
829 new_pred.clear();
830 }
831
832 if (lt.key_comp().requires_resort) {
833 order_type m;
834 for (i = lt.begin(); i != lt.end(); ++i)
835 m[(*i).first] = (*i).second;
836 lt.swap(m);
837 m.clear();
838 }
839
840 if (le.key_comp().requires_resort) {
841 order_type m;
842 for (i = le.begin(); i != le.end(); ++i)
843 m[(*i).first] = (*i).second;
844 le.swap(m);
845 m.clear();
846 }
847
848 if (eq.key_comp().requires_resort) {
849 order_type m;
850 for (i = eq.begin(); i != eq.end(); ++i)
851 m[(*i).first] = (*i).second;
852 eq.swap(m);
853 m.clear();
854 }
855
856 if (pending.key_comp().requires_resort) {
857 order_type m;
858 for (i = pending.begin(); i != pending.end(); ++i)
859 m[(*i).first] = (*i).second;
860 pending.swap(m);
861 m.clear();
862 }
863 }
864
865 order_sign compare(const indicator_term *a, const indicator_term *b);
866 void set_equal(const indicator_term *a, const indicator_term *b);
867 void unset_le(const indicator_term *a, const indicator_term *b);
868 void dec_pred(const indicator_term *it) {
869 if (--pred[it] == 0) {
870 pred.erase(it);
871 head.insert(it);
872 }
873 }
874 void inc_pred(const indicator_term *it) {
875 if (head.find(it) != head.end())
876 head.erase(it);
877 pred[it]++;
878 }
879
880 bool compared(const indicator_term* a, const indicator_term* b);
881 void add(const indicator_term* it, std::set<const indicator_term *> *filter);
882 void remove(const indicator_term* it);
883
884 void print(ostream& os, char **p);
885
886 /* replace references to orig to references to replacement */
887 void replace(const indicator_term* orig, indicator_term* replacement);
888 void sanity_check() const;
889 };
890
891 /* We actually replace the contents of orig by that of replacement,
892 * but we have to be careful since replacing the content changes
893 * the order of orig in the maps.
894 */
895 void partial_order::replace(const indicator_term* orig, indicator_term* replacement)
896 {
897 head_type::iterator k;
898 k = head.find(orig);
899 bool is_head = k != head.end();
900 int orig_pred;
901 if (is_head) {
902 head.erase(orig);
903 } else {
904 orig_pred = pred[orig];
905 pred.erase(orig);
906 }
907 vector<const indicator_term * > orig_lt;
908 vector<const indicator_term * > orig_le;
909 vector<const indicator_term * > orig_eq;
910 vector<const indicator_term * > orig_pending;
911 order_type::iterator i;
912 bool in_lt = ((i = lt.find(orig)) != lt.end());
913 if (in_lt) {
914 orig_lt = (*i).second;
915 lt.erase(orig);
916 }
917 bool in_le = ((i = le.find(orig)) != le.end());
918 if (in_le) {
919 orig_le = (*i).second;
920 le.erase(orig);
921 }
922 bool in_eq = ((i = eq.find(orig)) != eq.end());
923 if (in_eq) {
924 orig_eq = (*i).second;
925 eq.erase(orig);
926 }
927 bool in_pending = ((i = pending.find(orig)) != pending.end());
928 if (in_pending) {
929 orig_pending = (*i).second;
930 pending.erase(orig);
931 }
932 indicator_term *old = const_cast<indicator_term *>(orig);
933 old->swap(replacement);
934 if (is_head)
935 head.insert(old);
936 else
937 pred[old] = orig_pred;
938 if (in_lt)
939 lt[old] = orig_lt;
940 if (in_le)
941 le[old] = orig_le;
942 if (in_eq)
943 eq[old] = orig_eq;
944 if (in_pending)
945 pending[old] = orig_pending;
946 }
947
948 void partial_order::unset_le(const indicator_term *a, const indicator_term *b)
949 {
950 vector<const indicator_term *>::iterator i;
951 i = std::find(le[a].begin(), le[a].end(), b);
952 le[a].erase(i);
953 if (le[a].size() == 0)
954 le.erase(a);
955 dec_pred(b);
956 i = std::find(pending[a].begin(), pending[a].end(), b);
957 if (i != pending[a].end())
958 pending[a].erase(i);
959 }
960
961 void partial_order::set_equal(const indicator_term *a, const indicator_term *b)
962 {
963 if (eq[a].size() == 0)
964 eq[a].push_back(a);
965 if (eq[b].size() == 0)
966 eq[b].push_back(b);
967 a = eq[a][0];
968 b = eq[b][0];
969 assert(a != b);
970 if (pred.key_comp()(b, a)) {
971 const indicator_term *c = a;
972 a = b;
973 b = c;
974 }
975
976 const indicator_term *base = a;
977
978 order_type::iterator i;
979
980 for (int j = 0; j < eq[b].size(); ++j) {
981 eq[base].push_back(eq[b][j]);
982 eq[eq[b][j]][0] = base;
983 }
984 eq[b].resize(1);
985
986 i = lt.find(b);
987 if (i != lt.end()) {
988 for (int j = 0; j < lt[b].size(); ++j) {
989 if (std::find(eq[base].begin(), eq[base].end(), lt[b][j]) != eq[base].end())
990 dec_pred(lt[b][j]);
991 else if (std::find(lt[base].begin(), lt[base].end(), lt[b][j])
992 != lt[base].end())
993 dec_pred(lt[b][j]);
994 else
995 lt[base].push_back(lt[b][j]);
996 }
997 lt.erase(b);
998 }
999
1000 i = le.find(b);
1001 if (i != le.end()) {
1002 for (int j = 0; j < le[b].size(); ++j) {
1003 if (std::find(eq[base].begin(), eq[base].end(), le[b][j]) != eq[base].end())
1004 dec_pred(le[b][j]);
1005 else if (std::find(le[base].begin(), le[base].end(), le[b][j])
1006 != le[base].end())
1007 dec_pred(le[b][j]);
1008 else
1009 le[base].push_back(le[b][j]);
1010 }
1011 le.erase(b);
1012 }
1013
1014 i = pending.find(base);
1015 if (i != pending.end()) {
1016 vector<const indicator_term * > old = pending[base];
1017 pending[base].clear();
1018 for (int j = 0; j < old.size(); ++j) {
1019 if (std::find(eq[base].begin(), eq[base].end(), old[j]) == eq[base].end())
1020 pending[base].push_back(old[j]);
1021 }
1022 }
1023
1024 i = pending.find(b);
1025 if (i != pending.end()) {
1026 for (int j = 0; j < pending[b].size(); ++j) {
1027 if (std::find(eq[base].begin(), eq[base].end(), pending[b][j]) == eq[base].end())
1028 pending[base].push_back(pending[b][j]);
1029 }
1030 pending.erase(b);
1031 }
1032 }
1033
1034 void partial_order::copy(const partial_order& order,
1035 map< const indicator_term *, indicator_term * > old2new)
1036 {
1037 cache.copy(order.cache);
1038
1039 order_type::const_iterator i;
1040 pred_type::const_iterator j;
1041 head_type::const_iterator k;
1042
1043 for (k = order.head.begin(); k != order.head.end(); ++k)
1044 head.insert(old2new[*k]);
1045
1046 for (j = order.pred.begin(); j != order.pred.end(); ++j)
1047 pred[old2new[(*j).first]] = (*j).second;
1048
1049 for (i = order.lt.begin(); i != order.lt.end(); ++i) {
1050 for (int j = 0; j < (*i).second.size(); ++j)
1051 lt[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1052 }
1053 for (i = order.le.begin(); i != order.le.end(); ++i) {
1054 for (int j = 0; j < (*i).second.size(); ++j)
1055 le[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1056 }
1057 for (i = order.eq.begin(); i != order.eq.end(); ++i) {
1058 for (int j = 0; j < (*i).second.size(); ++j)
1059 eq[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1060 }
1061 for (i = order.pending.begin(); i != order.pending.end(); ++i) {
1062 for (int j = 0; j < (*i).second.size(); ++j)
1063 pending[old2new[(*i).first]].push_back(old2new[(*i).second[j]]);
1064 }
1065 }
1066
1067 struct max_term {
1068 EDomain *domain;
1069 vector<evalue *> max;
1070
1071 void print(ostream& os, char **p, barvinok_options *options) const;
1072 void substitute(Matrix *T, barvinok_options *options);
1073 Vector *eval(Value *val, unsigned MaxRays) const;
1074
1075 ~max_term() {
1076 for (int i = 0; i < max.size(); ++i) {
1077 free_evalue_refs(max[i]);
1078 delete max[i];
1079 }
1080 delete domain;
1081 }
1082 };
1083
1084 /*
1085 * Project on first dim dimensions
1086 */
1087 Polyhedron* Polyhedron_Project_Initial(Polyhedron *P, int dim)
1088 {
1089 int i;
1090 Matrix *T;
1091 Polyhedron *I;
1092
1093 if (P->Dimension == dim)
1094 return Polyhedron_Copy(P);
1095
1096 T = Matrix_Alloc(dim+1, P->Dimension+1);
1097 for (i = 0; i < dim; ++i)
1098 value_set_si(T->p[i][i], 1);
1099 value_set_si(T->p[dim][P->Dimension], 1);
1100 I = Polyhedron_Image(P, T, P->NbConstraints);
1101 Matrix_Free(T);
1102 return I;
1103 }
1104
1105 struct indicator {
1106 vector<indicator_term*> term;
1107 indicator_constructor& ic;
1108 partial_order order;
1109 EDomain *D;
1110 Polyhedron *P;
1111 Param_Domain *PD;
1112 lexmin_options *options;
1113 vector<evalue *> substitutions;
1114
1115 indicator(indicator_constructor& ic, Param_Domain *PD, EDomain *D,
1116 lexmin_options *options) :
1117 ic(ic), order(this), D(D), P(NULL), PD(PD), options(options) {}
1118 indicator(const indicator& ind, EDomain *D) :
1119 ic(ind.ic), order(this), D(NULL), P(Polyhedron_Copy(ind.P)),
1120 PD(ind.PD), options(ind.options) {
1121 map< const indicator_term *, indicator_term * > old2new;
1122 for (int i = 0; i < ind.term.size(); ++i) {
1123 indicator_term *it = new indicator_term(*ind.term[i]);
1124 old2new[ind.term[i]] = it;
1125 term.push_back(it);
1126 }
1127 order.copy(ind.order, old2new);
1128 set_domain(D);
1129 }
1130 ~indicator() {
1131 for (int i = 0; i < term.size(); ++i)
1132 delete term[i];
1133 if (D)
1134 delete D;
1135 if (P)
1136 Polyhedron_Free(P);
1137 }
1138
1139 void set_domain(EDomain *D) {
1140 order.cache.clear_transients();
1141 if (this->D)
1142 delete this->D;
1143 this->D = D;
1144 int nparam = ic.PP->Constraints->NbColumns-2 - ic.vertex.length();
1145 if (options->reduce) {
1146 Polyhedron *Q = Polyhedron_Project_Initial(D->D, nparam);
1147 Q = DomainConstraintSimplify(Q, options->verify->barvinok->MaxRays);
1148 if (!P || !PolyhedronIncludes(Q, P))
1149 reduce_in_domain(Q);
1150 if (P)
1151 Polyhedron_Free(P);
1152 P = Q;
1153 order.resort();
1154 }
1155 }
1156
1157 void add(const indicator_term* it);
1158 void remove(const indicator_term* it);
1159 void remove_initial_rational_vertices();
1160 void expand_rational_vertex(const indicator_term *initial);
1161
1162 void print(ostream& os, char **p);
1163 void simplify();
1164 void peel(int i, int j);
1165 void combine(const indicator_term *a, const indicator_term *b);
1166 void add_substitution(evalue *equation);
1167 void perform_pending_substitutions();
1168 void reduce_in_domain(Polyhedron *D);
1169 bool handle_equal_numerators(const indicator_term *base);
1170
1171 max_term* create_max_term(const indicator_term *it);
1172 private:
1173 void substitute(evalue *equation);
1174 };
1175
1176 void partial_order::sanity_check() const
1177 {
1178 order_type::const_iterator i;
1179 order_type::const_iterator prev;
1180 order_type::const_iterator l;
1181 pred_type::const_iterator k, prev_k;
1182
1183 for (k = pred.begin(); k != pred.end(); prev_k = k, ++k)
1184 if (k != pred.begin())
1185 assert(pred.key_comp()((*prev_k).first, (*k).first));
1186 for (i = lt.begin(); i != lt.end(); prev = i, ++i) {
1187 vec_ZZ i_v;
1188 if (ind->D->sample)
1189 i_v = (*i).first->eval(ind->D->sample->p);
1190 if (i != lt.begin())
1191 assert(lt.key_comp()((*prev).first, (*i).first));
1192 l = eq.find((*i).first);
1193 if (l != eq.end())
1194 assert((*l).second.size() > 1);
1195 assert(head.find((*i).first) != head.end() ||
1196 pred.find((*i).first) != pred.end());
1197 for (int j = 0; j < (*i).second.size(); ++j) {
1198 k = pred.find((*i).second[j]);
1199 assert(k != pred.end());
1200 assert((*k).second != 0);
1201 if ((*i).first->sign != 0 &&
1202 (*i).second[j]->sign != 0 && ind->D->sample) {
1203 vec_ZZ j_v = (*i).second[j]->eval(ind->D->sample->p);
1204 assert(lex_cmp(i_v, j_v) < 0);
1205 }
1206 }
1207 }
1208 for (i = le.begin(); i != le.end(); ++i) {
1209 assert((*i).second.size() > 0);
1210 assert(head.find((*i).first) != head.end() ||
1211 pred.find((*i).first) != pred.end());
1212 for (int j = 0; j < (*i).second.size(); ++j) {
1213 k = pred.find((*i).second[j]);
1214 assert(k != pred.end());
1215 assert((*k).second != 0);
1216 }
1217 }
1218 for (i = eq.begin(); i != eq.end(); ++i) {
1219 assert(head.find((*i).first) != head.end() ||
1220 pred.find((*i).first) != pred.end());
1221 assert((*i).second.size() >= 1);
1222 }
1223 for (i = pending.begin(); i != pending.end(); ++i) {
1224 assert(head.find((*i).first) != head.end() ||
1225 pred.find((*i).first) != pred.end());
1226 for (int j = 0; j < (*i).second.size(); ++j)
1227 assert(head.find((*i).second[j]) != head.end() ||
1228 pred.find((*i).second[j]) != pred.end());
1229 }
1230 }
1231
1232 max_term* indicator::create_max_term(const indicator_term *it)
1233 {
1234 int dim = it->den.NumCols();
1235 max_term *maximum = new max_term;
1236 maximum->domain = new EDomain(D);
1237 for (int j = 0; j < dim; ++j) {
1238 evalue *E = new evalue;
1239 value_init(E->d);
1240 evalue_copy(E, it->vertex[j]);
1241 if (evalue_frac2floor_in_domain(E, D->D))
1242 reduce_evalue(E);
1243 maximum->max.push_back(E);
1244 }
1245 return maximum;
1246 }
1247
1248 static order_sign evalue_sign(evalue *diff, EDomain *D, barvinok_options *options)
1249 {
1250 order_sign sign = order_eq;
1251 evalue mone;
1252 value_init(mone.d);
1253 evalue_set_si(&mone, -1, 1);
1254 int len = 1 + D->D->Dimension + 1;
1255 Vector *c = Vector_Alloc(len);
1256 Matrix *T = Matrix_Alloc(2, len-1);
1257
1258 int fract = evalue2constraint(D, diff, c->p, len);
1259 Vector_Copy(c->p+1, T->p[0], len-1);
1260 value_assign(T->p[1][len-2], c->p[0]);
1261
1262 order_sign upper_sign = polyhedron_affine_sign(D->D, T, options);
1263 if (upper_sign == order_lt || !fract)
1264 sign = upper_sign;
1265 else {
1266 emul(&mone, diff);
1267 evalue2constraint(D, diff, c->p, len);
1268 emul(&mone, diff);
1269 Vector_Copy(c->p+1, T->p[0], len-1);
1270 value_assign(T->p[1][len-2], c->p[0]);
1271
1272 order_sign neg_lower_sign = polyhedron_affine_sign(D->D, T, options);
1273
1274 if (neg_lower_sign == order_lt)
1275 sign = order_gt;
1276 else if (neg_lower_sign == order_eq || neg_lower_sign == order_le) {
1277 if (upper_sign == order_eq || upper_sign == order_le)
1278 sign = order_eq;
1279 else
1280 sign = order_ge;
1281 } else {
1282 if (upper_sign == order_lt || upper_sign == order_le ||
1283 upper_sign == order_eq)
1284 sign = order_le;
1285 else
1286 sign = order_unknown;
1287 }
1288 }
1289
1290 Matrix_Free(T);
1291 Vector_Free(c);
1292 free_evalue_refs(&mone);
1293
1294 return sign;
1295 }
1296
1297 /* An auxiliary class that keeps a reference to an evalue
1298 * and frees it when it goes out of scope.
1299 */
1300 struct temp_evalue {
1301 evalue *E;
1302 temp_evalue() : E(NULL) {}
1303 temp_evalue(evalue *e) : E(e) {}
1304 operator evalue* () const { return E; }
1305 evalue *operator=(evalue *e) {
1306 if (E) {
1307 free_evalue_refs(E);
1308 delete E;
1309 }
1310 E = e;
1311 return E;
1312 }
1313 ~temp_evalue() {
1314 if (E) {
1315 free_evalue_refs(E);
1316 delete E;
1317 }
1318 }
1319 };
1320
1321 static void substitute(vector<indicator_term*>& term, evalue *equation)
1322 {
1323 evalue *fract = NULL;
1324 evalue *val = new evalue;
1325 value_init(val->d);
1326 evalue_copy(val, equation);
1327
1328 evalue factor;
1329 value_init(factor.d);
1330 value_init(factor.x.n);
1331
1332 evalue *e;
1333 for (e = val; value_zero_p(e->d) && e->x.p->type != fractional;
1334 e = &e->x.p->arr[type_offset(e->x.p)])
1335 ;
1336
1337 if (value_zero_p(e->d) && e->x.p->type == fractional)
1338 fract = &e->x.p->arr[0];
1339 else {
1340 for (e = val; value_zero_p(e->d) && e->x.p->type != polynomial;
1341 e = &e->x.p->arr[type_offset(e->x.p)])
1342 ;
1343 assert(value_zero_p(e->d) && e->x.p->type == polynomial);
1344 }
1345
1346 int offset = type_offset(e->x.p);
1347
1348 assert(value_notzero_p(e->x.p->arr[offset+1].d));
1349 assert(value_notzero_p(e->x.p->arr[offset+1].x.n));
1350 if (value_neg_p(e->x.p->arr[offset+1].x.n)) {
1351 value_oppose(factor.d, e->x.p->arr[offset+1].x.n);
1352 value_assign(factor.x.n, e->x.p->arr[offset+1].d);
1353 } else {
1354 value_assign(factor.d, e->x.p->arr[offset+1].x.n);
1355 value_oppose(factor.x.n, e->x.p->arr[offset+1].d);
1356 }
1357
1358 free_evalue_refs(&e->x.p->arr[offset+1]);
1359 enode *p = e->x.p;
1360 value_clear(e->d);
1361 *e = e->x.p->arr[offset];
1362
1363 emul(&factor, val);
1364
1365 if (fract) {
1366 for (int i = 0; i < term.size(); ++i)
1367 term[i]->substitute(fract, val);
1368
1369 free_evalue_refs(&p->arr[0]);
1370 } else {
1371 for (int i = 0; i < term.size(); ++i)
1372 term[i]->substitute(p->pos, val);
1373 }
1374
1375 free_evalue_refs(&factor);
1376 free_evalue_refs(val);
1377 delete val;
1378
1379 free(p);
1380 }
1381
1382 order_sign partial_order::compare(const indicator_term *a, const indicator_term *b)
1383 {
1384 unsigned dim = a->den.NumCols();
1385 order_sign sign = order_eq;
1386 bool rational = a->sign == 0 || b->sign == 0;
1387
1388 order_sign cached_sign = order_eq;
1389 for (int k = 0; k < dim; ++k) {
1390 cached_sign = evalue_diff_constant_sign(a->vertex[k], b->vertex[k]);
1391 if (cached_sign != order_eq)
1392 break;
1393 }
1394 if (cached_sign != order_undefined)
1395 return cached_sign;
1396
1397 order_cache_el cache_el;
1398 cached_sign = order_undefined;
1399 if (!rational)
1400 cached_sign = cache.check(a, b, cache_el);
1401 if (cached_sign != order_undefined) {
1402 cache_el.free();
1403 return cached_sign;
1404 }
1405
1406 sign = order_eq;
1407
1408 vector<indicator_term *> term;
1409
1410 for (int k = 0; k < dim; ++k) {
1411 /* compute a->vertex[k] - b->vertex[k] */
1412 evalue *diff;
1413 if (cache_el.e.size() <= k) {
1414 diff = ediff(a->vertex[k], b->vertex[k]);
1415 cache_el.e.push_back(diff);
1416 } else
1417 diff = cache_el.e[k];
1418 temp_evalue tdiff;
1419 if (term.size())
1420 tdiff = diff = ediff(term[0]->vertex[k], term[1]->vertex[k]);
1421 order_sign diff_sign;
1422 if (eequal(a->vertex[k], b->vertex[k]))
1423 diff_sign = order_eq;
1424 else
1425 diff_sign = evalue_sign(diff, ind->D,
1426 ind->options->verify->barvinok);
1427
1428 if (diff_sign == order_undefined) {
1429 assert(sign == order_le || sign == order_ge);
1430 if (sign == order_le)
1431 sign = order_lt;
1432 else
1433 sign = order_gt;
1434 break;
1435 }
1436 if (diff_sign == order_lt) {
1437 if (sign == order_eq || sign == order_le)
1438 sign = order_lt;
1439 else
1440 sign = order_unknown;
1441 break;
1442 }
1443 if (diff_sign == order_gt) {
1444 if (sign == order_eq || sign == order_ge)
1445 sign = order_gt;
1446 else
1447 sign = order_unknown;
1448 break;
1449 }
1450 if (diff_sign == order_eq) {
1451 if (term.size() == 0 && !rational && !EVALUE_IS_ZERO(*diff))
1452 ind->add_substitution(diff);
1453 continue;
1454 }
1455 if ((diff_sign == order_unknown) ||
1456 ((diff_sign == order_lt || diff_sign == order_le) && sign == order_ge) ||
1457 ((diff_sign == order_gt || diff_sign == order_ge) && sign == order_le)) {
1458 sign = order_unknown;
1459 break;
1460 }
1461
1462 sign = diff_sign;
1463
1464 if (!term.size()) {
1465 term.push_back(new indicator_term(*a));
1466 term.push_back(new indicator_term(*b));
1467 }
1468 substitute(term, diff);
1469 }
1470
1471 if (!rational)
1472 cache.add(cache_el, sign);
1473 else
1474 cache_el.free();
1475
1476 if (term.size()) {
1477 delete term[0];
1478 delete term[1];
1479 }
1480
1481 return sign;
1482 }
1483
1484 bool partial_order::compared(const indicator_term* a, const indicator_term* b)
1485 {
1486 order_type::iterator j;
1487
1488 j = lt.find(a);
1489 if (j != lt.end() && std::find(lt[a].begin(), lt[a].end(), b) != lt[a].end())
1490 return true;
1491
1492 j = le.find(a);
1493 if (j != le.end() && std::find(le[a].begin(), le[a].end(), b) != le[a].end())
1494 return true;
1495
1496 return false;
1497 }
1498
1499 void partial_order::add(const indicator_term* it,
1500 std::set<const indicator_term *> *filter)
1501 {
1502 if (eq.find(it) != eq.end() && eq[it].size() == 1)
1503 return;
1504
1505 head_type head_copy(head);
1506
1507 if (!filter)
1508 head.insert(it);
1509
1510 head_type::iterator i;
1511 for (i = head_copy.begin(); i != head_copy.end(); ++i) {
1512 if (*i == it)
1513 continue;
1514 if (eq.find(*i) != eq.end() && eq[*i].size() == 1)
1515 continue;
1516 if (filter) {
1517 if (filter->find(*i) == filter->end())
1518 continue;
1519 if (compared(*i, it))
1520 continue;
1521 }
1522 order_sign sign = compare(it, *i);
1523 if (sign == order_lt) {
1524 lt[it].push_back(*i);
1525 inc_pred(*i);
1526 } else if (sign == order_le) {
1527 le[it].push_back(*i);
1528 inc_pred(*i);
1529 } else if (sign == order_eq) {
1530 set_equal(it, *i);
1531 return;
1532 } else if (sign == order_gt) {
1533 pending[*i].push_back(it);
1534 lt[*i].push_back(it);
1535 inc_pred(it);
1536 } else if (sign == order_ge) {
1537 pending[*i].push_back(it);
1538 le[*i].push_back(it);
1539 inc_pred(it);
1540 }
1541 }
1542 }
1543
1544 void partial_order::remove(const indicator_term* it)
1545 {
1546 std::set<const indicator_term *> filter;
1547 order_type::iterator i;
1548
1549 assert(head.find(it) != head.end());
1550
1551 i = eq.find(it);
1552 if (i != eq.end()) {
1553 assert(eq[it].size() >= 1);
1554 const indicator_term *base;
1555 if (eq[it].size() == 1) {
1556 base = eq[it][0];
1557 eq.erase(it);
1558
1559 vector<const indicator_term * >::iterator j;
1560 j = std::find(eq[base].begin(), eq[base].end(), it);
1561 assert(j != eq[base].end());
1562 eq[base].erase(j);
1563 } else {
1564 /* "it" may no longer be the smallest, since the order
1565 * structure may have been copied from another one.
1566 */
1567 std::sort(eq[it].begin()+1, eq[it].end(), pred.key_comp());
1568 assert(eq[it][0] == it);
1569 eq[it].erase(eq[it].begin());
1570 base = eq[it][0];
1571 eq[base] = eq[it];
1572 eq.erase(it);
1573
1574 for (int j = 1; j < eq[base].size(); ++j)
1575 eq[eq[base][j]][0] = base;
1576
1577 i = lt.find(it);
1578 if (i != lt.end()) {
1579 lt[base] = lt[it];
1580 lt.erase(it);
1581 }
1582
1583 i = le.find(it);
1584 if (i != le.end()) {
1585 le[base] = le[it];
1586 le.erase(it);
1587 }
1588
1589 i = pending.find(it);
1590 if (i != pending.end()) {
1591 pending[base] = pending[it];
1592 pending.erase(it);
1593 }
1594 }
1595
1596 if (eq[base].size() == 1)
1597 eq.erase(base);
1598
1599 head.erase(it);
1600
1601 return;
1602 }
1603
1604 i = lt.find(it);
1605 if (i != lt.end()) {
1606 for (int j = 0; j < lt[it].size(); ++j) {
1607 filter.insert(lt[it][j]);
1608 dec_pred(lt[it][j]);
1609 }
1610 lt.erase(it);
1611 }
1612
1613 i = le.find(it);
1614 if (i != le.end()) {
1615 for (int j = 0; j < le[it].size(); ++j) {
1616 filter.insert(le[it][j]);
1617 dec_pred(le[it][j]);
1618 }
1619 le.erase(it);
1620 }
1621
1622 head.erase(it);
1623
1624 i = pending.find(it);
1625 if (i != pending.end()) {
1626 vector<const indicator_term *> it_pending = pending[it];
1627 pending.erase(it);
1628 for (int j = 0; j < it_pending.size(); ++j) {
1629 filter.erase(it_pending[j]);
1630 add(it_pending[j], &filter);
1631 }
1632 }
1633 }
1634
1635 void partial_order::print(ostream& os, char **p)
1636 {
1637 order_type::iterator i;
1638 pred_type::iterator j;
1639 head_type::iterator k;
1640 for (k = head.begin(); k != head.end(); ++k) {
1641 (*k)->print(os, p);
1642 os << endl;
1643 }
1644 for (j = pred.begin(); j != pred.end(); ++j) {
1645 (*j).first->print(os, p);
1646 os << ": " << (*j).second << endl;
1647 }
1648 for (i = lt.begin(); i != lt.end(); ++i) {
1649 (*i).first->print(os, p);
1650 assert(head.find((*i).first) != head.end() ||
1651 pred.find((*i).first) != pred.end());
1652 if (pred.find((*i).first) != pred.end())
1653 os << "(" << pred[(*i).first] << ")";
1654 os << " < ";
1655 for (int j = 0; j < (*i).second.size(); ++j) {
1656 if (j)
1657 os << ", ";
1658 (*i).second[j]->print(os, p);
1659 assert(pred.find((*i).second[j]) != pred.end());
1660 os << "(" << pred[(*i).second[j]] << ")";
1661 }
1662 os << endl;
1663 }
1664 for (i = le.begin(); i != le.end(); ++i) {
1665 (*i).first->print(os, p);
1666 assert(head.find((*i).first) != head.end() ||
1667 pred.find((*i).first) != pred.end());
1668 if (pred.find((*i).first) != pred.end())
1669 os << "(" << pred[(*i).first] << ")";
1670 os << " <= ";
1671 for (int j = 0; j < (*i).second.size(); ++j) {
1672 if (j)
1673 os << ", ";
1674 (*i).second[j]->print(os, p);
1675 assert(pred.find((*i).second[j]) != pred.end());
1676 os << "(" << pred[(*i).second[j]] << ")";
1677 }
1678 os << endl;
1679 }
1680 for (i = eq.begin(); i != eq.end(); ++i) {
1681 if ((*i).second.size() <= 1)
1682 continue;
1683 (*i).first->print(os, p);
1684 assert(head.find((*i).first) != head.end() ||
1685 pred.find((*i).first) != pred.end());
1686 if (pred.find((*i).first) != pred.end())
1687 os << "(" << pred[(*i).first] << ")";
1688 for (int j = 1; j < (*i).second.size(); ++j) {
1689 if (j)
1690 os << " = ";
1691 (*i).second[j]->print(os, p);
1692 assert(head.find((*i).second[j]) != head.end() ||
1693 pred.find((*i).second[j]) != pred.end());
1694 if (pred.find((*i).second[j]) != pred.end())
1695 os << "(" << pred[(*i).second[j]] << ")";
1696 }
1697 os << endl;
1698 }
1699 for (i = pending.begin(); i != pending.end(); ++i) {
1700 os << "pending on ";
1701 (*i).first->print(os, p);
1702 assert(head.find((*i).first) != head.end() ||
1703 pred.find((*i).first) != pred.end());
1704 if (pred.find((*i).first) != pred.end())
1705 os << "(" << pred[(*i).first] << ")";
1706 os << ": ";
1707 for (int j = 0; j < (*i).second.size(); ++j) {
1708 if (j)
1709 os << ", ";
1710 (*i).second[j]->print(os, p);
1711 assert(pred.find((*i).second[j]) != pred.end());
1712 os << "(" << pred[(*i).second[j]] << ")";
1713 }
1714 os << endl;
1715 }
1716 }
1717
1718 void indicator::add(const indicator_term* it)
1719 {
1720 indicator_term *nt = new indicator_term(*it);
1721 if (options->reduce)
1722 nt->reduce_in_domain(P ? P : D->D);
1723 term.push_back(nt);
1724 order.add(nt, NULL);
1725 assert(term.size() == order.head.size() + order.pred.size());
1726 }
1727
1728 void indicator::remove(const indicator_term* it)
1729 {
1730 vector<indicator_term *>::iterator i;
1731 i = std::find(term.begin(), term.end(), it);
1732 assert(i!= term.end());
1733 order.remove(it);
1734 term.erase(i);
1735 assert(term.size() == order.head.size() + order.pred.size());
1736 delete it;
1737 }
1738
1739 void indicator::expand_rational_vertex(const indicator_term *initial)
1740 {
1741 int pos = initial->pos;
1742 remove(initial);
1743 if (ic.terms[pos].size() == 0) {
1744 Param_Vertices *V;
1745 FORALL_PVertex_in_ParamPolyhedron(V, PD, ic.PP) // _i is internal counter
1746 if (_i == pos) {
1747 ic.decompose_at_vertex(V, pos, options->verify->barvinok);
1748 break;
1749 }
1750 END_FORALL_PVertex_in_ParamPolyhedron;
1751 }
1752 for (int j = 0; j < ic.terms[pos].size(); ++j)
1753 add(ic.terms[pos][j]);
1754 }
1755
1756 void indicator::remove_initial_rational_vertices()
1757 {
1758 do {
1759 const indicator_term *initial = NULL;
1760 partial_order::head_type::iterator i;
1761 for (i = order.head.begin(); i != order.head.end(); ++i) {
1762 if ((*i)->sign != 0)
1763 continue;
1764 if (order.eq.find(*i) != order.eq.end() && order.eq[*i].size() <= 1)
1765 continue;
1766 initial = *i;
1767 break;
1768 }
1769 if (!initial)
1770 break;
1771 expand_rational_vertex(initial);
1772 } while(1);
1773 }
1774
1775 void indicator::reduce_in_domain(Polyhedron *D)
1776 {
1777 for (int i = 0; i < term.size(); ++i)
1778 term[i]->reduce_in_domain(D);
1779 }
1780
1781 void indicator::print(ostream& os, char **p)
1782 {
1783 assert(term.size() == order.head.size() + order.pred.size());
1784 for (int i = 0; i < term.size(); ++i) {
1785 term[i]->print(os, p);
1786 if (D->sample) {
1787 os << ": " << term[i]->eval(D->sample->p);
1788 }
1789 os << endl;
1790 }
1791 order.print(os, p);
1792 }
1793
1794 /* Remove pairs of opposite terms */
1795 void indicator::simplify()
1796 {
1797 for (int i = 0; i < term.size(); ++i) {
1798 for (int j = i+1; j < term.size(); ++j) {
1799 if (term[i]->sign + term[j]->sign != 0)
1800 continue;
1801 if (term[i]->den != term[j]->den)
1802 continue;
1803 int k;
1804 for (k = 0; k < term[i]->den.NumCols(); ++k)
1805 if (!eequal(term[i]->vertex[k], term[j]->vertex[k]))
1806 break;
1807 if (k < term[i]->den.NumCols())
1808 continue;
1809 delete term[i];
1810 delete term[j];
1811 term.erase(term.begin()+j);
1812 term.erase(term.begin()+i);
1813 --i;
1814 break;
1815 }
1816 }
1817 }
1818
1819 void indicator::peel(int i, int j)
1820 {
1821 if (j < i) {
1822 int tmp = i;
1823 i = j;
1824 j = tmp;
1825 }
1826 assert (i < j);
1827 int dim = term[i]->den.NumCols();
1828
1829 mat_ZZ common;
1830 mat_ZZ rest_i;
1831 mat_ZZ rest_j;
1832 int n_common = 0, n_i = 0, n_j = 0;
1833
1834 common.SetDims(min(term[i]->den.NumRows(), term[j]->den.NumRows()), dim);
1835 rest_i.SetDims(term[i]->den.NumRows(), dim);
1836 rest_j.SetDims(term[j]->den.NumRows(), dim);
1837
1838 int k, l;
1839 for (k = 0, l = 0; k < term[i]->den.NumRows() && l < term[j]->den.NumRows(); ) {
1840 int s = lex_cmp(term[i]->den[k], term[j]->den[l]);
1841 if (s == 0) {
1842 common[n_common++] = term[i]->den[k];
1843 ++k;
1844 ++l;
1845 } else if (s < 0)
1846 rest_i[n_i++] = term[i]->den[k++];
1847 else
1848 rest_j[n_j++] = term[j]->den[l++];
1849 }
1850 while (k < term[i]->den.NumRows())
1851 rest_i[n_i++] = term[i]->den[k++];
1852 while (l < term[j]->den.NumRows())
1853 rest_j[n_j++] = term[j]->den[l++];
1854 common.SetDims(n_common, dim);
1855 rest_i.SetDims(n_i, dim);
1856 rest_j.SetDims(n_j, dim);
1857
1858 for (k = 0; k <= n_i; ++k) {
1859 indicator_term *it = new indicator_term(*term[i]);
1860 it->den.SetDims(n_common + k, dim);
1861 for (l = 0; l < n_common; ++l)
1862 it->den[l] = common[l];
1863 for (l = 0; l < k; ++l)
1864 it->den[n_common+l] = rest_i[l];
1865 lex_order_rows(it->den);
1866 if (k)
1867 for (l = 0; l < dim; ++l)
1868 evalue_add_constant(it->vertex[l], rest_i[k-1][l]);
1869 term.push_back(it);
1870 }
1871
1872 for (k = 0; k <= n_j; ++k) {
1873 indicator_term *it = new indicator_term(*term[j]);
1874 it->den.SetDims(n_common + k, dim);
1875 for (l = 0; l < n_common; ++l)
1876 it->den[l] = common[l];
1877 for (l = 0; l < k; ++l)
1878 it->den[n_common+l] = rest_j[l];
1879 lex_order_rows(it->den);
1880 if (k)
1881 for (l = 0; l < dim; ++l)
1882 evalue_add_constant(it->vertex[l], rest_j[k-1][l]);
1883 term.push_back(it);
1884 }
1885 term.erase(term.begin()+j);
1886 term.erase(term.begin()+i);
1887 }
1888
1889 void indicator::combine(const indicator_term *a, const indicator_term *b)
1890 {
1891 int dim = a->den.NumCols();
1892
1893 mat_ZZ common;
1894 mat_ZZ rest_i; /* factors in a, but not in b */
1895 mat_ZZ rest_j; /* factors in b, but not in a */
1896 int n_common = 0, n_i = 0, n_j = 0;
1897
1898 common.SetDims(min(a->den.NumRows(), b->den.NumRows()), dim);
1899 rest_i.SetDims(a->den.NumRows(), dim);
1900 rest_j.SetDims(b->den.NumRows(), dim);
1901
1902 int k, l;
1903 for (k = 0, l = 0; k < a->den.NumRows() && l < b->den.NumRows(); ) {
1904 int s = lex_cmp(a->den[k], b->den[l]);
1905 if (s == 0) {
1906 common[n_common++] = a->den[k];
1907 ++k;
1908 ++l;
1909 } else if (s < 0)
1910 rest_i[n_i++] = a->den[k++];
1911 else
1912 rest_j[n_j++] = b->den[l++];
1913 }
1914 while (k < a->den.NumRows())
1915 rest_i[n_i++] = a->den[k++];
1916 while (l < b->den.NumRows())
1917 rest_j[n_j++] = b->den[l++];
1918 common.SetDims(n_common, dim);
1919 rest_i.SetDims(n_i, dim);
1920 rest_j.SetDims(n_j, dim);
1921
1922 assert(order.eq[a].size() > 1);
1923 indicator_term *prev;
1924
1925 prev = NULL;
1926 for (int k = n_i-1; k >= 0; --k) {
1927 indicator_term *it = new indicator_term(*b);
1928 it->den.SetDims(n_common + n_j + n_i-k, dim);
1929 for (int l = k; l < n_i; ++l)
1930 it->den[n_common+n_j+l-k] = rest_i[l];
1931 lex_order_rows(it->den);
1932 for (int m = 0; m < dim; ++m)
1933 evalue_add_constant(it->vertex[m], rest_i[k][m]);
1934 it->sign = -it->sign;
1935 if (prev) {
1936 order.pending[it].push_back(prev);
1937 order.lt[it].push_back(prev);
1938 order.inc_pred(prev);
1939 }
1940 term.push_back(it);
1941 order.head.insert(it);
1942 prev = it;
1943 }
1944 if (n_i) {
1945 indicator_term *it = new indicator_term(*b);
1946 it->den.SetDims(n_common + n_i + n_j, dim);
1947 for (l = 0; l < n_i; ++l)
1948 it->den[n_common+n_j+l] = rest_i[l];
1949 lex_order_rows(it->den);
1950 assert(prev);
1951 order.pending[a].push_back(prev);
1952 order.lt[a].push_back(prev);
1953 order.inc_pred(prev);
1954 order.replace(b, it);
1955 delete it;
1956 }
1957
1958 prev = NULL;
1959 for (int k = n_j-1; k >= 0; --k) {
1960 indicator_term *it = new indicator_term(*a);
1961 it->den.SetDims(n_common + n_i + n_j-k, dim);
1962 for (int l = k; l < n_j; ++l)
1963 it->den[n_common+n_i+l-k] = rest_j[l];
1964 lex_order_rows(it->den);
1965 for (int m = 0; m < dim; ++m)
1966 evalue_add_constant(it->vertex[m], rest_j[k][m]);
1967 it->sign = -it->sign;
1968 if (prev) {
1969 order.pending[it].push_back(prev);
1970 order.lt[it].push_back(prev);
1971 order.inc_pred(prev);
1972 }
1973 term.push_back(it);
1974 order.head.insert(it);
1975 prev = it;
1976 }
1977 if (n_j) {
1978 indicator_term *it = new indicator_term(*a);
1979 it->den.SetDims(n_common + n_i + n_j, dim);
1980 for (l = 0; l < n_j; ++l)
1981 it->den[n_common+n_i+l] = rest_j[l];
1982 lex_order_rows(it->den);
1983 assert(prev);
1984 order.pending[a].push_back(prev);
1985 order.lt[a].push_back(prev);
1986 order.inc_pred(prev);
1987 order.replace(a, it);
1988 delete it;
1989 }
1990
1991 assert(term.size() == order.head.size() + order.pred.size());
1992 }
1993
1994 bool indicator::handle_equal_numerators(const indicator_term *base)
1995 {
1996 for (int i = 0; i < order.eq[base].size(); ++i) {
1997 for (int j = i+1; j < order.eq[base].size(); ++j) {
1998 if (order.eq[base][i]->is_opposite(order.eq[base][j])) {
1999 remove(order.eq[base][j]);
2000 remove(i ? order.eq[base][i] : base);
2001 return true;
2002 }
2003 }
2004 }
2005 for (int j = 1; j < order.eq[base].size(); ++j)
2006 if (order.eq[base][j]->sign != base->sign) {
2007 combine(base, order.eq[base][j]);
2008 return true;
2009 }
2010 return false;
2011 }
2012
2013 void indicator::substitute(evalue *equation)
2014 {
2015 ::substitute(term, equation);
2016 }
2017
2018 void indicator::add_substitution(evalue *equation)
2019 {
2020 for (int i = 0; i < substitutions.size(); ++i)
2021 if (eequal(substitutions[i], equation))
2022 return;
2023 evalue *copy = new evalue();
2024 value_init(copy->d);
2025 evalue_copy(copy, equation);
2026 substitutions.push_back(copy);
2027 }
2028
2029 void indicator::perform_pending_substitutions()
2030 {
2031 if (substitutions.size() == 0)
2032 return;
2033
2034 for (int i = 0; i < substitutions.size(); ++i) {
2035 substitute(substitutions[i]);
2036 free_evalue_refs(substitutions[i]);
2037 delete substitutions[i];
2038 }
2039 substitutions.clear();
2040 order.resort();
2041 }
2042
2043 static void print_varlist(ostream& os, int n, char **names)
2044 {
2045 int i;
2046 os << "[";
2047 for (i = 0; i < n; ++i) {
2048 if (i)
2049 os << ",";
2050 os << names[i];
2051 }
2052 os << "]";
2053 }
2054
2055 void max_term::print(ostream& os, char **p, barvinok_options *options) const
2056 {
2057 os << "{ ";
2058 print_varlist(os, domain->dimension(), p);
2059 os << " -> ";
2060 os << "[";
2061 for (int i = 0; i < max.size(); ++i) {
2062 if (i)
2063 os << ",";
2064 evalue_print(os, max[i], p);
2065 }
2066 os << "]";
2067 os << " : ";
2068 domain->print_constraints(os, p, options);
2069 os << " }" << endl;
2070 }
2071
2072 /* T maps the compressed parameters to the original parameters,
2073 * while this max_term is based on the compressed parameters
2074 * and we want get the original parameters back.
2075 */
2076 void max_term::substitute(Matrix *T, barvinok_options *options)
2077 {
2078 assert(domain->dimension() == T->NbColumns-1);
2079 Matrix *Eq;
2080 Matrix *inv = left_inverse(T, &Eq);
2081
2082 evalue denom;
2083 value_init(denom.d);
2084 value_init(denom.x.n);
2085 value_set_si(denom.x.n, 1);
2086 value_assign(denom.d, inv->p[inv->NbRows-1][inv->NbColumns-1]);
2087
2088 vec_ZZ v;
2089 v.SetLength(inv->NbColumns);
2090 evalue **subs = new evalue *[inv->NbRows-1];
2091 for (int i = 0; i < inv->NbRows-1; ++i) {
2092 values2zz(inv->p[i], v, v.length());
2093 subs[i] = multi_monom(v);
2094 emul(&denom, subs[i]);
2095 }
2096 free_evalue_refs(&denom);
2097
2098 domain->substitute(subs, inv, Eq, options->MaxRays);
2099 Matrix_Free(Eq);
2100
2101 for (int i = 0; i < max.size(); ++i) {
2102 evalue_substitute(max[i], subs);
2103 reduce_evalue(max[i]);
2104 }
2105
2106 for (int i = 0; i < inv->NbRows-1; ++i) {
2107 free_evalue_refs(subs[i]);
2108 delete subs[i];
2109 }
2110 delete [] subs;
2111 Matrix_Free(inv);
2112 }
2113
2114 Vector *max_term::eval(Value *val, unsigned MaxRays) const
2115 {
2116 if (!domain->contains(val, domain->dimension()))
2117 return NULL;
2118 Vector *res = Vector_Alloc(max.size());
2119 for (int i = 0; i < max.size(); ++i) {
2120 compute_evalue(max[i], val, &res->p[i]);
2121 }
2122 return res;
2123 }
2124
2125 struct split {
2126 evalue *constraint;
2127 enum sign { le, ge, lge } sign;
2128
2129 split (evalue *c, enum sign s) : constraint(c), sign(s) {}
2130 };
2131
2132 static void split_on(const split& sp, EDomain *D,
2133 EDomain **Dlt, EDomain **Deq, EDomain **Dgt,
2134 lexmin_options *options)
2135 {
2136 EDomain *ED[3];
2137 *Dlt = NULL;
2138 *Deq = NULL;
2139 *Dgt = NULL;
2140 ge_constraint *ge = D->compute_ge_constraint(sp.constraint);
2141 if (sp.sign == split::lge || sp.sign == split::ge)
2142 ED[2] = EDomain::new_from_ge_constraint(ge, 1, options->verify->barvinok);
2143 else
2144 ED[2] = NULL;
2145 if (sp.sign == split::lge || sp.sign == split::le)
2146 ED[0] = EDomain::new_from_ge_constraint(ge, -1, options->verify->barvinok);
2147 else
2148 ED[0] = NULL;
2149
2150 assert(sp.sign == split::lge || sp.sign == split::ge || sp.sign == split::le);
2151 ED[1] = EDomain::new_from_ge_constraint(ge, 0, options->verify->barvinok);
2152
2153 delete ge;
2154
2155 for (int i = 0; i < 3; ++i) {
2156 if (!ED[i])
2157 continue;
2158 if (D->sample && ED[i]->contains(D->sample->p, D->sample->Size-1)) {
2159 ED[i]->sample = Vector_Alloc(D->sample->Size);
2160 Vector_Copy(D->sample->p, ED[i]->sample->p, D->sample->Size);
2161 } else if (emptyQ2(ED[i]->D) ||
2162 (options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2163 !(ED[i]->not_empty(options)))) {
2164 delete ED[i];
2165 ED[i] = NULL;
2166 }
2167 }
2168 *Dlt = ED[0];
2169 *Deq = ED[1];
2170 *Dgt = ED[2];
2171 }
2172
2173 ostream & operator<< (ostream & os, const vector<int> & v)
2174 {
2175 os << "[";
2176 for (int i = 0; i < v.size(); ++i) {
2177 if (i)
2178 os << ", ";
2179 os << v[i];
2180 }
2181 os << "]";
2182 return os;
2183 }
2184
2185 void construct_rational_vertices(Param_Polyhedron *PP, Matrix *T, unsigned dim,
2186 int nparam, vector<indicator_term *>& vertices)
2187 {
2188 int i;
2189 Param_Vertices *PV;
2190 Value lcm, tmp;
2191 value_init(lcm);
2192 value_init(tmp);
2193
2194 vec_ZZ v;
2195 v.SetLength(nparam+1);
2196
2197 evalue factor;
2198 value_init(factor.d);
2199 value_init(factor.x.n);
2200 value_set_si(factor.x.n, 1);
2201 value_set_si(factor.d, 1);
2202
2203 for (i = 0, PV = PP->V; PV; ++i, PV = PV->next) {
2204 indicator_term *term = new indicator_term(dim, i);
2205 vertices.push_back(term);
2206 Matrix *M = Matrix_Alloc(PV->Vertex->NbRows+nparam+1, nparam+1);
2207 value_set_si(lcm, 1);
2208 for (int j = 0; j < PV->Vertex->NbRows; ++j)
2209 value_lcm(lcm, lcm, PV->Vertex->p[j][nparam+1]);
2210 value_assign(M->p[M->NbRows-1][M->NbColumns-1], lcm);
2211 for (int j = 0; j < PV->Vertex->NbRows; ++j) {
2212 value_division(tmp, lcm, PV->Vertex->p[j][nparam+1]);
2213 Vector_Scale(PV->Vertex->p[j], M->p[j], tmp, nparam+1);
2214 }
2215 for (int j = 0; j < nparam; ++j)
2216 value_assign(M->p[PV->Vertex->NbRows+j][j], lcm);
2217 if (T) {
2218 Matrix *M2 = Matrix_Alloc(T->NbRows, M->NbColumns);
2219 Matrix_Product(T, M, M2);
2220 Matrix_Free(M);
2221 M = M2;
2222 }
2223 for (int j = 0; j < dim; ++j) {
2224 values2zz(M->p[j], v, nparam+1);
2225 term->vertex[j] = multi_monom(v);
2226 value_assign(factor.d, lcm);
2227 emul(&factor, term->vertex[j]);
2228 }
2229 Matrix_Free(M);
2230 }
2231 assert(i == PP->nbV);
2232 free_evalue_refs(&factor);
2233 value_clear(lcm);
2234 value_clear(tmp);
2235 }
2236
2237 static vector<max_term*> lexmin(indicator& ind, unsigned nparam,
2238 vector<int> loc)
2239 {
2240 vector<max_term*> maxima;
2241 partial_order::head_type::iterator i;
2242 vector<int> best_score;
2243 vector<int> second_score;
2244 vector<int> neg_score;
2245
2246 do {
2247 ind.perform_pending_substitutions();
2248 const indicator_term *best = NULL, *second = NULL, *neg = NULL,
2249 *neg_eq = NULL, *neg_le = NULL;
2250 for (i = ind.order.head.begin(); i != ind.order.head.end(); ++i) {
2251 vector<int> score;
2252 const indicator_term *term = *i;
2253 if (term->sign == 0) {
2254 ind.expand_rational_vertex(term);
2255 break;
2256 }
2257
2258 if (ind.order.eq.find(term) != ind.order.eq.end()) {
2259 int j;
2260 if (ind.order.eq[term].size() <= 1)
2261 continue;
2262 for (j = 1; j < ind.order.eq[term].size(); ++j)
2263 if (ind.order.pred.find(ind.order.eq[term][j]) !=
2264 ind.order.pred.end())
2265 break;
2266 if (j < ind.order.eq[term].size())
2267 continue;
2268 score.push_back(ind.order.eq[term].size());
2269 } else
2270 score.push_back(0);
2271 if (ind.order.le.find(term) != ind.order.le.end())
2272 score.push_back(ind.order.le[term].size());
2273 else
2274 score.push_back(0);
2275 if (ind.order.lt.find(term) != ind.order.lt.end())
2276 score.push_back(-ind.order.lt[term].size());
2277 else
2278 score.push_back(0);
2279
2280 if (term->sign > 0) {
2281 if (!best || score < best_score) {
2282 second = best;
2283 second_score = best_score;
2284 best = term;
2285 best_score = score;
2286 } else if (!second || score < second_score) {
2287 second = term;
2288 second_score = score;
2289 }
2290 } else {
2291 if (!neg_eq && ind.order.eq.find(term) != ind.order.eq.end()) {
2292 for (int j = 1; j < ind.order.eq[term].size(); ++j)
2293 if (ind.order.eq[term][j]->sign != term->sign) {
2294 neg_eq = term;
2295 break;
2296 }
2297 }
2298 if (!neg_le && ind.order.le.find(term) != ind.order.le.end())
2299 neg_le = term;
2300 if (!neg || score < neg_score) {
2301 neg = term;
2302 neg_score = score;
2303 }
2304 }
2305 }
2306 if (i != ind.order.head.end())
2307 continue;
2308
2309 if (!best && neg_eq) {
2310 assert(ind.order.eq[neg_eq].size() != 0);
2311 bool handled = ind.handle_equal_numerators(neg_eq);
2312 assert(handled);
2313 continue;
2314 }
2315
2316 if (!best && neg_le) {
2317 /* The smallest term is negative and <= some positive term */
2318 best = neg_le;
2319 neg = NULL;
2320 }
2321
2322 if (!best) {
2323 /* apparently there can be negative initial term on empty domains */
2324 if (ind.options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2325 ind.options->verify->barvinok->lp_solver == BV_LP_POLYLIB)
2326 assert(!neg);
2327 break;
2328 }
2329
2330 if (!second && !neg) {
2331 const indicator_term *rat = NULL;
2332 assert(best);
2333 if (ind.order.le.find(best) == ind.order.le.end()) {
2334 if (ind.order.eq.find(best) != ind.order.eq.end()) {
2335 bool handled = ind.handle_equal_numerators(best);
2336 if (ind.options->emptiness_check !=
2337 BV_LEXMIN_EMPTINESS_CHECK_NONE &&
2338 ind.options->verify->barvinok->lp_solver == BV_LP_POLYLIB)
2339 assert(handled);
2340 /* If !handled then the leading coefficient is bigger than one;
2341 * must be an empty domain
2342 */
2343 if (handled)
2344 continue;
2345 else
2346 break;
2347 }
2348 maxima.push_back(ind.create_max_term(best));
2349 break;
2350 }
2351 for (int j = 0; j < ind.order.le[best].size(); ++j) {
2352 if (ind.order.le[best][j]->sign == 0) {
2353 if (!rat && ind.order.pred[ind.order.le[best][j]] == 1)
2354 rat = ind.order.le[best][j];
2355 } else if (ind.order.le[best][j]->sign > 0) {
2356 second = ind.order.le[best][j];
2357 break;
2358 } else if (!neg)
2359 neg = ind.order.le[best][j];
2360 }
2361
2362 if (!second && !neg) {
2363 assert(rat);
2364 ind.order.unset_le(best, rat);
2365 ind.expand_rational_vertex(rat);
2366 continue;
2367 }
2368
2369 if (!second)
2370 second = neg;
2371
2372 ind.order.unset_le(best, second);
2373 }
2374
2375 if (!second)
2376 second = neg;
2377
2378 unsigned dim = best->den.NumCols();
2379 temp_evalue diff;
2380 order_sign sign;
2381 for (int k = 0; k < dim; ++k) {
2382 diff = ediff(best->vertex[k], second->vertex[k]);
2383 sign = evalue_sign(diff, ind.D, ind.options->verify->barvinok);
2384
2385 /* neg can never be smaller than best, unless it may still cancel.
2386 * This can happen if positive terms have been determined to
2387 * be equal or less than or equal to some negative term.
2388 */
2389 if (second == neg && !neg_eq && !neg_le) {
2390 if (sign == order_ge)
2391 sign = order_eq;
2392 if (sign == order_unknown)
2393 sign = order_le;
2394 }
2395
2396 if (sign != order_eq)
2397 break;
2398 if (!EVALUE_IS_ZERO(*diff)) {
2399 ind.add_substitution(diff);
2400 ind.perform_pending_substitutions();
2401 }
2402 }
2403 if (sign == order_eq) {
2404 ind.order.set_equal(best, second);
2405 continue;
2406 }
2407 if (sign == order_lt) {
2408 ind.order.lt[best].push_back(second);
2409 ind.order.inc_pred(second);
2410 continue;
2411 }
2412 if (sign == order_gt) {
2413 ind.order.lt[second].push_back(best);
2414 ind.order.inc_pred(best);
2415 continue;
2416 }
2417
2418 split sp(diff, sign == order_le ? split::le :
2419 sign == order_ge ? split::ge : split::lge);
2420
2421 EDomain *Dlt, *Deq, *Dgt;
2422 split_on(sp, ind.D, &Dlt, &Deq, &Dgt, ind.options);
2423 if (ind.options->emptiness_check != BV_LEXMIN_EMPTINESS_CHECK_NONE)
2424 assert(Dlt || Deq || Dgt);
2425 else if (!(Dlt || Deq || Dgt))
2426 /* Must have been empty all along */
2427 break;
2428
2429 if (Deq && (Dlt || Dgt)) {
2430 int locsize = loc.size();
2431 loc.push_back(0);
2432 indicator indeq(ind, Deq);
2433 Deq = NULL;
2434 indeq.add_substitution(diff);
2435 indeq.perform_pending_substitutions();
2436 vector<max_term*> maxeq = lexmin(indeq, nparam, loc);
2437 maxima.insert(maxima.end(), maxeq.begin(), maxeq.end());
2438 loc.resize(locsize);
2439 }
2440 if (Dgt && Dlt) {
2441 int locsize = loc.size();
2442 loc.push_back(1);
2443 indicator indgt(ind, Dgt);
2444 Dgt = NULL;
2445 /* we don't know the new location of these terms in indgt */
2446 /*
2447 indgt.order.lt[second].push_back(best);
2448 indgt.order.inc_pred(best);
2449 */
2450 vector<max_term*> maxgt = lexmin(indgt, nparam, loc);
2451 maxima.insert(maxima.end(), maxgt.begin(), maxgt.end());
2452 loc.resize(locsize);
2453 }
2454
2455 if (Deq) {
2456 loc.push_back(0);
2457 ind.set_domain(Deq);
2458 ind.add_substitution(diff);
2459 ind.perform_pending_substitutions();
2460 }
2461 if (Dlt) {
2462 loc.push_back(-1);
2463 ind.set_domain(Dlt);
2464 ind.order.lt[best].push_back(second);
2465 ind.order.inc_pred(second);
2466 }
2467 if (Dgt) {
2468 loc.push_back(1);
2469 ind.set_domain(Dgt);
2470 ind.order.lt[second].push_back(best);
2471 ind.order.inc_pred(best);
2472 }
2473 } while(1);
2474
2475 return maxima;
2476 }
2477
2478 static void lexmin_base(Polyhedron *P, Polyhedron *C,
2479 Matrix *CP, Matrix *T,
2480 vector<max_term*>& all_max,
2481 lexmin_options *options)
2482 {
2483 unsigned nparam = C->Dimension;
2484 Param_Polyhedron *PP = NULL;
2485
2486 PP = Polyhedron2Param_Polyhedron(P, C, options->verify->barvinok);
2487
2488 unsigned dim = P->Dimension - nparam;
2489
2490 int nd = -1;
2491
2492 indicator_constructor ic(dim, PP, T);
2493
2494 vector<indicator_term *> all_vertices;
2495 construct_rational_vertices(PP, T, T ? T->NbRows-nparam-1 : dim,
2496 nparam, all_vertices);
2497
2498 Polyhedron *TC = true_context(P, C, options->verify->barvinok->MaxRays);
2499 FORALL_REDUCED_DOMAIN(PP, TC, nd, options->verify->barvinok, i, D, rVD)
2500 Param_Vertices *V;
2501
2502 EDomain *epVD = new EDomain(rVD);
2503 indicator ind(ic, D, epVD, options);
2504
2505 FORALL_PVertex_in_ParamPolyhedron(V,D,PP) // _i is internal counter
2506 ind.add(all_vertices[_i]);
2507 END_FORALL_PVertex_in_ParamPolyhedron;
2508
2509 ind.remove_initial_rational_vertices();
2510
2511 vector<int> loc;
2512 vector<max_term*> maxima = lexmin(ind, nparam, loc);
2513 if (CP)
2514 for (int j = 0; j < maxima.size(); ++j)
2515 maxima[j]->substitute(CP, options->verify->barvinok);
2516 all_max.insert(all_max.end(), maxima.begin(), maxima.end());
2517
2518 Domain_Free(rVD);
2519 END_FORALL_REDUCED_DOMAIN
2520 Polyhedron_Free(TC);
2521 for (int i = 0; i < all_vertices.size(); ++i)
2522 delete all_vertices[i];
2523 Param_Polyhedron_Free(PP);
2524 }
2525
2526 static vector<max_term*> lexmin(Polyhedron *P, Polyhedron *C,
2527 lexmin_options *options)
2528 {
2529 unsigned nparam = C->Dimension;
2530 Matrix *T = NULL, *CP = NULL;
2531 Polyhedron *Porig = P;
2532 Polyhedron *Corig = C;
2533 vector<max_term*> all_max;
2534
2535 if (emptyQ2(P))
2536 return all_max;
2537
2538 POL_ENSURE_VERTICES(P);
2539
2540 if (emptyQ2(P))
2541 return all_max;
2542
2543 assert(P->NbBid == 0);
2544
2545 if (P->NbEq > 0)
2546 remove_all_equalities(&P, &C, &CP, &T, nparam,
2547 options->verify->barvinok->MaxRays);
2548 if (!emptyQ2(P))
2549 lexmin_base(P, C, CP, T, all_max, options);
2550 if (CP)
2551 Matrix_Free(CP);
2552 if (T)
2553 Matrix_Free(T);
2554 if (C != Corig)
2555 Polyhedron_Free(C);
2556 if (P != Porig)
2557 Polyhedron_Free(P);
2558
2559 return all_max;
2560 }
2561
2562 static void verify_results(Polyhedron *A, Polyhedron *C,
2563 vector<max_term*>& maxima,
2564 struct verify_options *options);
2565
2566 /* Turn the set dimensions of "context" into parameters and return
2567 * the corresponding parameter domain.
2568 */
2569 static struct isl_basic_set *to_parameter_domain(struct isl_basic_set *context)
2570 {
2571 context = isl_basic_set_move_dims(context, isl_dim_param, 0,
2572 isl_dim_set, 0, isl_basic_set_dim(context, isl_dim_set));
2573 context = isl_basic_set_params(context);
2574 return context;
2575 }
2576
2577 int main(int argc, char **argv)
2578 {
2579 isl_ctx *ctx;
2580 isl_basic_set *context, *bset;
2581 Polyhedron *A, *C;
2582 int neg_one, n;
2583 char s[1024];
2584 int urs_parms = 0;
2585 int urs_unknowns = 0;
2586 int print_solution = 1;
2587 struct lexmin_options *options = lexmin_options_new_with_defaults();
2588 int nparam;
2589 options->verify->barvinok->lookup_table = 0;
2590
2591 argc = lexmin_options_parse(options, argc, argv, ISL_ARG_ALL);
2592 ctx = isl_ctx_alloc_with_options(&lexmin_options_args, options);
2593
2594 context = isl_basic_set_read_from_file(ctx, stdin);
2595 assert(context);
2596 n = fscanf(stdin, "%d", &neg_one);
2597 assert(n == 1);
2598 assert(neg_one == -1);
2599 bset = isl_basic_set_read_from_file(ctx, stdin);
2600
2601 while (fgets(s, sizeof(s), stdin)) {
2602 if (strncasecmp(s, "Maximize", 8) == 0) {
2603 fprintf(stderr, "Maximize option not supported\n");
2604 abort();
2605 }
2606 if (strncasecmp(s, "Rational", 8) == 0) {
2607 fprintf(stderr, "Rational option not supported\n");
2608 abort();
2609 }
2610 if (strncasecmp(s, "Urs_parms", 9) == 0)
2611 urs_parms = 1;
2612 if (strncasecmp(s, "Urs_unknowns", 12) == 0)
2613 urs_unknowns = 1;
2614 }
2615 if (!urs_parms)
2616 context = isl_basic_set_intersect(context,
2617 isl_basic_set_positive_orthant(isl_basic_set_get_space(context)));
2618 context = to_parameter_domain(context);
2619 nparam = isl_basic_set_dim(context, isl_dim_param);
2620 if (nparam != isl_basic_set_dim(bset, isl_dim_param)) {
2621 int dim = isl_basic_set_dim(bset, isl_dim_set);
2622 bset = isl_basic_set_move_dims(bset, isl_dim_param, 0,
2623 isl_dim_set, dim - nparam, nparam);
2624 }
2625 if (!urs_unknowns)
2626 bset = isl_basic_set_intersect(bset,
2627 isl_basic_set_positive_orthant(isl_basic_set_get_space(bset)));
2628
2629 if (options->verify->verify)
2630 print_solution = 0;
2631
2632 A = isl_basic_set_to_polylib(bset);
2633 verify_options_set_range(options->verify, A->Dimension);
2634 C = isl_basic_set_to_polylib(context);
2635 vector<max_term*> maxima = lexmin(A, C, options);
2636 if (print_solution) {
2637 char **param_names;
2638 param_names = util_generate_names(C->Dimension, "p");
2639 for (int i = 0; i < maxima.size(); ++i)
2640 maxima[i]->print(cout, param_names,
2641 options->verify->barvinok);
2642 util_free_names(C->Dimension, param_names);
2643 }
2644
2645 if (options->verify->verify)
2646 verify_results(A, C, maxima, options->verify);
2647
2648 for (int i = 0; i < maxima.size(); ++i)
2649 delete maxima[i];
2650
2651 Polyhedron_Free(A);
2652 Polyhedron_Free(C);
2653
2654 isl_basic_set_free(bset);
2655 isl_basic_set_free(context);
2656 isl_ctx_free(ctx);
2657
2658 return 0;
2659 }
2660
2661 static bool lexmin(int pos, Polyhedron *P, Value *context)
2662 {
2663 Value LB, UB, k;
2664 int lu_flags;
2665 bool found = false;
2666
2667 if (emptyQ(P))
2668 return false;
2669
2670 value_init(LB); value_init(UB); value_init(k);
2671 value_set_si(LB,0);
2672 value_set_si(UB,0);
2673 lu_flags = lower_upper_bounds(pos,P,context,&LB,&UB);
2674 assert(!(lu_flags & LB_INFINITY));
2675
2676 value_set_si(context[pos],0);
2677 if (!lu_flags && value_lt(UB,LB)) {
2678 value_clear(LB); value_clear(UB); value_clear(k);
2679 return false;
2680 }
2681 if (!P->next) {
2682 value_assign(context[pos], LB);
2683 value_clear(LB); value_clear(UB); value_clear(k);
2684 return true;
2685 }
2686 for (value_assign(k,LB); lu_flags || value_le(k,UB); value_increment(k,k)) {
2687 value_assign(context[pos],k);
2688 if ((found = lexmin(pos+1, P->next, context)))
2689 break;
2690 }
2691 if (!found)
2692 value_set_si(context[pos],0);
2693 value_clear(LB); value_clear(UB); value_clear(k);
2694 return found;
2695 }
2696
2697 static void print_list(FILE *out, Value *z, const char* brackets, int len)
2698 {
2699 fprintf(out, "%c", brackets[0]);
2700 value_print(out, VALUE_FMT,z[0]);
2701 for (int k = 1; k < len; ++k) {
2702 fprintf(out, ", ");
2703 value_print(out, VALUE_FMT,z[k]);
2704 }
2705 fprintf(out, "%c", brackets[1]);
2706 }
2707
2708 static int check_poly_lexmin(const struct check_poly_data *data,
2709 int nparam, Value *z,
2710 const struct verify_options *options);
2711
2712 struct check_poly_lexmin_data : public check_poly_data {
2713 Polyhedron *S;
2714 vector<max_term*>& maxima;
2715
2716 check_poly_lexmin_data(Polyhedron *S, Value *z,
2717 vector<max_term*>& maxima) : S(S), maxima(maxima) {
2718 this->z = z;
2719 this->check = &check_poly_lexmin;
2720 }
2721 };
2722
2723 static int check_poly_lexmin(const struct check_poly_data *data,
2724 int nparam, Value *z,
2725 const struct verify_options *options)
2726 {
2727 const check_poly_lexmin_data *lexmin_data;
2728 lexmin_data = static_cast<const check_poly_lexmin_data *>(data);
2729 Polyhedron *S = lexmin_data->S;
2730 vector<max_term*>& maxima = lexmin_data->maxima;
2731 int k;
2732 bool found = lexmin(1, S, lexmin_data->z);
2733
2734 if (options->print_all) {
2735 printf("lexmin");
2736 print_list(stdout, z, "()", nparam);
2737 printf(" = ");
2738 if (found)
2739 print_list(stdout, lexmin_data->z+1, "[]", S->Dimension-nparam);
2740 printf(" ");
2741 }
2742
2743 Vector *min = NULL;
2744 for (int i = 0; i < maxima.size(); ++i)
2745 if ((min = maxima[i]->eval(z, options->barvinok->MaxRays)))
2746 break;
2747
2748 int ok = !(found ^ !!min);
2749 if (found && min)
2750 for (int i = 0; i < S->Dimension-nparam; ++i)
2751 if (value_ne(lexmin_data->z[1+i], min->p[i])) {
2752 ok = 0;
2753 break;
2754 }
2755 if (!ok) {
2756 printf("\n");
2757 fflush(stdout);
2758 fprintf(stderr, "Error !\n");
2759 fprintf(stderr, "lexmin");
2760 print_list(stderr, z, "()", nparam);
2761 fprintf(stderr, " should be ");
2762 if (found)
2763 print_list(stderr, lexmin_data->z+1, "[]", S->Dimension-nparam);
2764 fprintf(stderr, " while digging gives ");
2765 if (min)
2766 print_list(stderr, min->p, "[]", S->Dimension-nparam);
2767 fprintf(stderr, ".\n");
2768 return 0;
2769 } else if (options->print_all)
2770 printf("OK.\n");
2771 if (min)
2772 Vector_Free(min);
2773
2774 for (k = 1; k <= S->Dimension-nparam; ++k)
2775 value_set_si(lexmin_data->z[k], 0);
2776
2777 return 0;
2778 }
2779
2780 void verify_results(Polyhedron *A, Polyhedron *C, vector<max_term*>& maxima,
2781 struct verify_options *options)
2782 {
2783 Polyhedron *CS, *S;
2784 unsigned nparam = C->Dimension;
2785 unsigned MaxRays = options->barvinok->MaxRays;
2786 Vector *p;
2787
2788 CS = check_poly_context_scan(A, &C, nparam, options);
2789
2790 p = Vector_Alloc(A->Dimension+2);
2791 value_set_si(p->p[A->Dimension+1], 1);
2792
2793 S = Polyhedron_Scan(A, C, MaxRays & POL_NO_DUAL ? 0 : MaxRays);
2794
2795 check_poly_init(C, options);
2796
2797 if (S) {
2798 if (!(CS && emptyQ2(CS))) {
2799 check_poly_lexmin_data data(S, p->p, maxima);
2800 check_poly(CS, &data, nparam, 0, p->p+S->Dimension-nparam+1, options);
2801 }
2802 Domain_Free(S);
2803 }
2804
2805 if (!options->print_all)
2806 printf("\n");
2807
2808 Vector_Free(p);
2809 if (CS) {
2810 Domain_Free(CS);
2811 Polyhedron_Free(C);
2812 }
2813 }
lattice point enumeration library