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