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