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