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