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