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