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