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