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