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