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