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