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