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