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