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