6 #define partition STL_PARTITION
10 #include <NTL/vec_ZZ.h>
11 #include <NTL/mat_ZZ.h>
12 #include <isl_set_polylib.h>
13 #include <barvinok/barvinok.h>
14 #include <barvinok/evalue.h>
15 #include <barvinok/options.h>
16 #include <barvinok/util.h>
17 #include "conversion.h"
18 #include "decomposer.h"
19 #include "lattice_point.h"
20 #include "reduce_domain.h"
23 #include "evalue_util.h"
24 #include "remove_equalities.h"
28 #include "param_util.h"
30 #undef CS /* for Solaris 10 */
41 #define ALLOC(type) (type*)malloc(sizeof(type))
42 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
44 static int type_offset(enode
*p
)
46 return p
->type
== fractional
? 1 :
47 p
->type
== flooring
? 1 : 0;
50 void compute_evalue(evalue
*e
, Value
*val
, Value
*res
)
52 double d
= compute_evalue(e
, val
);
57 value_set_double(*res
, d
);
60 struct indicator_term
{
62 int pos
; /* number of rational vertex */
63 int n
; /* number of cone associated to given rational vertex */
67 indicator_term(unsigned dim
, int pos
) {
69 vertex
= new evalue
* [dim
];
74 indicator_term(unsigned dim
, int pos
, int n
) {
75 den
.SetDims(dim
, dim
);
76 vertex
= new evalue
* [dim
];
80 indicator_term(const indicator_term
& src
) {
85 unsigned dim
= den
.NumCols();
86 vertex
= new evalue
* [dim
];
87 for (int i
= 0; i
< dim
; ++i
) {
88 vertex
[i
] = ALLOC(evalue
);
89 value_init(vertex
[i
]->d
);
90 evalue_copy(vertex
[i
], src
.vertex
[i
]);
93 void swap(indicator_term
*other
) {
95 tmp
= sign
; sign
= other
->sign
; other
->sign
= tmp
;
96 tmp
= pos
; pos
= other
->pos
; other
->pos
= tmp
;
97 tmp
= n
; n
= other
->n
; other
->n
= tmp
;
98 mat_ZZ tmp_den
= den
; den
= other
->den
; other
->den
= tmp_den
;
99 unsigned dim
= den
.NumCols();
100 for (int i
= 0; i
< dim
; ++i
) {
101 evalue
*tmp
= vertex
[i
];
102 vertex
[i
] = other
->vertex
[i
];
103 other
->vertex
[i
] = tmp
;
107 unsigned dim
= den
.NumCols();
108 for (int i
= 0; i
< dim
; ++i
)
109 evalue_free(vertex
[i
]);
112 void print(ostream
& os
, char **p
) const;
113 void substitute(Matrix
*T
);
115 void substitute(evalue
*fract
, evalue
*val
);
116 void substitute(int pos
, evalue
*val
);
117 void reduce_in_domain(Polyhedron
*D
);
118 bool is_opposite(const indicator_term
*neg
) const;
119 vec_ZZ
eval(Value
*val
) const {
121 unsigned dim
= den
.NumCols();
125 for (int i
= 0; i
< dim
; ++i
) {
126 compute_evalue(vertex
[i
], val
, &tmp
);
134 static int evalue_rational_cmp(const evalue
*e1
, const evalue
*e2
)
142 assert(value_notzero_p(e1
->d
));
143 assert(value_notzero_p(e2
->d
));
144 value_multiply(m
, e1
->x
.n
, e2
->d
);
145 value_multiply(m2
, e2
->x
.n
, e1
->d
);
148 else if (value_gt(m
, m2
))
158 static int evalue_cmp(const evalue
*e1
, const evalue
*e2
)
160 if (value_notzero_p(e1
->d
)) {
161 if (value_zero_p(e2
->d
))
163 return evalue_rational_cmp(e1
, e2
);
165 if (value_notzero_p(e2
->d
))
167 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
168 return e1
->x
.p
->type
- e2
->x
.p
->type
;
169 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
170 return e1
->x
.p
->size
- e2
->x
.p
->size
;
171 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
172 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
173 assert(e1
->x
.p
->type
== polynomial
||
174 e1
->x
.p
->type
== fractional
||
175 e1
->x
.p
->type
== flooring
);
176 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
) {
177 int s
= evalue_cmp(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]);
184 void evalue_length(evalue
*e
, int len
[2])
189 while (value_zero_p(e
->d
)) {
190 assert(e
->x
.p
->type
== polynomial
||
191 e
->x
.p
->type
== fractional
||
192 e
->x
.p
->type
== flooring
);
193 if (e
->x
.p
->type
== polynomial
)
197 int offset
= type_offset(e
->x
.p
);
198 assert(e
->x
.p
->size
== offset
+2);
199 e
= &e
->x
.p
->arr
[offset
];
203 static bool it_smaller(const indicator_term
* it1
, const indicator_term
* it2
)
207 int len1
[2], len2
[2];
208 unsigned dim
= it1
->den
.NumCols();
209 for (int i
= 0; i
< dim
; ++i
) {
210 evalue_length(it1
->vertex
[i
], len1
);
211 evalue_length(it2
->vertex
[i
], len2
);
212 if (len1
[0] != len2
[0])
213 return len1
[0] < len2
[0];
214 if (len1
[1] != len2
[1])
215 return len1
[1] < len2
[1];
217 if (it1
->pos
!= it2
->pos
)
218 return it1
->pos
< it2
->pos
;
219 if (it1
->n
!= it2
->n
)
220 return it1
->n
< it2
->n
;
221 int s
= lex_cmp(it1
->den
, it2
->den
);
224 for (int i
= 0; i
< dim
; ++i
) {
225 s
= evalue_cmp(it1
->vertex
[i
], it2
->vertex
[i
]);
229 assert(it1
->sign
!= 0);
230 assert(it2
->sign
!= 0);
231 if (it1
->sign
!= it2
->sign
)
232 return it1
->sign
> 0;
237 static const int requires_resort
;
238 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
239 return it_smaller(it1
, it2
);
242 const int smaller_it::requires_resort
= 1;
244 struct smaller_it_p
{
245 static const int requires_resort
;
246 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
250 const int smaller_it_p::requires_resort
= 0;
252 /* Returns true if this and neg are opposite using the knowledge
253 * that they have the same numerator.
254 * In particular, we check that the signs are different and that
255 * the denominator is the same.
257 bool indicator_term::is_opposite(const indicator_term
*neg
) const
259 if (sign
+ neg
->sign
!= 0)
266 void indicator_term::reduce_in_domain(Polyhedron
*D
)
268 for (int k
= 0; k
< den
.NumCols(); ++k
) {
269 reduce_evalue_in_domain(vertex
[k
], D
);
270 if (evalue_range_reduction_in_domain(vertex
[k
], D
))
271 reduce_evalue(vertex
[k
]);
275 void indicator_term::print(ostream
& os
, char **p
) const
277 unsigned dim
= den
.NumCols();
278 unsigned factors
= den
.NumRows();
286 for (int i
= 0; i
< dim
; ++i
) {
289 evalue_print(os
, vertex
[i
], p
);
292 for (int i
= 0; i
< factors
; ++i
) {
293 os
<< " + t" << i
<< "*[";
294 for (int j
= 0; j
< dim
; ++j
) {
301 os
<< " ((" << pos
<< ", " << n
<< ", " << (void*)this << "))";
304 /* Perform the substitution specified by T on the variables.
305 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
306 * The substitution is performed as in gen_fun::substitute
308 void indicator_term::substitute(Matrix
*T
)
310 unsigned dim
= den
.NumCols();
311 unsigned nparam
= T
->NbColumns
- dim
- 1;
312 unsigned newdim
= T
->NbRows
- nparam
- 1;
315 matrix2zz(T
, trans
, newdim
, dim
);
316 trans
= transpose(trans
);
318 newvertex
= new evalue
* [newdim
];
321 v
.SetLength(nparam
+1);
324 value_init(factor
.d
);
325 value_set_si(factor
.d
, 1);
326 value_init(factor
.x
.n
);
327 for (int i
= 0; i
< newdim
; ++i
) {
328 values2zz(T
->p
[i
]+dim
, v
, nparam
+1);
329 newvertex
[i
] = multi_monom(v
);
331 for (int j
= 0; j
< dim
; ++j
) {
332 if (value_zero_p(T
->p
[i
][j
]))
336 evalue_copy(&term
, vertex
[j
]);
337 value_assign(factor
.x
.n
, T
->p
[i
][j
]);
338 emul(&factor
, &term
);
339 eadd(&term
, newvertex
[i
]);
340 free_evalue_refs(&term
);
343 free_evalue_refs(&factor
);
344 for (int i
= 0; i
< dim
; ++i
)
345 evalue_free(vertex
[i
]);
350 static void evalue_add_constant(evalue
*e
, ZZ v
)
355 /* go down to constant term */
356 while (value_zero_p(e
->d
))
357 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
360 value_multiply(tmp
, tmp
, e
->d
);
361 value_addto(e
->x
.n
, e
->x
.n
, tmp
);
366 /* Make all powers in denominator lexico-positive */
367 void indicator_term::normalize()
370 extra_vertex
.SetLength(den
.NumCols());
371 for (int r
= 0; r
< den
.NumRows(); ++r
) {
372 for (int k
= 0; k
< den
.NumCols(); ++k
) {
379 extra_vertex
+= den
[r
];
383 for (int k
= 0; k
< extra_vertex
.length(); ++k
)
384 if (extra_vertex
[k
] != 0)
385 evalue_add_constant(vertex
[k
], extra_vertex
[k
]);
388 static void substitute(evalue
*e
, evalue
*fract
, evalue
*val
)
392 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
393 if (t
->x
.p
->type
== fractional
&& eequal(&t
->x
.p
->arr
[0], fract
))
396 if (value_notzero_p(t
->d
))
399 free_evalue_refs(&t
->x
.p
->arr
[0]);
400 evalue
*term
= &t
->x
.p
->arr
[2];
407 free_evalue_refs(term
);
413 void indicator_term::substitute(evalue
*fract
, evalue
*val
)
415 unsigned dim
= den
.NumCols();
416 for (int i
= 0; i
< dim
; ++i
) {
417 ::substitute(vertex
[i
], fract
, val
);
421 static void substitute(evalue
*e
, int pos
, evalue
*val
)
425 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
426 if (t
->x
.p
->type
== polynomial
&& t
->x
.p
->pos
== pos
)
429 if (value_notzero_p(t
->d
))
432 evalue
*term
= &t
->x
.p
->arr
[1];
439 free_evalue_refs(term
);
445 void indicator_term::substitute(int pos
, evalue
*val
)
447 unsigned dim
= den
.NumCols();
448 for (int i
= 0; i
< dim
; ++i
) {
449 ::substitute(vertex
[i
], pos
, val
);
453 struct indicator_constructor
: public signed_cone_consumer
,
454 public vertex_decomposer
{
456 vector
<indicator_term
*> *terms
;
457 Matrix
*T
; /* Transformation to original space */
462 indicator_constructor(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
,
464 vertex_decomposer(PP
, *this), T(T
), nbV(PP
->nbV
) {
465 vertex
.SetLength(dim
);
466 terms
= new vector
<indicator_term
*>[PP
->nbV
];
468 ~indicator_constructor() {
469 for (int i
= 0; i
< nbV
; ++i
)
470 for (int j
= 0; j
< terms
[i
].size(); ++j
)
475 void print(ostream
& os
, char **p
);
477 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
478 void decompose_at_vertex(Param_Vertices
*V
, int _i
,
479 barvinok_options
*options
) {
482 vertex_decomposer::decompose_at_vertex(V
, _i
, options
);
486 void indicator_constructor::handle(const signed_cone
& sc
, barvinok_options
*options
)
489 unsigned dim
= vertex
.length();
491 assert(sc
.rays
.NumRows() == dim
);
493 indicator_term
*term
= new indicator_term(dim
, pos
, n
++);
494 term
->sign
= sc
.sign
;
495 terms
[vert
].push_back(term
);
497 lattice_point(V
, sc
.rays
, vertex
, term
->vertex
, options
);
500 for (int r
= 0; r
< dim
; ++r
) {
501 for (int j
= 0; j
< dim
; ++j
) {
502 if (term
->den
[r
][j
] == 0)
504 if (term
->den
[r
][j
] > 0)
506 term
->sign
= -term
->sign
;
507 term
->den
[r
] = -term
->den
[r
];
508 vertex
+= term
->den
[r
];
513 for (int i
= 0; i
< dim
; ++i
) {
514 if (!term
->vertex
[i
]) {
515 term
->vertex
[i
] = ALLOC(evalue
);
516 value_init(term
->vertex
[i
]->d
);
517 value_init(term
->vertex
[i
]->x
.n
);
518 zz2value(vertex
[i
], term
->vertex
[i
]->x
.n
);
519 value_set_si(term
->vertex
[i
]->d
, 1);
524 evalue_add_constant(term
->vertex
[i
], vertex
[i
]);
532 lex_order_rows(term
->den
);
535 void indicator_constructor::print(ostream
& os
, char **p
)
537 for (int i
= 0; i
< PP
->nbV
; ++i
)
538 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
539 os
<< "i: " << i
<< ", j: " << j
<< endl
;
540 terms
[i
][j
]->print(os
, p
);
545 void indicator_constructor::normalize()
547 for (int i
= 0; i
< PP
->nbV
; ++i
)
548 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
550 vertex
.SetLength(terms
[i
][j
]->den
.NumCols());
551 for (int r
= 0; r
< terms
[i
][j
]->den
.NumRows(); ++r
) {
552 for (int k
= 0; k
< terms
[i
][j
]->den
.NumCols(); ++k
) {
553 if (terms
[i
][j
]->den
[r
][k
] == 0)
555 if (terms
[i
][j
]->den
[r
][k
] > 0)
557 terms
[i
][j
]->sign
= -terms
[i
][j
]->sign
;
558 terms
[i
][j
]->den
[r
] = -terms
[i
][j
]->den
[r
];
559 vertex
+= terms
[i
][j
]->den
[r
];
563 lex_order_rows(terms
[i
][j
]->den
);
564 for (int k
= 0; k
< vertex
.length(); ++k
)
566 evalue_add_constant(terms
[i
][j
]->vertex
[k
], vertex
[k
]);
570 struct order_cache_el
{
572 order_cache_el
copy() const {
574 for (int i
= 0; i
< e
.size(); ++i
) {
575 evalue
*c
= new evalue
;
577 evalue_copy(c
, e
[i
]);
583 for (int i
= 0; i
< e
.size(); ++i
) {
584 free_evalue_refs(e
[i
]);
591 evalue_set_si(&mone
, -1, 1);
592 for (int i
= 0; i
< e
.size(); ++i
)
594 free_evalue_refs(&mone
);
596 void print(ostream
& os
, char **p
);
599 void order_cache_el::print(ostream
& os
, char **p
)
602 for (int i
= 0; i
< e
.size(); ++i
) {
605 evalue_print(os
, e
[i
], p
);
611 vector
<order_cache_el
> lt
;
612 vector
<order_cache_el
> le
;
613 vector
<order_cache_el
> unknown
;
615 void clear_transients() {
616 for (int i
= 0; i
< le
.size(); ++i
)
618 for (int i
= 0; i
< unknown
.size(); ++i
)
625 for (int i
= 0; i
< lt
.size(); ++i
)
629 void add(order_cache_el
& cache_el
, order_sign sign
);
630 order_sign
check_lt(vector
<order_cache_el
>* list
,
631 const indicator_term
*a
, const indicator_term
*b
,
632 order_cache_el
& cache_el
);
633 order_sign
check_lt(const indicator_term
*a
, const indicator_term
*b
,
634 order_cache_el
& cache_el
);
635 order_sign
check_direct(const indicator_term
*a
, const indicator_term
*b
,
636 order_cache_el
& cache_el
);
637 order_sign
check(const indicator_term
*a
, const indicator_term
*b
,
638 order_cache_el
& cache_el
);
639 void copy(const order_cache
& cache
);
640 void print(ostream
& os
, char **p
);
643 void order_cache::copy(const order_cache
& cache
)
645 for (int i
= 0; i
< cache
.lt
.size(); ++i
) {
646 order_cache_el n
= cache
.lt
[i
].copy();
651 void order_cache::add(order_cache_el
& cache_el
, order_sign sign
)
653 if (sign
== order_lt
) {
654 lt
.push_back(cache_el
);
655 } else if (sign
== order_gt
) {
657 lt
.push_back(cache_el
);
658 } else if (sign
== order_le
) {
659 le
.push_back(cache_el
);
660 } else if (sign
== order_ge
) {
662 le
.push_back(cache_el
);
663 } else if (sign
== order_unknown
) {
664 unknown
.push_back(cache_el
);
666 assert(sign
== order_eq
);
673 static evalue
*ediff(const evalue
*a
, const evalue
*b
)
677 evalue_set_si(&mone
, -1, 1);
678 evalue
*diff
= new evalue
;
680 evalue_copy(diff
, b
);
684 free_evalue_refs(&mone
);
688 static bool evalue_first_difference(const evalue
*e1
, const evalue
*e2
,
689 const evalue
**d1
, const evalue
**d2
)
694 if (value_ne(e1
->d
, e2
->d
))
697 if (value_notzero_p(e1
->d
)) {
698 if (value_eq(e1
->x
.n
, e2
->x
.n
))
702 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
704 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
706 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
709 assert(e1
->x
.p
->type
== polynomial
||
710 e1
->x
.p
->type
== fractional
||
711 e1
->x
.p
->type
== flooring
);
712 int offset
= type_offset(e1
->x
.p
);
713 assert(e1
->x
.p
->size
== offset
+2);
714 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
)
715 if (i
!= type_offset(e1
->x
.p
) &&
716 !eequal(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]))
719 return evalue_first_difference(&e1
->x
.p
->arr
[offset
],
720 &e2
->x
.p
->arr
[offset
], d1
, d2
);
723 static order_sign
evalue_diff_constant_sign(const evalue
*e1
, const evalue
*e2
)
725 if (!evalue_first_difference(e1
, e2
, &e1
, &e2
))
727 if (value_zero_p(e1
->d
) || value_zero_p(e2
->d
))
728 return order_undefined
;
729 int s
= evalue_rational_cmp(e1
, e2
);
738 order_sign
order_cache::check_lt(vector
<order_cache_el
>* list
,
739 const indicator_term
*a
, const indicator_term
*b
,
740 order_cache_el
& cache_el
)
742 order_sign sign
= order_undefined
;
743 for (int i
= 0; i
< list
->size(); ++i
) {
745 for (j
= cache_el
.e
.size(); j
< (*list
)[i
].e
.size(); ++j
)
746 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
747 for (j
= 0; j
< (*list
)[i
].e
.size(); ++j
) {
748 order_sign diff_sign
;
749 diff_sign
= evalue_diff_constant_sign((*list
)[i
].e
[j
], cache_el
.e
[j
]);
750 if (diff_sign
== order_gt
) {
753 } else if (diff_sign
== order_lt
)
755 else if (diff_sign
== order_undefined
)
758 assert(diff_sign
== order_eq
);
760 if (j
== (*list
)[i
].e
.size())
761 sign
= list
== <
? order_lt
: order_le
;
762 if (sign
!= order_undefined
)
768 order_sign
order_cache::check_direct(const indicator_term
*a
,
769 const indicator_term
*b
,
770 order_cache_el
& cache_el
)
772 order_sign sign
= check_lt(<
, a
, b
, cache_el
);
773 if (sign
!= order_undefined
)
775 sign
= check_lt(&le
, a
, b
, cache_el
);
776 if (sign
!= order_undefined
)
779 for (int i
= 0; i
< unknown
.size(); ++i
) {
781 for (j
= cache_el
.e
.size(); j
< unknown
[i
].e
.size(); ++j
)
782 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
783 for (j
= 0; j
< unknown
[i
].e
.size(); ++j
) {
784 if (!eequal(unknown
[i
].e
[j
], cache_el
.e
[j
]))
787 if (j
== unknown
[i
].e
.size()) {
788 sign
= order_unknown
;
795 order_sign
order_cache::check(const indicator_term
*a
, const indicator_term
*b
,
796 order_cache_el
& cache_el
)
798 order_sign sign
= check_direct(a
, b
, cache_el
);
799 if (sign
!= order_undefined
)
801 int size
= cache_el
.e
.size();
803 sign
= check_direct(a
, b
, cache_el
);
805 assert(cache_el
.e
.size() == size
);
806 if (sign
== order_undefined
)
808 if (sign
== order_lt
)
810 else if (sign
== order_le
)
813 assert(sign
== order_unknown
);
819 struct partial_order
{
822 typedef std::set
<const indicator_term
*, smaller_it
> head_type
;
824 typedef map
<const indicator_term
*, int, smaller_it
> pred_type
;
826 typedef map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> order_type
;
835 partial_order(indicator
*ind
) : ind(ind
) {}
836 void copy(const partial_order
& order
,
837 map
< const indicator_term
*, indicator_term
* > old2new
);
839 order_type::iterator i
;
840 pred_type::iterator j
;
841 head_type::iterator k
;
843 if (head
.key_comp().requires_resort
) {
845 for (k
= head
.begin(); k
!= head
.end(); ++k
)
851 if (pred
.key_comp().requires_resort
) {
853 for (j
= pred
.begin(); j
!= pred
.end(); ++j
)
854 new_pred
[(*j
).first
] = (*j
).second
;
859 if (lt
.key_comp().requires_resort
) {
861 for (i
= lt
.begin(); i
!= lt
.end(); ++i
)
862 m
[(*i
).first
] = (*i
).second
;
867 if (le
.key_comp().requires_resort
) {
869 for (i
= le
.begin(); i
!= le
.end(); ++i
)
870 m
[(*i
).first
] = (*i
).second
;
875 if (eq
.key_comp().requires_resort
) {
877 for (i
= eq
.begin(); i
!= eq
.end(); ++i
)
878 m
[(*i
).first
] = (*i
).second
;
883 if (pending
.key_comp().requires_resort
) {
885 for (i
= pending
.begin(); i
!= pending
.end(); ++i
)
886 m
[(*i
).first
] = (*i
).second
;
892 order_sign
compare(const indicator_term
*a
, const indicator_term
*b
);
893 void set_equal(const indicator_term
*a
, const indicator_term
*b
);
894 void unset_le(const indicator_term
*a
, const indicator_term
*b
);
895 void dec_pred(const indicator_term
*it
) {
896 if (--pred
[it
] == 0) {
901 void inc_pred(const indicator_term
*it
) {
902 if (head
.find(it
) != head
.end())
907 bool compared(const indicator_term
* a
, const indicator_term
* b
);
908 void add(const indicator_term
* it
, std::set
<const indicator_term
*> *filter
);
909 void remove(const indicator_term
* it
);
911 void print(ostream
& os
, char **p
);
913 /* replace references to orig to references to replacement */
914 void replace(const indicator_term
* orig
, indicator_term
* replacement
);
915 void sanity_check() const;
918 /* We actually replace the contents of orig by that of replacement,
919 * but we have to be careful since replacing the content changes
920 * the order of orig in the maps.
922 void partial_order::replace(const indicator_term
* orig
, indicator_term
* replacement
)
924 head_type::iterator k
;
926 bool is_head
= k
!= head
.end();
931 orig_pred
= pred
[orig
];
934 vector
<const indicator_term
* > orig_lt
;
935 vector
<const indicator_term
* > orig_le
;
936 vector
<const indicator_term
* > orig_eq
;
937 vector
<const indicator_term
* > orig_pending
;
938 order_type::iterator i
;
939 bool in_lt
= ((i
= lt
.find(orig
)) != lt
.end());
941 orig_lt
= (*i
).second
;
944 bool in_le
= ((i
= le
.find(orig
)) != le
.end());
946 orig_le
= (*i
).second
;
949 bool in_eq
= ((i
= eq
.find(orig
)) != eq
.end());
951 orig_eq
= (*i
).second
;
954 bool in_pending
= ((i
= pending
.find(orig
)) != pending
.end());
956 orig_pending
= (*i
).second
;
959 indicator_term
*old
= const_cast<indicator_term
*>(orig
);
960 old
->swap(replacement
);
964 pred
[old
] = orig_pred
;
972 pending
[old
] = orig_pending
;
975 void partial_order::unset_le(const indicator_term
*a
, const indicator_term
*b
)
977 vector
<const indicator_term
*>::iterator i
;
978 i
= std::find(le
[a
].begin(), le
[a
].end(), b
);
980 if (le
[a
].size() == 0)
983 i
= std::find(pending
[a
].begin(), pending
[a
].end(), b
);
984 if (i
!= pending
[a
].end())
988 void partial_order::set_equal(const indicator_term
*a
, const indicator_term
*b
)
990 if (eq
[a
].size() == 0)
992 if (eq
[b
].size() == 0)
997 if (pred
.key_comp()(b
, a
)) {
998 const indicator_term
*c
= a
;
1003 const indicator_term
*base
= a
;
1005 order_type::iterator i
;
1007 for (int j
= 0; j
< eq
[b
].size(); ++j
) {
1008 eq
[base
].push_back(eq
[b
][j
]);
1009 eq
[eq
[b
][j
]][0] = base
;
1014 if (i
!= lt
.end()) {
1015 for (int j
= 0; j
< lt
[b
].size(); ++j
) {
1016 if (std::find(eq
[base
].begin(), eq
[base
].end(), lt
[b
][j
]) != eq
[base
].end())
1018 else if (std::find(lt
[base
].begin(), lt
[base
].end(), lt
[b
][j
])
1022 lt
[base
].push_back(lt
[b
][j
]);
1028 if (i
!= le
.end()) {
1029 for (int j
= 0; j
< le
[b
].size(); ++j
) {
1030 if (std::find(eq
[base
].begin(), eq
[base
].end(), le
[b
][j
]) != eq
[base
].end())
1032 else if (std::find(le
[base
].begin(), le
[base
].end(), le
[b
][j
])
1036 le
[base
].push_back(le
[b
][j
]);
1041 i
= pending
.find(base
);
1042 if (i
!= pending
.end()) {
1043 vector
<const indicator_term
* > old
= pending
[base
];
1044 pending
[base
].clear();
1045 for (int j
= 0; j
< old
.size(); ++j
) {
1046 if (std::find(eq
[base
].begin(), eq
[base
].end(), old
[j
]) == eq
[base
].end())
1047 pending
[base
].push_back(old
[j
]);
1051 i
= pending
.find(b
);
1052 if (i
!= pending
.end()) {
1053 for (int j
= 0; j
< pending
[b
].size(); ++j
) {
1054 if (std::find(eq
[base
].begin(), eq
[base
].end(), pending
[b
][j
]) == eq
[base
].end())
1055 pending
[base
].push_back(pending
[b
][j
]);
1061 void partial_order::copy(const partial_order
& order
,
1062 map
< const indicator_term
*, indicator_term
* > old2new
)
1064 cache
.copy(order
.cache
);
1066 order_type::const_iterator i
;
1067 pred_type::const_iterator j
;
1068 head_type::const_iterator k
;
1070 for (k
= order
.head
.begin(); k
!= order
.head
.end(); ++k
)
1071 head
.insert(old2new
[*k
]);
1073 for (j
= order
.pred
.begin(); j
!= order
.pred
.end(); ++j
)
1074 pred
[old2new
[(*j
).first
]] = (*j
).second
;
1076 for (i
= order
.lt
.begin(); i
!= order
.lt
.end(); ++i
) {
1077 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1078 lt
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1080 for (i
= order
.le
.begin(); i
!= order
.le
.end(); ++i
) {
1081 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1082 le
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1084 for (i
= order
.eq
.begin(); i
!= order
.eq
.end(); ++i
) {
1085 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1086 eq
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1088 for (i
= order
.pending
.begin(); i
!= order
.pending
.end(); ++i
) {
1089 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1090 pending
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1096 vector
<evalue
*> max
;
1098 void print(ostream
& os
, char **p
, barvinok_options
*options
) const;
1099 void substitute(Matrix
*T
, barvinok_options
*options
);
1100 Vector
*eval(Value
*val
, unsigned MaxRays
) const;
1103 for (int i
= 0; i
< max
.size(); ++i
) {
1104 free_evalue_refs(max
[i
]);
1112 * Project on first dim dimensions
1114 Polyhedron
* Polyhedron_Project_Initial(Polyhedron
*P
, int dim
)
1120 if (P
->Dimension
== dim
)
1121 return Polyhedron_Copy(P
);
1123 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
1124 for (i
= 0; i
< dim
; ++i
)
1125 value_set_si(T
->p
[i
][i
], 1);
1126 value_set_si(T
->p
[dim
][P
->Dimension
], 1);
1127 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
1133 vector
<indicator_term
*> term
;
1134 indicator_constructor
& ic
;
1135 partial_order order
;
1139 lexmin_options
*options
;
1140 vector
<evalue
*> substitutions
;
1142 indicator(indicator_constructor
& ic
, Param_Domain
*PD
, EDomain
*D
,
1143 lexmin_options
*options
) :
1144 ic(ic
), PD(PD
), D(D
), order(this), options(options
), P(NULL
) {}
1145 indicator(const indicator
& ind
, EDomain
*D
) :
1146 ic(ind
.ic
), PD(ind
.PD
), D(NULL
), order(this), options(ind
.options
),
1147 P(Polyhedron_Copy(ind
.P
)) {
1148 map
< const indicator_term
*, indicator_term
* > old2new
;
1149 for (int i
= 0; i
< ind
.term
.size(); ++i
) {
1150 indicator_term
*it
= new indicator_term(*ind
.term
[i
]);
1151 old2new
[ind
.term
[i
]] = it
;
1154 order
.copy(ind
.order
, old2new
);
1158 for (int i
= 0; i
< term
.size(); ++i
)
1166 void set_domain(EDomain
*D
) {
1167 order
.cache
.clear_transients();
1171 int nparam
= ic
.PP
->Constraints
->NbColumns
-2 - ic
.vertex
.length();
1172 if (options
->reduce
) {
1173 Polyhedron
*Q
= Polyhedron_Project_Initial(D
->D
, nparam
);
1174 Q
= DomainConstraintSimplify(Q
, options
->verify
->barvinok
->MaxRays
);
1175 if (!P
|| !PolyhedronIncludes(Q
, P
))
1176 reduce_in_domain(Q
);
1184 void add(const indicator_term
* it
);
1185 void remove(const indicator_term
* it
);
1186 void remove_initial_rational_vertices();
1187 void expand_rational_vertex(const indicator_term
*initial
);
1189 void print(ostream
& os
, char **p
);
1191 void peel(int i
, int j
);
1192 void combine(const indicator_term
*a
, const indicator_term
*b
);
1193 void add_substitution(evalue
*equation
);
1194 void perform_pending_substitutions();
1195 void reduce_in_domain(Polyhedron
*D
);
1196 bool handle_equal_numerators(const indicator_term
*base
);
1198 max_term
* create_max_term(const indicator_term
*it
);
1200 void substitute(evalue
*equation
);
1203 void partial_order::sanity_check() const
1205 order_type::const_iterator i
;
1206 order_type::const_iterator prev
;
1207 order_type::const_iterator l
;
1208 pred_type::const_iterator k
, prev_k
;
1210 for (k
= pred
.begin(); k
!= pred
.end(); prev_k
= k
, ++k
)
1211 if (k
!= pred
.begin())
1212 assert(pred
.key_comp()((*prev_k
).first
, (*k
).first
));
1213 for (i
= lt
.begin(); i
!= lt
.end(); prev
= i
, ++i
) {
1216 i_v
= (*i
).first
->eval(ind
->D
->sample
->p
);
1217 if (i
!= lt
.begin())
1218 assert(lt
.key_comp()((*prev
).first
, (*i
).first
));
1219 l
= eq
.find((*i
).first
);
1221 assert((*l
).second
.size() > 1);
1222 assert(head
.find((*i
).first
) != head
.end() ||
1223 pred
.find((*i
).first
) != pred
.end());
1224 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1225 k
= pred
.find((*i
).second
[j
]);
1226 assert(k
!= pred
.end());
1227 assert((*k
).second
!= 0);
1228 if ((*i
).first
->sign
!= 0 &&
1229 (*i
).second
[j
]->sign
!= 0 && ind
->D
->sample
) {
1230 vec_ZZ j_v
= (*i
).second
[j
]->eval(ind
->D
->sample
->p
);
1231 assert(lex_cmp(i_v
, j_v
) < 0);
1235 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1236 assert((*i
).second
.size() > 0);
1237 assert(head
.find((*i
).first
) != head
.end() ||
1238 pred
.find((*i
).first
) != pred
.end());
1239 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1240 k
= pred
.find((*i
).second
[j
]);
1241 assert(k
!= pred
.end());
1242 assert((*k
).second
!= 0);
1245 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1246 assert(head
.find((*i
).first
) != head
.end() ||
1247 pred
.find((*i
).first
) != pred
.end());
1248 assert((*i
).second
.size() >= 1);
1250 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1251 assert(head
.find((*i
).first
) != head
.end() ||
1252 pred
.find((*i
).first
) != pred
.end());
1253 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1254 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1255 pred
.find((*i
).second
[j
]) != pred
.end());
1259 max_term
* indicator::create_max_term(const indicator_term
*it
)
1261 int dim
= it
->den
.NumCols();
1262 int nparam
= ic
.PP
->Constraints
->NbColumns
-2 - ic
.vertex
.length();
1263 max_term
*maximum
= new max_term
;
1264 maximum
->domain
= new EDomain(D
);
1265 for (int j
= 0; j
< dim
; ++j
) {
1266 evalue
*E
= new evalue
;
1268 evalue_copy(E
, it
->vertex
[j
]);
1269 if (evalue_frac2floor_in_domain(E
, D
->D
))
1271 maximum
->max
.push_back(E
);
1276 static order_sign
evalue_sign(evalue
*diff
, EDomain
*D
, barvinok_options
*options
)
1278 order_sign sign
= order_eq
;
1281 evalue_set_si(&mone
, -1, 1);
1282 int len
= 1 + D
->D
->Dimension
+ 1;
1283 Vector
*c
= Vector_Alloc(len
);
1284 Matrix
*T
= Matrix_Alloc(2, len
-1);
1286 int fract
= evalue2constraint(D
, diff
, c
->p
, len
);
1287 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1288 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1290 order_sign upper_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1291 if (upper_sign
== order_lt
|| !fract
)
1295 evalue2constraint(D
, diff
, c
->p
, len
);
1297 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1298 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1300 order_sign neg_lower_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1302 if (neg_lower_sign
== order_lt
)
1304 else if (neg_lower_sign
== order_eq
|| neg_lower_sign
== order_le
) {
1305 if (upper_sign
== order_eq
|| upper_sign
== order_le
)
1310 if (upper_sign
== order_lt
|| upper_sign
== order_le
||
1311 upper_sign
== order_eq
)
1314 sign
= order_unknown
;
1320 free_evalue_refs(&mone
);
1325 /* An auxiliary class that keeps a reference to an evalue
1326 * and frees it when it goes out of scope.
1328 struct temp_evalue
{
1330 temp_evalue() : E(NULL
) {}
1331 temp_evalue(evalue
*e
) : E(e
) {}
1332 operator evalue
* () const { return E
; }
1333 evalue
*operator=(evalue
*e
) {
1335 free_evalue_refs(E
);
1343 free_evalue_refs(E
);
1349 static void substitute(vector
<indicator_term
*>& term
, evalue
*equation
)
1351 evalue
*fract
= NULL
;
1352 evalue
*val
= new evalue
;
1354 evalue_copy(val
, equation
);
1357 value_init(factor
.d
);
1358 value_init(factor
.x
.n
);
1361 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= fractional
;
1362 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1365 if (value_zero_p(e
->d
) && e
->x
.p
->type
== fractional
)
1366 fract
= &e
->x
.p
->arr
[0];
1368 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= polynomial
;
1369 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1371 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== polynomial
);
1374 int offset
= type_offset(e
->x
.p
);
1376 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].d
));
1377 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].x
.n
));
1378 if (value_neg_p(e
->x
.p
->arr
[offset
+1].x
.n
)) {
1379 value_oppose(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1380 value_assign(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1382 value_assign(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1383 value_oppose(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1386 free_evalue_refs(&e
->x
.p
->arr
[offset
+1]);
1389 *e
= e
->x
.p
->arr
[offset
];
1394 for (int i
= 0; i
< term
.size(); ++i
)
1395 term
[i
]->substitute(fract
, val
);
1397 free_evalue_refs(&p
->arr
[0]);
1399 for (int i
= 0; i
< term
.size(); ++i
)
1400 term
[i
]->substitute(p
->pos
, val
);
1403 free_evalue_refs(&factor
);
1404 free_evalue_refs(val
);
1410 order_sign
partial_order::compare(const indicator_term
*a
, const indicator_term
*b
)
1412 unsigned dim
= a
->den
.NumCols();
1413 order_sign sign
= order_eq
;
1414 EDomain
*D
= ind
->D
;
1415 unsigned MaxRays
= ind
->options
->verify
->barvinok
->MaxRays
;
1416 bool rational
= a
->sign
== 0 || b
->sign
== 0;
1418 order_sign cached_sign
= order_eq
;
1419 for (int k
= 0; k
< dim
; ++k
) {
1420 cached_sign
= evalue_diff_constant_sign(a
->vertex
[k
], b
->vertex
[k
]);
1421 if (cached_sign
!= order_eq
)
1424 if (cached_sign
!= order_undefined
)
1427 order_cache_el cache_el
;
1428 cached_sign
= order_undefined
;
1430 cached_sign
= cache
.check(a
, b
, cache_el
);
1431 if (cached_sign
!= order_undefined
) {
1436 if (rational
&& POL_ISSET(MaxRays
, POL_INTEGER
)) {
1437 ind
->options
->verify
->barvinok
->MaxRays
&= ~POL_INTEGER
;
1438 if (ind
->options
->verify
->barvinok
->MaxRays
)
1439 ind
->options
->verify
->barvinok
->MaxRays
|= POL_HIGH_BIT
;
1444 vector
<indicator_term
*> term
;
1446 for (int k
= 0; k
< dim
; ++k
) {
1447 /* compute a->vertex[k] - b->vertex[k] */
1449 if (cache_el
.e
.size() <= k
) {
1450 diff
= ediff(a
->vertex
[k
], b
->vertex
[k
]);
1451 cache_el
.e
.push_back(diff
);
1453 diff
= cache_el
.e
[k
];
1456 tdiff
= diff
= ediff(term
[0]->vertex
[k
], term
[1]->vertex
[k
]);
1457 order_sign diff_sign
;
1459 diff_sign
= order_undefined
;
1460 else if (eequal(a
->vertex
[k
], b
->vertex
[k
]))
1461 diff_sign
= order_eq
;
1463 diff_sign
= evalue_sign(diff
, D
, ind
->options
->verify
->barvinok
);
1465 if (diff_sign
== order_undefined
) {
1466 assert(sign
== order_le
|| sign
== order_ge
);
1467 if (sign
== order_le
)
1473 if (diff_sign
== order_lt
) {
1474 if (sign
== order_eq
|| sign
== order_le
)
1477 sign
= order_unknown
;
1480 if (diff_sign
== order_gt
) {
1481 if (sign
== order_eq
|| sign
== order_ge
)
1484 sign
= order_unknown
;
1487 if (diff_sign
== order_eq
) {
1488 if (D
== ind
->D
&& term
.size() == 0 && !rational
&&
1489 !EVALUE_IS_ZERO(*diff
))
1490 ind
->add_substitution(diff
);
1493 if ((diff_sign
== order_unknown
) ||
1494 ((diff_sign
== order_lt
|| diff_sign
== order_le
) && sign
== order_ge
) ||
1495 ((diff_sign
== order_gt
|| diff_sign
== order_ge
) && sign
== order_le
)) {
1496 sign
= order_unknown
;
1503 term
.push_back(new indicator_term(*a
));
1504 term
.push_back(new indicator_term(*b
));
1506 substitute(term
, diff
);
1510 cache
.add(cache_el
, sign
);
1514 if (D
&& D
!= ind
->D
)
1522 ind
->options
->verify
->barvinok
->MaxRays
= MaxRays
;
1526 bool partial_order::compared(const indicator_term
* a
, const indicator_term
* b
)
1528 order_type::iterator j
;
1531 if (j
!= lt
.end() && std::find(lt
[a
].begin(), lt
[a
].end(), b
) != lt
[a
].end())
1535 if (j
!= le
.end() && std::find(le
[a
].begin(), le
[a
].end(), b
) != le
[a
].end())
1541 void partial_order::add(const indicator_term
* it
,
1542 std::set
<const indicator_term
*> *filter
)
1544 if (eq
.find(it
) != eq
.end() && eq
[it
].size() == 1)
1547 head_type
head_copy(head
);
1552 head_type::iterator i
;
1553 for (i
= head_copy
.begin(); i
!= head_copy
.end(); ++i
) {
1556 if (eq
.find(*i
) != eq
.end() && eq
[*i
].size() == 1)
1559 if (filter
->find(*i
) == filter
->end())
1561 if (compared(*i
, it
))
1564 order_sign sign
= compare(it
, *i
);
1565 if (sign
== order_lt
) {
1566 lt
[it
].push_back(*i
);
1568 } else if (sign
== order_le
) {
1569 le
[it
].push_back(*i
);
1571 } else if (sign
== order_eq
) {
1574 } else if (sign
== order_gt
) {
1575 pending
[*i
].push_back(it
);
1576 lt
[*i
].push_back(it
);
1578 } else if (sign
== order_ge
) {
1579 pending
[*i
].push_back(it
);
1580 le
[*i
].push_back(it
);
1586 void partial_order::remove(const indicator_term
* it
)
1588 std::set
<const indicator_term
*> filter
;
1589 order_type::iterator i
;
1591 assert(head
.find(it
) != head
.end());
1594 if (i
!= eq
.end()) {
1595 assert(eq
[it
].size() >= 1);
1596 const indicator_term
*base
;
1597 if (eq
[it
].size() == 1) {
1601 vector
<const indicator_term
* >::iterator j
;
1602 j
= std::find(eq
[base
].begin(), eq
[base
].end(), it
);
1603 assert(j
!= eq
[base
].end());
1606 /* "it" may no longer be the smallest, since the order
1607 * structure may have been copied from another one.
1609 std::sort(eq
[it
].begin()+1, eq
[it
].end(), pred
.key_comp());
1610 assert(eq
[it
][0] == it
);
1611 eq
[it
].erase(eq
[it
].begin());
1616 for (int j
= 1; j
< eq
[base
].size(); ++j
)
1617 eq
[eq
[base
][j
]][0] = base
;
1620 if (i
!= lt
.end()) {
1626 if (i
!= le
.end()) {
1631 i
= pending
.find(it
);
1632 if (i
!= pending
.end()) {
1633 pending
[base
] = pending
[it
];
1638 if (eq
[base
].size() == 1)
1647 if (i
!= lt
.end()) {
1648 for (int j
= 0; j
< lt
[it
].size(); ++j
) {
1649 filter
.insert(lt
[it
][j
]);
1650 dec_pred(lt
[it
][j
]);
1656 if (i
!= le
.end()) {
1657 for (int j
= 0; j
< le
[it
].size(); ++j
) {
1658 filter
.insert(le
[it
][j
]);
1659 dec_pred(le
[it
][j
]);
1666 i
= pending
.find(it
);
1667 if (i
!= pending
.end()) {
1668 vector
<const indicator_term
*> it_pending
= pending
[it
];
1670 for (int j
= 0; j
< it_pending
.size(); ++j
) {
1671 filter
.erase(it_pending
[j
]);
1672 add(it_pending
[j
], &filter
);
1677 void partial_order::print(ostream
& os
, char **p
)
1679 order_type::iterator i
;
1680 pred_type::iterator j
;
1681 head_type::iterator k
;
1682 for (k
= head
.begin(); k
!= head
.end(); ++k
) {
1686 for (j
= pred
.begin(); j
!= pred
.end(); ++j
) {
1687 (*j
).first
->print(os
, p
);
1688 os
<< ": " << (*j
).second
<< endl
;
1690 for (i
= lt
.begin(); i
!= lt
.end(); ++i
) {
1691 (*i
).first
->print(os
, p
);
1692 assert(head
.find((*i
).first
) != head
.end() ||
1693 pred
.find((*i
).first
) != pred
.end());
1694 if (pred
.find((*i
).first
) != pred
.end())
1695 os
<< "(" << pred
[(*i
).first
] << ")";
1697 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1700 (*i
).second
[j
]->print(os
, p
);
1701 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1702 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1706 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1707 (*i
).first
->print(os
, p
);
1708 assert(head
.find((*i
).first
) != head
.end() ||
1709 pred
.find((*i
).first
) != pred
.end());
1710 if (pred
.find((*i
).first
) != pred
.end())
1711 os
<< "(" << pred
[(*i
).first
] << ")";
1713 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1716 (*i
).second
[j
]->print(os
, p
);
1717 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1718 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1722 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1723 if ((*i
).second
.size() <= 1)
1725 (*i
).first
->print(os
, p
);
1726 assert(head
.find((*i
).first
) != head
.end() ||
1727 pred
.find((*i
).first
) != pred
.end());
1728 if (pred
.find((*i
).first
) != pred
.end())
1729 os
<< "(" << pred
[(*i
).first
] << ")";
1730 for (int j
= 1; j
< (*i
).second
.size(); ++j
) {
1733 (*i
).second
[j
]->print(os
, p
);
1734 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1735 pred
.find((*i
).second
[j
]) != pred
.end());
1736 if (pred
.find((*i
).second
[j
]) != pred
.end())
1737 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1741 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1742 os
<< "pending on ";
1743 (*i
).first
->print(os
, p
);
1744 assert(head
.find((*i
).first
) != head
.end() ||
1745 pred
.find((*i
).first
) != pred
.end());
1746 if (pred
.find((*i
).first
) != pred
.end())
1747 os
<< "(" << pred
[(*i
).first
] << ")";
1749 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1752 (*i
).second
[j
]->print(os
, p
);
1753 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1754 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1760 void indicator::add(const indicator_term
* it
)
1762 indicator_term
*nt
= new indicator_term(*it
);
1763 if (options
->reduce
)
1764 nt
->reduce_in_domain(P
? P
: D
->D
);
1766 order
.add(nt
, NULL
);
1767 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1770 void indicator::remove(const indicator_term
* it
)
1772 vector
<indicator_term
*>::iterator i
;
1773 i
= std::find(term
.begin(), term
.end(), it
);
1774 assert(i
!= term
.end());
1777 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1781 void indicator::expand_rational_vertex(const indicator_term
*initial
)
1783 int pos
= initial
->pos
;
1785 if (ic
.terms
[pos
].size() == 0) {
1787 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, ic
.PP
) // _i is internal counter
1789 ic
.decompose_at_vertex(V
, pos
, options
->verify
->barvinok
);
1792 END_FORALL_PVertex_in_ParamPolyhedron
;
1794 for (int j
= 0; j
< ic
.terms
[pos
].size(); ++j
)
1795 add(ic
.terms
[pos
][j
]);
1798 void indicator::remove_initial_rational_vertices()
1801 const indicator_term
*initial
= NULL
;
1802 partial_order::head_type::iterator i
;
1803 for (i
= order
.head
.begin(); i
!= order
.head
.end(); ++i
) {
1804 if ((*i
)->sign
!= 0)
1806 if (order
.eq
.find(*i
) != order
.eq
.end() && order
.eq
[*i
].size() <= 1)
1813 expand_rational_vertex(initial
);
1817 void indicator::reduce_in_domain(Polyhedron
*D
)
1819 for (int i
= 0; i
< term
.size(); ++i
)
1820 term
[i
]->reduce_in_domain(D
);
1823 void indicator::print(ostream
& os
, char **p
)
1825 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1826 for (int i
= 0; i
< term
.size(); ++i
) {
1827 term
[i
]->print(os
, p
);
1829 os
<< ": " << term
[i
]->eval(D
->sample
->p
);
1836 /* Remove pairs of opposite terms */
1837 void indicator::simplify()
1839 for (int i
= 0; i
< term
.size(); ++i
) {
1840 for (int j
= i
+1; j
< term
.size(); ++j
) {
1841 if (term
[i
]->sign
+ term
[j
]->sign
!= 0)
1843 if (term
[i
]->den
!= term
[j
]->den
)
1846 for (k
= 0; k
< term
[i
]->den
.NumCols(); ++k
)
1847 if (!eequal(term
[i
]->vertex
[k
], term
[j
]->vertex
[k
]))
1849 if (k
< term
[i
]->den
.NumCols())
1853 term
.erase(term
.begin()+j
);
1854 term
.erase(term
.begin()+i
);
1861 void indicator::peel(int i
, int j
)
1869 int dim
= term
[i
]->den
.NumCols();
1874 int n_common
= 0, n_i
= 0, n_j
= 0;
1876 common
.SetDims(min(term
[i
]->den
.NumRows(), term
[j
]->den
.NumRows()), dim
);
1877 rest_i
.SetDims(term
[i
]->den
.NumRows(), dim
);
1878 rest_j
.SetDims(term
[j
]->den
.NumRows(), dim
);
1881 for (k
= 0, l
= 0; k
< term
[i
]->den
.NumRows() && l
< term
[j
]->den
.NumRows(); ) {
1882 int s
= lex_cmp(term
[i
]->den
[k
], term
[j
]->den
[l
]);
1884 common
[n_common
++] = term
[i
]->den
[k
];
1888 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1890 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1892 while (k
< term
[i
]->den
.NumRows())
1893 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1894 while (l
< term
[j
]->den
.NumRows())
1895 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1896 common
.SetDims(n_common
, dim
);
1897 rest_i
.SetDims(n_i
, dim
);
1898 rest_j
.SetDims(n_j
, dim
);
1900 for (k
= 0; k
<= n_i
; ++k
) {
1901 indicator_term
*it
= new indicator_term(*term
[i
]);
1902 it
->den
.SetDims(n_common
+ k
, dim
);
1903 for (l
= 0; l
< n_common
; ++l
)
1904 it
->den
[l
] = common
[l
];
1905 for (l
= 0; l
< k
; ++l
)
1906 it
->den
[n_common
+l
] = rest_i
[l
];
1907 lex_order_rows(it
->den
);
1909 for (l
= 0; l
< dim
; ++l
)
1910 evalue_add_constant(it
->vertex
[l
], rest_i
[k
-1][l
]);
1914 for (k
= 0; k
<= n_j
; ++k
) {
1915 indicator_term
*it
= new indicator_term(*term
[j
]);
1916 it
->den
.SetDims(n_common
+ k
, dim
);
1917 for (l
= 0; l
< n_common
; ++l
)
1918 it
->den
[l
] = common
[l
];
1919 for (l
= 0; l
< k
; ++l
)
1920 it
->den
[n_common
+l
] = rest_j
[l
];
1921 lex_order_rows(it
->den
);
1923 for (l
= 0; l
< dim
; ++l
)
1924 evalue_add_constant(it
->vertex
[l
], rest_j
[k
-1][l
]);
1927 term
.erase(term
.begin()+j
);
1928 term
.erase(term
.begin()+i
);
1931 void indicator::combine(const indicator_term
*a
, const indicator_term
*b
)
1933 int dim
= a
->den
.NumCols();
1936 mat_ZZ rest_i
; /* factors in a, but not in b */
1937 mat_ZZ rest_j
; /* factors in b, but not in a */
1938 int n_common
= 0, n_i
= 0, n_j
= 0;
1940 common
.SetDims(min(a
->den
.NumRows(), b
->den
.NumRows()), dim
);
1941 rest_i
.SetDims(a
->den
.NumRows(), dim
);
1942 rest_j
.SetDims(b
->den
.NumRows(), dim
);
1945 for (k
= 0, l
= 0; k
< a
->den
.NumRows() && l
< b
->den
.NumRows(); ) {
1946 int s
= lex_cmp(a
->den
[k
], b
->den
[l
]);
1948 common
[n_common
++] = a
->den
[k
];
1952 rest_i
[n_i
++] = a
->den
[k
++];
1954 rest_j
[n_j
++] = b
->den
[l
++];
1956 while (k
< a
->den
.NumRows())
1957 rest_i
[n_i
++] = a
->den
[k
++];
1958 while (l
< b
->den
.NumRows())
1959 rest_j
[n_j
++] = b
->den
[l
++];
1960 common
.SetDims(n_common
, dim
);
1961 rest_i
.SetDims(n_i
, dim
);
1962 rest_j
.SetDims(n_j
, dim
);
1964 assert(order
.eq
[a
].size() > 1);
1965 indicator_term
*prev
;
1968 for (int k
= n_i
-1; k
>= 0; --k
) {
1969 indicator_term
*it
= new indicator_term(*b
);
1970 it
->den
.SetDims(n_common
+ n_j
+ n_i
-k
, dim
);
1971 for (int l
= k
; l
< n_i
; ++l
)
1972 it
->den
[n_common
+n_j
+l
-k
] = rest_i
[l
];
1973 lex_order_rows(it
->den
);
1974 for (int m
= 0; m
< dim
; ++m
)
1975 evalue_add_constant(it
->vertex
[m
], rest_i
[k
][m
]);
1976 it
->sign
= -it
->sign
;
1978 order
.pending
[it
].push_back(prev
);
1979 order
.lt
[it
].push_back(prev
);
1980 order
.inc_pred(prev
);
1983 order
.head
.insert(it
);
1987 indicator_term
*it
= new indicator_term(*b
);
1988 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
1989 for (l
= 0; l
< n_i
; ++l
)
1990 it
->den
[n_common
+n_j
+l
] = rest_i
[l
];
1991 lex_order_rows(it
->den
);
1993 order
.pending
[a
].push_back(prev
);
1994 order
.lt
[a
].push_back(prev
);
1995 order
.inc_pred(prev
);
1996 order
.replace(b
, it
);
2001 for (int k
= n_j
-1; k
>= 0; --k
) {
2002 indicator_term
*it
= new indicator_term(*a
);
2003 it
->den
.SetDims(n_common
+ n_i
+ n_j
-k
, dim
);
2004 for (int l
= k
; l
< n_j
; ++l
)
2005 it
->den
[n_common
+n_i
+l
-k
] = rest_j
[l
];
2006 lex_order_rows(it
->den
);
2007 for (int m
= 0; m
< dim
; ++m
)
2008 evalue_add_constant(it
->vertex
[m
], rest_j
[k
][m
]);
2009 it
->sign
= -it
->sign
;
2011 order
.pending
[it
].push_back(prev
);
2012 order
.lt
[it
].push_back(prev
);
2013 order
.inc_pred(prev
);
2016 order
.head
.insert(it
);
2020 indicator_term
*it
= new indicator_term(*a
);
2021 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
2022 for (l
= 0; l
< n_j
; ++l
)
2023 it
->den
[n_common
+n_i
+l
] = rest_j
[l
];
2024 lex_order_rows(it
->den
);
2026 order
.pending
[a
].push_back(prev
);
2027 order
.lt
[a
].push_back(prev
);
2028 order
.inc_pred(prev
);
2029 order
.replace(a
, it
);
2033 assert(term
.size() == order
.head
.size() + order
.pred
.size());
2036 bool indicator::handle_equal_numerators(const indicator_term
*base
)
2038 for (int i
= 0; i
< order
.eq
[base
].size(); ++i
) {
2039 for (int j
= i
+1; j
< order
.eq
[base
].size(); ++j
) {
2040 if (order
.eq
[base
][i
]->is_opposite(order
.eq
[base
][j
])) {
2041 remove(order
.eq
[base
][j
]);
2042 remove(i
? order
.eq
[base
][i
] : base
);
2047 for (int j
= 1; j
< order
.eq
[base
].size(); ++j
)
2048 if (order
.eq
[base
][j
]->sign
!= base
->sign
) {
2049 combine(base
, order
.eq
[base
][j
]);
2055 void indicator::substitute(evalue
*equation
)
2057 ::substitute(term
, equation
);
2060 void indicator::add_substitution(evalue
*equation
)
2062 for (int i
= 0; i
< substitutions
.size(); ++i
)
2063 if (eequal(substitutions
[i
], equation
))
2065 evalue
*copy
= new evalue();
2066 value_init(copy
->d
);
2067 evalue_copy(copy
, equation
);
2068 substitutions
.push_back(copy
);
2071 void indicator::perform_pending_substitutions()
2073 if (substitutions
.size() == 0)
2076 for (int i
= 0; i
< substitutions
.size(); ++i
) {
2077 substitute(substitutions
[i
]);
2078 free_evalue_refs(substitutions
[i
]);
2079 delete substitutions
[i
];
2081 substitutions
.clear();
2085 static void print_varlist(ostream
& os
, int n
, char **names
)
2089 for (i
= 0; i
< n
; ++i
) {
2097 void max_term::print(ostream
& os
, char **p
, barvinok_options
*options
) const
2100 print_varlist(os
, domain
->dimension(), p
);
2103 for (int i
= 0; i
< max
.size(); ++i
) {
2106 evalue_print(os
, max
[i
], p
);
2110 domain
->print_constraints(os
, p
, options
);
2114 /* T maps the compressed parameters to the original parameters,
2115 * while this max_term is based on the compressed parameters
2116 * and we want get the original parameters back.
2118 void max_term::substitute(Matrix
*T
, barvinok_options
*options
)
2120 assert(domain
->dimension() == T
->NbColumns
-1);
2121 int nexist
= domain
->D
->Dimension
- (T
->NbColumns
-1);
2123 Matrix
*inv
= left_inverse(T
, &Eq
);
2126 value_init(denom
.d
);
2127 value_init(denom
.x
.n
);
2128 value_set_si(denom
.x
.n
, 1);
2129 value_assign(denom
.d
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
2132 v
.SetLength(inv
->NbColumns
);
2133 evalue
**subs
= new evalue
*[inv
->NbRows
-1];
2134 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2135 values2zz(inv
->p
[i
], v
, v
.length());
2136 subs
[i
] = multi_monom(v
);
2137 emul(&denom
, subs
[i
]);
2139 free_evalue_refs(&denom
);
2141 domain
->substitute(subs
, inv
, Eq
, options
->MaxRays
);
2144 for (int i
= 0; i
< max
.size(); ++i
) {
2145 evalue_substitute(max
[i
], subs
);
2146 reduce_evalue(max
[i
]);
2149 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2150 free_evalue_refs(subs
[i
]);
2157 Vector
*max_term::eval(Value
*val
, unsigned MaxRays
) const
2159 if (!domain
->contains(val
, domain
->dimension()))
2161 Vector
*res
= Vector_Alloc(max
.size());
2162 for (int i
= 0; i
< max
.size(); ++i
) {
2163 compute_evalue(max
[i
], val
, &res
->p
[i
]);
2170 enum sign
{ le
, ge
, lge
} sign
;
2172 split (evalue
*c
, enum sign s
) : constraint(c
), sign(s
) {}
2175 static void split_on(const split
& sp
, EDomain
*D
,
2176 EDomain
**Dlt
, EDomain
**Deq
, EDomain
**Dgt
,
2177 lexmin_options
*options
)
2183 ge_constraint
*ge
= D
->compute_ge_constraint(sp
.constraint
);
2184 if (sp
.sign
== split::lge
|| sp
.sign
== split::ge
)
2185 ED
[2] = EDomain::new_from_ge_constraint(ge
, 1, options
->verify
->barvinok
);
2188 if (sp
.sign
== split::lge
|| sp
.sign
== split::le
)
2189 ED
[0] = EDomain::new_from_ge_constraint(ge
, -1, options
->verify
->barvinok
);
2193 assert(sp
.sign
== split::lge
|| sp
.sign
== split::ge
|| sp
.sign
== split::le
);
2194 ED
[1] = EDomain::new_from_ge_constraint(ge
, 0, options
->verify
->barvinok
);
2198 for (int i
= 0; i
< 3; ++i
) {
2201 if (D
->sample
&& ED
[i
]->contains(D
->sample
->p
, D
->sample
->Size
-1)) {
2202 ED
[i
]->sample
= Vector_Alloc(D
->sample
->Size
);
2203 Vector_Copy(D
->sample
->p
, ED
[i
]->sample
->p
, D
->sample
->Size
);
2204 } else if (emptyQ2(ED
[i
]->D
) ||
2205 (options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2206 !(ED
[i
]->not_empty(options
)))) {
2216 ostream
& operator<< (ostream
& os
, const vector
<int> & v
)
2219 for (int i
= 0; i
< v
.size(); ++i
) {
2228 static bool isTranslation(Matrix
*M
)
2231 if (M
->NbRows
!= M
->NbColumns
)
2234 for (i
= 0;i
< M
->NbRows
; i
++)
2235 for (j
= 0; j
< M
->NbColumns
-1; j
++)
2237 if(value_notone_p(M
->p
[i
][j
]))
2240 if(value_notzero_p(M
->p
[i
][j
]))
2243 return value_one_p(M
->p
[M
->NbRows
-1][M
->NbColumns
-1]);
2246 static Matrix
*compress_parameters(Polyhedron
**P
, Polyhedron
**C
,
2247 unsigned nparam
, unsigned MaxRays
)
2251 /* compress_parms doesn't like equalities that only involve parameters */
2252 for (int i
= 0; i
< (*P
)->NbEq
; ++i
)
2253 assert(First_Non_Zero((*P
)->Constraint
[i
]+1, (*P
)->Dimension
-nparam
) != -1);
2255 M
= Matrix_Alloc((*P
)->NbEq
, (*P
)->Dimension
+2);
2256 Vector_Copy((*P
)->Constraint
[0], M
->p
[0], (*P
)->NbEq
* ((*P
)->Dimension
+2));
2257 CP
= compress_parms(M
, nparam
);
2260 if (isTranslation(CP
)) {
2265 T
= align_matrix(CP
, (*P
)->Dimension
+1);
2266 *P
= Polyhedron_Preimage(*P
, T
, MaxRays
);
2269 *C
= Polyhedron_Preimage(*C
, CP
, MaxRays
);
2274 void construct_rational_vertices(Param_Polyhedron
*PP
, Matrix
*T
, unsigned dim
,
2275 int nparam
, vector
<indicator_term
*>& vertices
)
2284 v
.SetLength(nparam
+1);
2287 value_init(factor
.d
);
2288 value_init(factor
.x
.n
);
2289 value_set_si(factor
.x
.n
, 1);
2290 value_set_si(factor
.d
, 1);
2292 for (i
= 0, PV
= PP
->V
; PV
; ++i
, PV
= PV
->next
) {
2293 indicator_term
*term
= new indicator_term(dim
, i
);
2294 vertices
.push_back(term
);
2295 Matrix
*M
= Matrix_Alloc(PV
->Vertex
->NbRows
+nparam
+1, nparam
+1);
2296 value_set_si(lcm
, 1);
2297 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
)
2298 value_lcm(lcm
, lcm
, PV
->Vertex
->p
[j
][nparam
+1]);
2299 value_assign(M
->p
[M
->NbRows
-1][M
->NbColumns
-1], lcm
);
2300 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
) {
2301 value_division(tmp
, lcm
, PV
->Vertex
->p
[j
][nparam
+1]);
2302 Vector_Scale(PV
->Vertex
->p
[j
], M
->p
[j
], tmp
, nparam
+1);
2304 for (int j
= 0; j
< nparam
; ++j
)
2305 value_assign(M
->p
[PV
->Vertex
->NbRows
+j
][j
], lcm
);
2307 Matrix
*M2
= Matrix_Alloc(T
->NbRows
, M
->NbColumns
);
2308 Matrix_Product(T
, M
, M2
);
2312 for (int j
= 0; j
< dim
; ++j
) {
2313 values2zz(M
->p
[j
], v
, nparam
+1);
2314 term
->vertex
[j
] = multi_monom(v
);
2315 value_assign(factor
.d
, lcm
);
2316 emul(&factor
, term
->vertex
[j
]);
2320 assert(i
== PP
->nbV
);
2321 free_evalue_refs(&factor
);
2326 static vector
<max_term
*> lexmin(indicator
& ind
, unsigned nparam
,
2329 vector
<max_term
*> maxima
;
2330 partial_order::head_type::iterator i
;
2331 vector
<int> best_score
;
2332 vector
<int> second_score
;
2333 vector
<int> neg_score
;
2336 ind
.perform_pending_substitutions();
2337 const indicator_term
*best
= NULL
, *second
= NULL
, *neg
= NULL
,
2338 *neg_eq
= NULL
, *neg_le
= NULL
;
2339 for (i
= ind
.order
.head
.begin(); i
!= ind
.order
.head
.end(); ++i
) {
2341 const indicator_term
*term
= *i
;
2342 if (term
->sign
== 0) {
2343 ind
.expand_rational_vertex(term
);
2347 if (ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2349 if (ind
.order
.eq
[term
].size() <= 1)
2351 for (j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2352 if (ind
.order
.pred
.find(ind
.order
.eq
[term
][j
]) !=
2353 ind
.order
.pred
.end())
2355 if (j
< ind
.order
.eq
[term
].size())
2357 score
.push_back(ind
.order
.eq
[term
].size());
2360 if (ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2361 score
.push_back(ind
.order
.le
[term
].size());
2364 if (ind
.order
.lt
.find(term
) != ind
.order
.lt
.end())
2365 score
.push_back(-ind
.order
.lt
[term
].size());
2369 if (term
->sign
> 0) {
2370 if (!best
|| score
< best_score
) {
2372 second_score
= best_score
;
2375 } else if (!second
|| score
< second_score
) {
2377 second_score
= score
;
2380 if (!neg_eq
&& ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2381 for (int j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2382 if (ind
.order
.eq
[term
][j
]->sign
!= term
->sign
) {
2387 if (!neg_le
&& ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2389 if (!neg
|| score
< neg_score
) {
2395 if (i
!= ind
.order
.head
.end())
2398 if (!best
&& neg_eq
) {
2399 assert(ind
.order
.eq
[neg_eq
].size() != 0);
2400 bool handled
= ind
.handle_equal_numerators(neg_eq
);
2405 if (!best
&& neg_le
) {
2406 /* The smallest term is negative and <= some positive term */
2412 /* apparently there can be negative initial term on empty domains */
2413 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2414 ind
.options
->verify
->barvinok
->lp_solver
== BV_LP_POLYLIB
)
2419 if (!second
&& !neg
) {
2420 const indicator_term
*rat
= NULL
;
2422 if (ind
.order
.le
.find(best
) == ind
.order
.le
.end()) {
2423 if (ind
.order
.eq
.find(best
) != ind
.order
.eq
.end()) {
2424 bool handled
= ind
.handle_equal_numerators(best
);
2425 if (ind
.options
->emptiness_check
!=
2426 BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2427 ind
.options
->verify
->barvinok
->lp_solver
== BV_LP_POLYLIB
)
2429 /* If !handled then the leading coefficient is bigger than one;
2430 * must be an empty domain
2437 maxima
.push_back(ind
.create_max_term(best
));
2440 for (int j
= 0; j
< ind
.order
.le
[best
].size(); ++j
) {
2441 if (ind
.order
.le
[best
][j
]->sign
== 0) {
2442 if (!rat
&& ind
.order
.pred
[ind
.order
.le
[best
][j
]] == 1)
2443 rat
= ind
.order
.le
[best
][j
];
2444 } else if (ind
.order
.le
[best
][j
]->sign
> 0) {
2445 second
= ind
.order
.le
[best
][j
];
2448 neg
= ind
.order
.le
[best
][j
];
2451 if (!second
&& !neg
) {
2453 ind
.order
.unset_le(best
, rat
);
2454 ind
.expand_rational_vertex(rat
);
2461 ind
.order
.unset_le(best
, second
);
2467 unsigned dim
= best
->den
.NumCols();
2470 for (int k
= 0; k
< dim
; ++k
) {
2471 diff
= ediff(best
->vertex
[k
], second
->vertex
[k
]);
2472 sign
= evalue_sign(diff
, ind
.D
, ind
.options
->verify
->barvinok
);
2474 /* neg can never be smaller than best, unless it may still cancel.
2475 * This can happen if positive terms have been determined to
2476 * be equal or less than or equal to some negative term.
2478 if (second
== neg
&& !neg_eq
&& !neg_le
) {
2479 if (sign
== order_ge
)
2481 if (sign
== order_unknown
)
2485 if (sign
!= order_eq
)
2487 if (!EVALUE_IS_ZERO(*diff
)) {
2488 ind
.add_substitution(diff
);
2489 ind
.perform_pending_substitutions();
2492 if (sign
== order_eq
) {
2493 ind
.order
.set_equal(best
, second
);
2496 if (sign
== order_lt
) {
2497 ind
.order
.lt
[best
].push_back(second
);
2498 ind
.order
.inc_pred(second
);
2501 if (sign
== order_gt
) {
2502 ind
.order
.lt
[second
].push_back(best
);
2503 ind
.order
.inc_pred(best
);
2507 split
sp(diff
, sign
== order_le
? split::le
:
2508 sign
== order_ge
? split::ge
: split::lge
);
2510 EDomain
*Dlt
, *Deq
, *Dgt
;
2511 split_on(sp
, ind
.D
, &Dlt
, &Deq
, &Dgt
, ind
.options
);
2512 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
)
2513 assert(Dlt
|| Deq
|| Dgt
);
2514 else if (!(Dlt
|| Deq
|| Dgt
))
2515 /* Must have been empty all along */
2518 if (Deq
&& (Dlt
|| Dgt
)) {
2519 int locsize
= loc
.size();
2521 indicator
indeq(ind
, Deq
);
2523 indeq
.add_substitution(diff
);
2524 indeq
.perform_pending_substitutions();
2525 vector
<max_term
*> maxeq
= lexmin(indeq
, nparam
, loc
);
2526 maxima
.insert(maxima
.end(), maxeq
.begin(), maxeq
.end());
2527 loc
.resize(locsize
);
2530 int locsize
= loc
.size();
2532 indicator
indgt(ind
, Dgt
);
2534 /* we don't know the new location of these terms in indgt */
2536 indgt.order.lt[second].push_back(best);
2537 indgt.order.inc_pred(best);
2539 vector
<max_term
*> maxgt
= lexmin(indgt
, nparam
, loc
);
2540 maxima
.insert(maxima
.end(), maxgt
.begin(), maxgt
.end());
2541 loc
.resize(locsize
);
2546 ind
.set_domain(Deq
);
2547 ind
.add_substitution(diff
);
2548 ind
.perform_pending_substitutions();
2552 ind
.set_domain(Dlt
);
2553 ind
.order
.lt
[best
].push_back(second
);
2554 ind
.order
.inc_pred(second
);
2558 ind
.set_domain(Dgt
);
2559 ind
.order
.lt
[second
].push_back(best
);
2560 ind
.order
.inc_pred(best
);
2567 static void lexmin_base(Polyhedron
*P
, Polyhedron
*C
,
2568 Matrix
*CP
, Matrix
*T
,
2569 vector
<max_term
*>& all_max
,
2570 lexmin_options
*options
)
2572 unsigned nparam
= C
->Dimension
;
2573 Param_Polyhedron
*PP
= NULL
;
2575 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
->verify
->barvinok
);
2577 unsigned dim
= P
->Dimension
- nparam
;
2581 indicator_constructor
ic(P
, dim
, PP
, T
);
2583 vector
<indicator_term
*> all_vertices
;
2584 construct_rational_vertices(PP
, T
, T
? T
->NbRows
-nparam
-1 : dim
,
2585 nparam
, all_vertices
);
2587 Polyhedron
*TC
= true_context(P
, C
, options
->verify
->barvinok
->MaxRays
);
2588 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
->verify
->barvinok
, i
, D
, rVD
)
2591 EDomain
*epVD
= new EDomain(rVD
);
2592 indicator
ind(ic
, D
, epVD
, options
);
2594 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
2595 ind
.add(all_vertices
[_i
]);
2596 END_FORALL_PVertex_in_ParamPolyhedron
;
2598 ind
.remove_initial_rational_vertices();
2601 vector
<max_term
*> maxima
= lexmin(ind
, nparam
, loc
);
2603 for (int j
= 0; j
< maxima
.size(); ++j
)
2604 maxima
[j
]->substitute(CP
, options
->verify
->barvinok
);
2605 all_max
.insert(all_max
.end(), maxima
.begin(), maxima
.end());
2608 END_FORALL_REDUCED_DOMAIN
2609 Polyhedron_Free(TC
);
2610 for (int i
= 0; i
< all_vertices
.size(); ++i
)
2611 delete all_vertices
[i
];
2612 Param_Polyhedron_Free(PP
);
2615 static vector
<max_term
*> lexmin(Polyhedron
*P
, Polyhedron
*C
,
2616 lexmin_options
*options
)
2618 unsigned nparam
= C
->Dimension
;
2619 Matrix
*T
= NULL
, *CP
= NULL
;
2620 Polyhedron
*Porig
= P
;
2621 Polyhedron
*Corig
= C
;
2622 vector
<max_term
*> all_max
;
2627 POL_ENSURE_VERTICES(P
);
2632 assert(P
->NbBid
== 0);
2635 remove_all_equalities(&P
, &C
, &CP
, &T
, nparam
,
2636 options
->verify
->barvinok
->MaxRays
);
2638 lexmin_base(P
, C
, CP
, T
, all_max
, options
);
2651 static void verify_results(Polyhedron
*A
, Polyhedron
*C
,
2652 vector
<max_term
*>& maxima
,
2653 struct verify_options
*options
);
2655 /* Turn the set dimensions of "context" into parameters and return
2656 * the corresponding parameter domain.
2658 static struct isl_basic_set
*to_parameter_domain(struct isl_basic_set
*context
)
2660 context
= isl_basic_set_move_dims(context
, isl_dim_param
, 0,
2661 isl_dim_set
, 0, isl_basic_set_dim(context
, isl_dim_set
));
2662 context
= isl_basic_set_params(context
);
2666 int main(int argc
, char **argv
)
2669 isl_basic_set
*context
, *bset
;
2674 int urs_unknowns
= 0;
2675 int print_solution
= 1;
2676 struct lexmin_options
*options
= lexmin_options_new_with_defaults();
2678 options
->verify
->barvinok
->lookup_table
= 0;
2680 argc
= lexmin_options_parse(options
, argc
, argv
, ISL_ARG_ALL
);
2681 ctx
= isl_ctx_alloc_with_options(&lexmin_options_args
, options
);
2683 context
= isl_basic_set_read_from_file(ctx
, stdin
);
2685 n
= fscanf(stdin
, "%d", &neg_one
);
2687 assert(neg_one
== -1);
2688 bset
= isl_basic_set_read_from_file(ctx
, stdin
);
2690 while (fgets(s
, sizeof(s
), stdin
)) {
2691 if (strncasecmp(s
, "Maximize", 8) == 0) {
2692 fprintf(stderr
, "Maximize option not supported\n");
2695 if (strncasecmp(s
, "Rational", 8) == 0) {
2696 fprintf(stderr
, "Rational option not supported\n");
2699 if (strncasecmp(s
, "Urs_parms", 9) == 0)
2701 if (strncasecmp(s
, "Urs_unknowns", 12) == 0)
2705 context
= isl_basic_set_intersect(context
,
2706 isl_basic_set_positive_orthant(isl_basic_set_get_space(context
)));
2707 context
= to_parameter_domain(context
);
2708 nparam
= isl_basic_set_dim(context
, isl_dim_param
);
2709 if (nparam
!= isl_basic_set_dim(bset
, isl_dim_param
)) {
2710 int dim
= isl_basic_set_dim(bset
, isl_dim_set
);
2711 bset
= isl_basic_set_move_dims(bset
, isl_dim_param
, 0,
2712 isl_dim_set
, dim
- nparam
, nparam
);
2715 bset
= isl_basic_set_intersect(bset
,
2716 isl_basic_set_positive_orthant(isl_basic_set_get_space(bset
)));
2718 if (options
->verify
->verify
)
2721 A
= isl_basic_set_to_polylib(bset
);
2722 verify_options_set_range(options
->verify
, A
->Dimension
);
2723 C
= isl_basic_set_to_polylib(context
);
2724 vector
<max_term
*> maxima
= lexmin(A
, C
, options
);
2725 if (print_solution
) {
2727 param_names
= util_generate_names(C
->Dimension
, "p");
2728 for (int i
= 0; i
< maxima
.size(); ++i
)
2729 maxima
[i
]->print(cout
, param_names
,
2730 options
->verify
->barvinok
);
2731 util_free_names(C
->Dimension
, param_names
);
2734 if (options
->verify
->verify
)
2735 verify_results(A
, C
, maxima
, options
->verify
);
2737 for (int i
= 0; i
< maxima
.size(); ++i
)
2743 isl_basic_set_free(bset
);
2744 isl_basic_set_free(context
);
2750 static bool lexmin(int pos
, Polyhedron
*P
, Value
*context
)
2759 value_init(LB
); value_init(UB
); value_init(k
);
2762 lu_flags
= lower_upper_bounds(pos
,P
,context
,&LB
,&UB
);
2763 assert(!(lu_flags
& LB_INFINITY
));
2765 value_set_si(context
[pos
],0);
2766 if (!lu_flags
&& value_lt(UB
,LB
)) {
2767 value_clear(LB
); value_clear(UB
); value_clear(k
);
2771 value_assign(context
[pos
], LB
);
2772 value_clear(LB
); value_clear(UB
); value_clear(k
);
2775 for (value_assign(k
,LB
); lu_flags
|| value_le(k
,UB
); value_increment(k
,k
)) {
2776 value_assign(context
[pos
],k
);
2777 if ((found
= lexmin(pos
+1, P
->next
, context
)))
2781 value_set_si(context
[pos
],0);
2782 value_clear(LB
); value_clear(UB
); value_clear(k
);
2786 static void print_list(FILE *out
, Value
*z
, const char* brackets
, int len
)
2788 fprintf(out
, "%c", brackets
[0]);
2789 value_print(out
, VALUE_FMT
,z
[0]);
2790 for (int k
= 1; k
< len
; ++k
) {
2792 value_print(out
, VALUE_FMT
,z
[k
]);
2794 fprintf(out
, "%c", brackets
[1]);
2797 static int check_poly_lexmin(const struct check_poly_data
*data
,
2798 int nparam
, Value
*z
,
2799 const struct verify_options
*options
);
2801 struct check_poly_lexmin_data
: public check_poly_data
{
2803 vector
<max_term
*>& maxima
;
2805 check_poly_lexmin_data(Polyhedron
*S
, Value
*z
,
2806 vector
<max_term
*>& maxima
) : S(S
), maxima(maxima
) {
2808 this->check
= &check_poly_lexmin
;
2812 static int check_poly_lexmin(const struct check_poly_data
*data
,
2813 int nparam
, Value
*z
,
2814 const struct verify_options
*options
)
2816 const check_poly_lexmin_data
*lexmin_data
;
2817 lexmin_data
= static_cast<const check_poly_lexmin_data
*>(data
);
2818 Polyhedron
*S
= lexmin_data
->S
;
2819 vector
<max_term
*>& maxima
= lexmin_data
->maxima
;
2821 bool found
= lexmin(1, S
, lexmin_data
->z
);
2823 if (options
->print_all
) {
2825 print_list(stdout
, z
, "()", nparam
);
2828 print_list(stdout
, lexmin_data
->z
+1, "[]", S
->Dimension
-nparam
);
2833 for (int i
= 0; i
< maxima
.size(); ++i
)
2834 if ((min
= maxima
[i
]->eval(z
, options
->barvinok
->MaxRays
)))
2837 int ok
= !(found
^ !!min
);
2839 for (int i
= 0; i
< S
->Dimension
-nparam
; ++i
)
2840 if (value_ne(lexmin_data
->z
[1+i
], min
->p
[i
])) {
2847 fprintf(stderr
, "Error !\n");
2848 fprintf(stderr
, "lexmin");
2849 print_list(stderr
, z
, "()", nparam
);
2850 fprintf(stderr
, " should be ");
2852 print_list(stderr
, lexmin_data
->z
+1, "[]", S
->Dimension
-nparam
);
2853 fprintf(stderr
, " while digging gives ");
2855 print_list(stderr
, min
->p
, "[]", S
->Dimension
-nparam
);
2856 fprintf(stderr
, ".\n");
2858 } else if (options
->print_all
)
2863 for (k
= 1; k
<= S
->Dimension
-nparam
; ++k
)
2864 value_set_si(lexmin_data
->z
[k
], 0);
2869 void verify_results(Polyhedron
*A
, Polyhedron
*C
, vector
<max_term
*>& maxima
,
2870 struct verify_options
*options
)
2873 unsigned nparam
= C
->Dimension
;
2874 unsigned MaxRays
= options
->barvinok
->MaxRays
;
2879 CS
= check_poly_context_scan(A
, &C
, nparam
, options
);
2881 p
= Vector_Alloc(A
->Dimension
+2);
2882 value_set_si(p
->p
[A
->Dimension
+1], 1);
2884 S
= Polyhedron_Scan(A
, C
, MaxRays
& POL_NO_DUAL
? 0 : MaxRays
);
2886 check_poly_init(C
, options
);
2889 if (!(CS
&& emptyQ2(CS
))) {
2890 check_poly_lexmin_data
data(S
, p
->p
, maxima
);
2891 check_poly(CS
, &data
, nparam
, 0, p
->p
+S
->Dimension
-nparam
+1, options
);
2896 if (!options
->print_all
)