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