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