6 #include <NTL/vec_ZZ.h>
7 #include <NTL/mat_ZZ.h>
8 #include <barvinok/barvinok.h>
9 #include <barvinok/evalue.h>
10 #include <barvinok/options.h>
11 #include <barvinok/util.h>
13 #include "conversion.h"
14 #include "decomposer.h"
15 #include "lattice_point.h"
16 #include "reduce_domain.h"
20 #include "evalue_util.h"
21 #include "remove_equalities.h"
37 #define EMPTINESS_CHECK (BV_OPT_LAST+1)
38 #define NO_REDUCTION (BV_OPT_LAST+2)
39 #define POLYSIGN (BV_OPT_LAST+3)
41 struct argp_option argp_options
[] = {
42 { "emptiness-check", EMPTINESS_CHECK
, "[none|count]", 0 },
43 { "no-reduction", NO_REDUCTION
, 0, 0 },
44 { "polysign", POLYSIGN
, "[cdd|cddf]", 0 },
48 static error_t
parse_opt(int key
, char *arg
, struct argp_state
*state
)
50 struct lexmin_options
*options
= (struct lexmin_options
*)(state
->input
);
51 struct barvinok_options
*bv_options
= options
->barvinok
;
55 state
->child_inputs
[0] = options
->barvinok
;
56 state
->child_inputs
[1] = &options
->verify
;
57 options
->emptiness_check
= BV_LEXMIN_EMPTINESS_CHECK_SAMPLE
;
59 options
->polysign
= BV_LEXMIN_POLYSIGN_POLYLIB
;
62 if (!strcmp(arg
, "none"))
63 options
->emptiness_check
= BV_LEXMIN_EMPTINESS_CHECK_NONE
;
64 else if (!strcmp(arg
, "count")) {
65 options
->emptiness_check
= BV_LEXMIN_EMPTINESS_CHECK_COUNT
;
66 bv_options
->count_sample_infinite
= 0;
73 if (!strcmp(arg
, "cddf"))
74 options
->polysign
= BV_LEXMIN_POLYSIGN_CDDF
;
75 else if (!strcmp(arg
, "cdd"))
76 options
->polysign
= BV_LEXMIN_POLYSIGN_CDD
;
79 return ARGP_ERR_UNKNOWN
;
84 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
86 static int type_offset(enode
*p
)
88 return p
->type
== fractional
? 1 :
89 p
->type
== flooring
? 1 : 0;
92 void compute_evalue(evalue
*e
, Value
*val
, Value
*res
)
94 double d
= compute_evalue(e
, val
);
99 value_set_double(*res
, d
);
102 struct indicator_term
{
104 int pos
; /* number of rational vertex */
105 int n
; /* number of cone associated to given rational vertex */
109 indicator_term(unsigned dim
, int pos
) {
111 vertex
= new evalue
* [dim
];
116 indicator_term(unsigned dim
, int pos
, int n
) {
117 den
.SetDims(dim
, dim
);
118 vertex
= new evalue
* [dim
];
122 indicator_term(const indicator_term
& src
) {
127 unsigned dim
= den
.NumCols();
128 vertex
= new evalue
* [dim
];
129 for (int i
= 0; i
< dim
; ++i
) {
130 vertex
[i
] = new evalue();
131 value_init(vertex
[i
]->d
);
132 evalue_copy(vertex
[i
], src
.vertex
[i
]);
135 void swap(indicator_term
*other
) {
137 tmp
= sign
; sign
= other
->sign
; other
->sign
= tmp
;
138 tmp
= pos
; pos
= other
->pos
; other
->pos
= tmp
;
139 tmp
= n
; n
= other
->n
; other
->n
= tmp
;
140 mat_ZZ tmp_den
= den
; den
= other
->den
; other
->den
= tmp_den
;
141 unsigned dim
= den
.NumCols();
142 for (int i
= 0; i
< dim
; ++i
) {
143 evalue
*tmp
= vertex
[i
];
144 vertex
[i
] = other
->vertex
[i
];
145 other
->vertex
[i
] = tmp
;
149 unsigned dim
= den
.NumCols();
150 for (int i
= 0; i
< dim
; ++i
) {
151 free_evalue_refs(vertex
[i
]);
156 void print(ostream
& os
, char **p
) const;
157 void substitute(Matrix
*T
);
159 void substitute(evalue
*fract
, evalue
*val
);
160 void substitute(int pos
, evalue
*val
);
161 void reduce_in_domain(Polyhedron
*D
);
162 bool is_opposite(const indicator_term
*neg
) const;
163 vec_ZZ
eval(Value
*val
) const {
165 unsigned dim
= den
.NumCols();
169 for (int i
= 0; i
< dim
; ++i
) {
170 compute_evalue(vertex
[i
], val
, &tmp
);
178 static int evalue_rational_cmp(const evalue
*e1
, const evalue
*e2
)
186 assert(value_notzero_p(e1
->d
));
187 assert(value_notzero_p(e2
->d
));
188 value_multiply(m
, e1
->x
.n
, e2
->d
);
189 value_multiply(m2
, e2
->x
.n
, e1
->d
);
192 else if (value_gt(m
, m2
))
202 static int evalue_cmp(const evalue
*e1
, const evalue
*e2
)
204 if (value_notzero_p(e1
->d
)) {
205 if (value_zero_p(e2
->d
))
207 return evalue_rational_cmp(e1
, e2
);
209 if (value_notzero_p(e2
->d
))
211 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
212 return e1
->x
.p
->type
- e2
->x
.p
->type
;
213 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
214 return e1
->x
.p
->size
- e2
->x
.p
->size
;
215 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
216 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
217 assert(e1
->x
.p
->type
== polynomial
||
218 e1
->x
.p
->type
== fractional
||
219 e1
->x
.p
->type
== flooring
);
220 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
) {
221 int s
= evalue_cmp(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]);
228 void evalue_length(evalue
*e
, int len
[2])
233 while (value_zero_p(e
->d
)) {
234 assert(e
->x
.p
->type
== polynomial
||
235 e
->x
.p
->type
== fractional
||
236 e
->x
.p
->type
== flooring
);
237 if (e
->x
.p
->type
== polynomial
)
241 int offset
= type_offset(e
->x
.p
);
242 assert(e
->x
.p
->size
== offset
+2);
243 e
= &e
->x
.p
->arr
[offset
];
247 static bool it_smaller(const indicator_term
* it1
, const indicator_term
* it2
)
251 int len1
[2], len2
[2];
252 unsigned dim
= it1
->den
.NumCols();
253 for (int i
= 0; i
< dim
; ++i
) {
254 evalue_length(it1
->vertex
[i
], len1
);
255 evalue_length(it2
->vertex
[i
], len2
);
256 if (len1
[0] != len2
[0])
257 return len1
[0] < len2
[0];
258 if (len1
[1] != len2
[1])
259 return len1
[1] < len2
[1];
261 if (it1
->pos
!= it2
->pos
)
262 return it1
->pos
< it2
->pos
;
263 if (it1
->n
!= it2
->n
)
264 return it1
->n
< it2
->n
;
265 int s
= lex_cmp(it1
->den
, it2
->den
);
268 for (int i
= 0; i
< dim
; ++i
) {
269 s
= evalue_cmp(it1
->vertex
[i
], it2
->vertex
[i
]);
273 assert(it1
->sign
!= 0);
274 assert(it2
->sign
!= 0);
275 if (it1
->sign
!= it2
->sign
)
276 return it1
->sign
> 0;
281 static const int requires_resort
;
282 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
283 return it_smaller(it1
, it2
);
286 const int smaller_it::requires_resort
= 1;
288 struct smaller_it_p
{
289 static const int requires_resort
;
290 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
294 const int smaller_it_p::requires_resort
= 0;
296 /* Returns true if this and neg are opposite using the knowledge
297 * that they have the same numerator.
298 * In particular, we check that the signs are different and that
299 * the denominator is the same.
301 bool indicator_term::is_opposite(const indicator_term
*neg
) const
303 if (sign
+ neg
->sign
!= 0)
310 void indicator_term::reduce_in_domain(Polyhedron
*D
)
312 for (int k
= 0; k
< den
.NumCols(); ++k
) {
313 reduce_evalue_in_domain(vertex
[k
], D
);
314 if (evalue_range_reduction_in_domain(vertex
[k
], D
))
315 reduce_evalue(vertex
[k
]);
319 void indicator_term::print(ostream
& os
, char **p
) const
321 unsigned dim
= den
.NumCols();
322 unsigned factors
= den
.NumRows();
330 for (int i
= 0; i
< dim
; ++i
) {
333 evalue_print(os
, vertex
[i
], p
);
336 for (int i
= 0; i
< factors
; ++i
) {
337 os
<< " + t" << i
<< "*[";
338 for (int j
= 0; j
< dim
; ++j
) {
345 os
<< " ((" << pos
<< ", " << n
<< ", " << (void*)this << "))";
348 /* Perform the substitution specified by T on the variables.
349 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
350 * The substitution is performed as in gen_fun::substitute
352 void indicator_term::substitute(Matrix
*T
)
354 unsigned dim
= den
.NumCols();
355 unsigned nparam
= T
->NbColumns
- dim
- 1;
356 unsigned newdim
= T
->NbRows
- nparam
- 1;
359 matrix2zz(T
, trans
, newdim
, dim
);
360 trans
= transpose(trans
);
362 newvertex
= new evalue
* [newdim
];
365 v
.SetLength(nparam
+1);
368 value_init(factor
.d
);
369 value_set_si(factor
.d
, 1);
370 value_init(factor
.x
.n
);
371 for (int i
= 0; i
< newdim
; ++i
) {
372 values2zz(T
->p
[i
]+dim
, v
, nparam
+1);
373 newvertex
[i
] = multi_monom(v
);
375 for (int j
= 0; j
< dim
; ++j
) {
376 if (value_zero_p(T
->p
[i
][j
]))
380 evalue_copy(&term
, vertex
[j
]);
381 value_assign(factor
.x
.n
, T
->p
[i
][j
]);
382 emul(&factor
, &term
);
383 eadd(&term
, newvertex
[i
]);
384 free_evalue_refs(&term
);
387 free_evalue_refs(&factor
);
388 for (int i
= 0; i
< dim
; ++i
) {
389 free_evalue_refs(vertex
[i
]);
396 static void evalue_add_constant(evalue
*e
, ZZ v
)
401 /* go down to constant term */
402 while (value_zero_p(e
->d
))
403 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
406 value_multiply(tmp
, tmp
, e
->d
);
407 value_addto(e
->x
.n
, e
->x
.n
, tmp
);
412 /* Make all powers in denominator lexico-positive */
413 void indicator_term::normalize()
416 extra_vertex
.SetLength(den
.NumCols());
417 for (int r
= 0; r
< den
.NumRows(); ++r
) {
418 for (int k
= 0; k
< den
.NumCols(); ++k
) {
425 extra_vertex
+= den
[r
];
429 for (int k
= 0; k
< extra_vertex
.length(); ++k
)
430 if (extra_vertex
[k
] != 0)
431 evalue_add_constant(vertex
[k
], extra_vertex
[k
]);
434 static void substitute(evalue
*e
, evalue
*fract
, evalue
*val
)
438 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
439 if (t
->x
.p
->type
== fractional
&& eequal(&t
->x
.p
->arr
[0], fract
))
442 if (value_notzero_p(t
->d
))
445 free_evalue_refs(&t
->x
.p
->arr
[0]);
446 evalue
*term
= &t
->x
.p
->arr
[2];
453 free_evalue_refs(term
);
459 void indicator_term::substitute(evalue
*fract
, evalue
*val
)
461 unsigned dim
= den
.NumCols();
462 for (int i
= 0; i
< dim
; ++i
) {
463 ::substitute(vertex
[i
], fract
, val
);
467 static void substitute(evalue
*e
, int pos
, evalue
*val
)
471 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
472 if (t
->x
.p
->type
== polynomial
&& t
->x
.p
->pos
== pos
)
475 if (value_notzero_p(t
->d
))
478 evalue
*term
= &t
->x
.p
->arr
[1];
485 free_evalue_refs(term
);
491 void indicator_term::substitute(int pos
, evalue
*val
)
493 unsigned dim
= den
.NumCols();
494 for (int i
= 0; i
< dim
; ++i
) {
495 ::substitute(vertex
[i
], pos
, val
);
499 struct indicator_constructor
: public signed_cone_consumer
,
500 public vertex_decomposer
{
502 vector
<indicator_term
*> *terms
;
503 Matrix
*T
; /* Transformation to original space */
504 Param_Polyhedron
*PP
;
508 indicator_constructor(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
,
510 vertex_decomposer(P
, PP
->nbV
, *this), T(T
), PP(PP
) {
511 vertex
.SetLength(dim
);
512 terms
= new vector
<indicator_term
*>[nbV
];
514 ~indicator_constructor() {
515 for (int i
= 0; i
< nbV
; ++i
)
516 for (int j
= 0; j
< terms
[i
].size(); ++j
)
520 void substitute(Matrix
*T
);
522 void print(ostream
& os
, char **p
);
524 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
525 void decompose_at_vertex(Param_Vertices
*V
, int _i
,
526 barvinok_options
*options
) {
529 vertex_decomposer::decompose_at_vertex(V
, _i
, options
);
533 void indicator_constructor::handle(const signed_cone
& sc
, barvinok_options
*options
)
536 unsigned dim
= vertex
.length();
538 assert(sc
.rays
.NumRows() == dim
);
540 indicator_term
*term
= new indicator_term(dim
, pos
, n
++);
541 term
->sign
= sc
.sign
;
542 terms
[vert
].push_back(term
);
544 lattice_point(V
, sc
.rays
, vertex
, term
->vertex
, options
);
547 for (int r
= 0; r
< dim
; ++r
) {
548 for (int j
= 0; j
< dim
; ++j
) {
549 if (term
->den
[r
][j
] == 0)
551 if (term
->den
[r
][j
] > 0)
553 term
->sign
= -term
->sign
;
554 term
->den
[r
] = -term
->den
[r
];
555 vertex
+= term
->den
[r
];
560 for (int i
= 0; i
< dim
; ++i
) {
561 if (!term
->vertex
[i
]) {
562 term
->vertex
[i
] = new evalue();
563 value_init(term
->vertex
[i
]->d
);
564 value_init(term
->vertex
[i
]->x
.n
);
565 zz2value(vertex
[i
], term
->vertex
[i
]->x
.n
);
566 value_set_si(term
->vertex
[i
]->d
, 1);
571 evalue_add_constant(term
->vertex
[i
], vertex
[i
]);
579 lex_order_rows(term
->den
);
582 void indicator_constructor::print(ostream
& os
, char **p
)
584 for (int i
= 0; i
< nbV
; ++i
)
585 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
586 os
<< "i: " << i
<< ", j: " << j
<< endl
;
587 terms
[i
][j
]->print(os
, p
);
592 void indicator_constructor::normalize()
594 for (int i
= 0; i
< nbV
; ++i
)
595 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
597 vertex
.SetLength(terms
[i
][j
]->den
.NumCols());
598 for (int r
= 0; r
< terms
[i
][j
]->den
.NumRows(); ++r
) {
599 for (int k
= 0; k
< terms
[i
][j
]->den
.NumCols(); ++k
) {
600 if (terms
[i
][j
]->den
[r
][k
] == 0)
602 if (terms
[i
][j
]->den
[r
][k
] > 0)
604 terms
[i
][j
]->sign
= -terms
[i
][j
]->sign
;
605 terms
[i
][j
]->den
[r
] = -terms
[i
][j
]->den
[r
];
606 vertex
+= terms
[i
][j
]->den
[r
];
610 lex_order_rows(terms
[i
][j
]->den
);
611 for (int k
= 0; k
< vertex
.length(); ++k
)
613 evalue_add_constant(terms
[i
][j
]->vertex
[k
], vertex
[k
]);
617 struct order_cache_el
{
619 order_cache_el
copy() const {
621 for (int i
= 0; i
< e
.size(); ++i
) {
622 evalue
*c
= new evalue
;
624 evalue_copy(c
, e
[i
]);
630 for (int i
= 0; i
< e
.size(); ++i
) {
631 free_evalue_refs(e
[i
]);
638 evalue_set_si(&mone
, -1, 1);
639 for (int i
= 0; i
< e
.size(); ++i
)
641 free_evalue_refs(&mone
);
643 void print(ostream
& os
, char **p
);
646 void order_cache_el::print(ostream
& os
, char **p
)
649 for (int i
= 0; i
< e
.size(); ++i
) {
652 evalue_print(os
, e
[i
], p
);
658 vector
<order_cache_el
> lt
;
659 vector
<order_cache_el
> le
;
660 vector
<order_cache_el
> unknown
;
662 void clear_transients() {
663 for (int i
= 0; i
< le
.size(); ++i
)
665 for (int i
= 0; i
< unknown
.size(); ++i
)
672 for (int i
= 0; i
< lt
.size(); ++i
)
676 void add(order_cache_el
& cache_el
, order_sign sign
);
677 order_sign
check_lt(vector
<order_cache_el
>* list
,
678 const indicator_term
*a
, const indicator_term
*b
,
679 order_cache_el
& cache_el
);
680 order_sign
check_lt(const indicator_term
*a
, const indicator_term
*b
,
681 order_cache_el
& cache_el
);
682 order_sign
check_direct(const indicator_term
*a
, const indicator_term
*b
,
683 order_cache_el
& cache_el
);
684 order_sign
check(const indicator_term
*a
, const indicator_term
*b
,
685 order_cache_el
& cache_el
);
686 void copy(const order_cache
& cache
);
687 void print(ostream
& os
, char **p
);
690 void order_cache::copy(const order_cache
& cache
)
692 for (int i
= 0; i
< cache
.lt
.size(); ++i
) {
693 order_cache_el n
= cache
.lt
[i
].copy();
698 void order_cache::add(order_cache_el
& cache_el
, order_sign sign
)
700 if (sign
== order_lt
) {
701 lt
.push_back(cache_el
);
702 } else if (sign
== order_gt
) {
704 lt
.push_back(cache_el
);
705 } else if (sign
== order_le
) {
706 le
.push_back(cache_el
);
707 } else if (sign
== order_ge
) {
709 le
.push_back(cache_el
);
710 } else if (sign
== order_unknown
) {
711 unknown
.push_back(cache_el
);
713 assert(sign
== order_eq
);
720 static evalue
*ediff(const evalue
*a
, const evalue
*b
)
724 evalue_set_si(&mone
, -1, 1);
725 evalue
*diff
= new evalue
;
727 evalue_copy(diff
, b
);
731 free_evalue_refs(&mone
);
735 static bool evalue_first_difference(const evalue
*e1
, const evalue
*e2
,
736 const evalue
**d1
, const evalue
**d2
)
741 if (value_ne(e1
->d
, e2
->d
))
744 if (value_notzero_p(e1
->d
)) {
745 if (value_eq(e1
->x
.n
, e2
->x
.n
))
749 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
751 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
753 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
756 assert(e1
->x
.p
->type
== polynomial
||
757 e1
->x
.p
->type
== fractional
||
758 e1
->x
.p
->type
== flooring
);
759 int offset
= type_offset(e1
->x
.p
);
760 assert(e1
->x
.p
->size
== offset
+2);
761 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
)
762 if (i
!= type_offset(e1
->x
.p
) &&
763 !eequal(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]))
766 return evalue_first_difference(&e1
->x
.p
->arr
[offset
],
767 &e2
->x
.p
->arr
[offset
], d1
, d2
);
770 static order_sign
evalue_diff_constant_sign(const evalue
*e1
, const evalue
*e2
)
772 if (!evalue_first_difference(e1
, e2
, &e1
, &e2
))
774 if (value_zero_p(e1
->d
) || value_zero_p(e2
->d
))
775 return order_undefined
;
776 int s
= evalue_rational_cmp(e1
, e2
);
785 order_sign
order_cache::check_lt(vector
<order_cache_el
>* list
,
786 const indicator_term
*a
, const indicator_term
*b
,
787 order_cache_el
& cache_el
)
789 order_sign sign
= order_undefined
;
790 for (int i
= 0; i
< list
->size(); ++i
) {
792 for (j
= cache_el
.e
.size(); j
< (*list
)[i
].e
.size(); ++j
)
793 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
794 for (j
= 0; j
< (*list
)[i
].e
.size(); ++j
) {
795 order_sign diff_sign
;
796 diff_sign
= evalue_diff_constant_sign((*list
)[i
].e
[j
], cache_el
.e
[j
]);
797 if (diff_sign
== order_gt
) {
800 } else if (diff_sign
== order_lt
)
802 else if (diff_sign
== order_undefined
)
805 assert(diff_sign
== order_eq
);
807 if (j
== (*list
)[i
].e
.size())
808 sign
= list
== <
? order_lt
: order_le
;
809 if (sign
!= order_undefined
)
815 order_sign
order_cache::check_direct(const indicator_term
*a
,
816 const indicator_term
*b
,
817 order_cache_el
& cache_el
)
819 order_sign sign
= check_lt(<
, a
, b
, cache_el
);
820 if (sign
!= order_undefined
)
822 sign
= check_lt(&le
, a
, b
, cache_el
);
823 if (sign
!= order_undefined
)
826 for (int i
= 0; i
< unknown
.size(); ++i
) {
828 for (j
= cache_el
.e
.size(); j
< unknown
[i
].e
.size(); ++j
)
829 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
830 for (j
= 0; j
< unknown
[i
].e
.size(); ++j
) {
831 if (!eequal(unknown
[i
].e
[j
], cache_el
.e
[j
]))
834 if (j
== unknown
[i
].e
.size()) {
835 sign
= order_unknown
;
842 order_sign
order_cache::check(const indicator_term
*a
, const indicator_term
*b
,
843 order_cache_el
& cache_el
)
845 order_sign sign
= check_direct(a
, b
, cache_el
);
846 if (sign
!= order_undefined
)
848 int size
= cache_el
.e
.size();
850 sign
= check_direct(a
, b
, cache_el
);
852 assert(cache_el
.e
.size() == size
);
853 if (sign
== order_undefined
)
855 if (sign
== order_lt
)
857 else if (sign
== order_le
)
860 assert(sign
== order_unknown
);
866 struct partial_order
{
869 std::set
<const indicator_term
*, smaller_it
> head
;
870 map
<const indicator_term
*, int, smaller_it
> pred
;
871 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> lt
;
872 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> le
;
873 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> eq
;
875 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> pending
;
879 partial_order(indicator
*ind
) : ind(ind
) {}
880 void copy(const partial_order
& order
,
881 map
< const indicator_term
*, indicator_term
* > old2new
);
883 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
884 map
<const indicator_term
*, int >::iterator j
;
885 std::set
<const indicator_term
*>::iterator k
;
887 if (head
.key_comp().requires_resort
) {
888 typeof(head
) new_head
;
889 for (k
= head
.begin(); k
!= head
.end(); ++k
)
895 if (pred
.key_comp().requires_resort
) {
896 typeof(pred
) new_pred
;
897 for (j
= pred
.begin(); j
!= pred
.end(); ++j
)
898 new_pred
[(*j
).first
] = (*j
).second
;
903 if (lt
.key_comp().requires_resort
) {
905 for (i
= lt
.begin(); i
!= lt
.end(); ++i
)
906 m
[(*i
).first
] = (*i
).second
;
911 if (le
.key_comp().requires_resort
) {
913 for (i
= le
.begin(); i
!= le
.end(); ++i
)
914 m
[(*i
).first
] = (*i
).second
;
919 if (eq
.key_comp().requires_resort
) {
921 for (i
= eq
.begin(); i
!= eq
.end(); ++i
)
922 m
[(*i
).first
] = (*i
).second
;
927 if (pending
.key_comp().requires_resort
) {
929 for (i
= pending
.begin(); i
!= pending
.end(); ++i
)
930 m
[(*i
).first
] = (*i
).second
;
936 order_sign
compare(const indicator_term
*a
, const indicator_term
*b
);
937 void set_equal(const indicator_term
*a
, const indicator_term
*b
);
938 void unset_le(const indicator_term
*a
, const indicator_term
*b
);
939 void dec_pred(const indicator_term
*it
) {
940 if (--pred
[it
] == 0) {
945 void inc_pred(const indicator_term
*it
) {
946 if (head
.find(it
) != head
.end())
951 bool compared(const indicator_term
* a
, const indicator_term
* b
);
952 void add(const indicator_term
* it
, std::set
<const indicator_term
*> *filter
);
953 void remove(const indicator_term
* it
);
955 void print(ostream
& os
, char **p
);
957 /* replace references to orig to references to replacement */
958 void replace(const indicator_term
* orig
, indicator_term
* replacement
);
959 void sanity_check() const;
962 /* We actually replace the contents of orig by that of replacement,
963 * but we have to be careful since replacing the content changes
964 * the order of orig in the maps.
966 void partial_order::replace(const indicator_term
* orig
, indicator_term
* replacement
)
968 std::set
<const indicator_term
*>::iterator k
;
970 bool is_head
= k
!= head
.end();
975 orig_pred
= pred
[orig
];
978 vector
<const indicator_term
* > orig_lt
;
979 vector
<const indicator_term
* > orig_le
;
980 vector
<const indicator_term
* > orig_eq
;
981 vector
<const indicator_term
* > orig_pending
;
982 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
983 bool in_lt
= ((i
= lt
.find(orig
)) != lt
.end());
985 orig_lt
= (*i
).second
;
988 bool in_le
= ((i
= le
.find(orig
)) != le
.end());
990 orig_le
= (*i
).second
;
993 bool in_eq
= ((i
= eq
.find(orig
)) != eq
.end());
995 orig_eq
= (*i
).second
;
998 bool in_pending
= ((i
= pending
.find(orig
)) != pending
.end());
1000 orig_pending
= (*i
).second
;
1001 pending
.erase(orig
);
1003 indicator_term
*old
= const_cast<indicator_term
*>(orig
);
1004 old
->swap(replacement
);
1008 pred
[old
] = orig_pred
;
1016 pending
[old
] = orig_pending
;
1019 void partial_order::unset_le(const indicator_term
*a
, const indicator_term
*b
)
1021 vector
<const indicator_term
*>::iterator i
;
1022 i
= find(le
[a
].begin(), le
[a
].end(), b
);
1024 if (le
[a
].size() == 0)
1027 i
= find(pending
[a
].begin(), pending
[a
].end(), b
);
1028 if (i
!= pending
[a
].end())
1029 pending
[a
].erase(i
);
1032 void partial_order::set_equal(const indicator_term
*a
, const indicator_term
*b
)
1034 if (eq
[a
].size() == 0)
1036 if (eq
[b
].size() == 0)
1041 if (pred
.key_comp()(b
, a
)) {
1042 const indicator_term
*c
= a
;
1047 const indicator_term
*base
= a
;
1049 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1051 for (int j
= 0; j
< eq
[b
].size(); ++j
) {
1052 eq
[base
].push_back(eq
[b
][j
]);
1053 eq
[eq
[b
][j
]][0] = base
;
1058 if (i
!= lt
.end()) {
1059 for (int j
= 0; j
< lt
[b
].size(); ++j
) {
1060 if (find(eq
[base
].begin(), eq
[base
].end(), lt
[b
][j
]) != eq
[base
].end())
1062 else if (find(lt
[base
].begin(), lt
[base
].end(), lt
[b
][j
])
1066 lt
[base
].push_back(lt
[b
][j
]);
1072 if (i
!= le
.end()) {
1073 for (int j
= 0; j
< le
[b
].size(); ++j
) {
1074 if (find(eq
[base
].begin(), eq
[base
].end(), le
[b
][j
]) != eq
[base
].end())
1076 else if (find(le
[base
].begin(), le
[base
].end(), le
[b
][j
])
1080 le
[base
].push_back(le
[b
][j
]);
1085 i
= pending
.find(base
);
1086 if (i
!= pending
.end()) {
1087 vector
<const indicator_term
* > old
= pending
[base
];
1088 pending
[base
].clear();
1089 for (int j
= 0; j
< old
.size(); ++j
) {
1090 if (find(eq
[base
].begin(), eq
[base
].end(), old
[j
]) == eq
[base
].end())
1091 pending
[base
].push_back(old
[j
]);
1095 i
= pending
.find(b
);
1096 if (i
!= pending
.end()) {
1097 for (int j
= 0; j
< pending
[b
].size(); ++j
) {
1098 if (find(eq
[base
].begin(), eq
[base
].end(), pending
[b
][j
]) == eq
[base
].end())
1099 pending
[base
].push_back(pending
[b
][j
]);
1105 void partial_order::copy(const partial_order
& order
,
1106 map
< const indicator_term
*, indicator_term
* > old2new
)
1108 cache
.copy(order
.cache
);
1110 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator i
;
1111 map
<const indicator_term
*, int >::const_iterator j
;
1112 std::set
<const indicator_term
*>::const_iterator k
;
1114 for (k
= order
.head
.begin(); k
!= order
.head
.end(); ++k
)
1115 head
.insert(old2new
[*k
]);
1117 for (j
= order
.pred
.begin(); j
!= order
.pred
.end(); ++j
)
1118 pred
[old2new
[(*j
).first
]] = (*j
).second
;
1120 for (i
= order
.lt
.begin(); i
!= order
.lt
.end(); ++i
) {
1121 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1122 lt
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1124 for (i
= order
.le
.begin(); i
!= order
.le
.end(); ++i
) {
1125 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1126 le
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1128 for (i
= order
.eq
.begin(); i
!= order
.eq
.end(); ++i
) {
1129 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1130 eq
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1132 for (i
= order
.pending
.begin(); i
!= order
.pending
.end(); ++i
) {
1133 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1134 pending
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1140 vector
<evalue
*> max
;
1142 void print(ostream
& os
, char **p
, barvinok_options
*options
) const;
1143 void substitute(Matrix
*T
, barvinok_options
*options
);
1144 Vector
*eval(Value
*val
, unsigned MaxRays
) const;
1147 for (int i
= 0; i
< max
.size(); ++i
) {
1148 free_evalue_refs(max
[i
]);
1156 * Project on first dim dimensions
1158 Polyhedron
* Polyhedron_Project_Initial(Polyhedron
*P
, int dim
)
1164 if (P
->Dimension
== dim
)
1165 return Polyhedron_Copy(P
);
1167 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
1168 for (i
= 0; i
< dim
; ++i
)
1169 value_set_si(T
->p
[i
][i
], 1);
1170 value_set_si(T
->p
[dim
][P
->Dimension
], 1);
1171 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
1177 vector
<indicator_term
*> term
;
1178 indicator_constructor
& ic
;
1179 partial_order order
;
1183 lexmin_options
*options
;
1184 vector
<evalue
*> substitutions
;
1186 indicator(indicator_constructor
& ic
, Param_Domain
*PD
, EDomain
*D
,
1187 lexmin_options
*options
) :
1188 ic(ic
), PD(PD
), D(D
), order(this), options(options
), P(NULL
) {}
1189 indicator(const indicator
& ind
, EDomain
*D
) :
1190 ic(ind
.ic
), PD(ind
.PD
), D(NULL
), order(this), options(ind
.options
),
1191 P(Polyhedron_Copy(ind
.P
)) {
1192 map
< const indicator_term
*, indicator_term
* > old2new
;
1193 for (int i
= 0; i
< ind
.term
.size(); ++i
) {
1194 indicator_term
*it
= new indicator_term(*ind
.term
[i
]);
1195 old2new
[ind
.term
[i
]] = it
;
1198 order
.copy(ind
.order
, old2new
);
1202 for (int i
= 0; i
< term
.size(); ++i
)
1210 void set_domain(EDomain
*D
) {
1211 order
.cache
.clear_transients();
1215 int nparam
= ic
.P
->Dimension
- ic
.vertex
.length();
1216 if (options
->reduce
) {
1217 Polyhedron
*Q
= Polyhedron_Project_Initial(D
->D
, nparam
);
1218 Q
= DomainConstraintSimplify(Q
, options
->barvinok
->MaxRays
);
1219 if (!P
|| !PolyhedronIncludes(Q
, P
))
1220 reduce_in_domain(Q
);
1228 void add(const indicator_term
* it
);
1229 void remove(const indicator_term
* it
);
1230 void remove_initial_rational_vertices();
1231 void expand_rational_vertex(const indicator_term
*initial
);
1233 void print(ostream
& os
, char **p
);
1235 void peel(int i
, int j
);
1236 void combine(const indicator_term
*a
, const indicator_term
*b
);
1237 void add_substitution(evalue
*equation
);
1238 void perform_pending_substitutions();
1239 void reduce_in_domain(Polyhedron
*D
);
1240 bool handle_equal_numerators(const indicator_term
*base
);
1242 max_term
* create_max_term(const indicator_term
*it
);
1244 void substitute(evalue
*equation
);
1247 void partial_order::sanity_check() const
1249 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator i
;
1250 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator prev
;
1251 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator l
;
1252 map
<const indicator_term
*, int >::const_iterator k
, prev_k
;
1254 for (k
= pred
.begin(); k
!= pred
.end(); prev_k
= k
, ++k
)
1255 if (k
!= pred
.begin())
1256 assert(pred
.key_comp()((*prev_k
).first
, (*k
).first
));
1257 for (i
= lt
.begin(); i
!= lt
.end(); prev
= i
, ++i
) {
1260 i_v
= (*i
).first
->eval(ind
->D
->sample
->p
);
1261 if (i
!= lt
.begin())
1262 assert(lt
.key_comp()((*prev
).first
, (*i
).first
));
1263 l
= eq
.find((*i
).first
);
1265 assert((*l
).second
.size() > 1);
1266 assert(head
.find((*i
).first
) != head
.end() ||
1267 pred
.find((*i
).first
) != pred
.end());
1268 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1269 k
= pred
.find((*i
).second
[j
]);
1270 assert(k
!= pred
.end());
1271 assert((*k
).second
!= 0);
1272 if ((*i
).first
->sign
!= 0 &&
1273 (*i
).second
[j
]->sign
!= 0 && ind
->D
->sample
) {
1274 vec_ZZ j_v
= (*i
).second
[j
]->eval(ind
->D
->sample
->p
);
1275 assert(lex_cmp(i_v
, j_v
) < 0);
1279 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1280 assert((*i
).second
.size() > 0);
1281 assert(head
.find((*i
).first
) != head
.end() ||
1282 pred
.find((*i
).first
) != pred
.end());
1283 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1284 k
= pred
.find((*i
).second
[j
]);
1285 assert(k
!= pred
.end());
1286 assert((*k
).second
!= 0);
1289 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1290 assert(head
.find((*i
).first
) != head
.end() ||
1291 pred
.find((*i
).first
) != pred
.end());
1292 assert((*i
).second
.size() >= 1);
1294 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1295 assert(head
.find((*i
).first
) != head
.end() ||
1296 pred
.find((*i
).first
) != pred
.end());
1297 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1298 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1299 pred
.find((*i
).second
[j
]) != pred
.end());
1303 max_term
* indicator::create_max_term(const indicator_term
*it
)
1305 int dim
= it
->den
.NumCols();
1306 int nparam
= ic
.P
->Dimension
- ic
.vertex
.length();
1307 max_term
*maximum
= new max_term
;
1308 maximum
->domain
= new EDomain(D
);
1309 for (int j
= 0; j
< dim
; ++j
) {
1310 evalue
*E
= new evalue
;
1312 evalue_copy(E
, it
->vertex
[j
]);
1313 if (evalue_frac2floor_in_domain3(E
, D
->D
, 0))
1315 maximum
->max
.push_back(E
);
1320 static order_sign
evalue_sign(evalue
*diff
, EDomain
*D
, lexmin_options
*options
)
1322 order_sign sign
= order_eq
;
1325 evalue_set_si(&mone
, -1, 1);
1326 int len
= 1 + D
->D
->Dimension
+ 1;
1327 Vector
*c
= Vector_Alloc(len
);
1328 Matrix
*T
= Matrix_Alloc(2, len
-1);
1330 int fract
= evalue2constraint(D
, diff
, c
->p
, len
);
1331 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1332 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1334 order_sign upper_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1335 if (upper_sign
== order_lt
|| !fract
)
1339 evalue2constraint(D
, diff
, c
->p
, len
);
1341 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1342 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1344 order_sign neg_lower_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1346 if (neg_lower_sign
== order_lt
)
1348 else if (neg_lower_sign
== order_eq
|| neg_lower_sign
== order_le
) {
1349 if (upper_sign
== order_eq
|| upper_sign
== order_le
)
1354 if (upper_sign
== order_lt
|| upper_sign
== order_le
||
1355 upper_sign
== order_eq
)
1358 sign
= order_unknown
;
1364 free_evalue_refs(&mone
);
1369 /* An auxiliary class that keeps a reference to an evalue
1370 * and frees it when it goes out of scope.
1372 struct temp_evalue
{
1374 temp_evalue() : E(NULL
) {}
1375 temp_evalue(evalue
*e
) : E(e
) {}
1376 operator evalue
* () const { return E
; }
1377 evalue
*operator=(evalue
*e
) {
1379 free_evalue_refs(E
);
1387 free_evalue_refs(E
);
1393 static void substitute(vector
<indicator_term
*>& term
, evalue
*equation
)
1395 evalue
*fract
= NULL
;
1396 evalue
*val
= new evalue
;
1398 evalue_copy(val
, equation
);
1401 value_init(factor
.d
);
1402 value_init(factor
.x
.n
);
1405 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= fractional
;
1406 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1409 if (value_zero_p(e
->d
) && e
->x
.p
->type
== fractional
)
1410 fract
= &e
->x
.p
->arr
[0];
1412 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= polynomial
;
1413 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1415 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== polynomial
);
1418 int offset
= type_offset(e
->x
.p
);
1420 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].d
));
1421 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].x
.n
));
1422 if (value_neg_p(e
->x
.p
->arr
[offset
+1].x
.n
)) {
1423 value_oppose(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1424 value_assign(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1426 value_assign(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1427 value_oppose(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1430 free_evalue_refs(&e
->x
.p
->arr
[offset
+1]);
1433 *e
= e
->x
.p
->arr
[offset
];
1438 for (int i
= 0; i
< term
.size(); ++i
)
1439 term
[i
]->substitute(fract
, val
);
1441 free_evalue_refs(&p
->arr
[0]);
1443 for (int i
= 0; i
< term
.size(); ++i
)
1444 term
[i
]->substitute(p
->pos
, val
);
1447 free_evalue_refs(&factor
);
1448 free_evalue_refs(val
);
1454 order_sign
partial_order::compare(const indicator_term
*a
, const indicator_term
*b
)
1456 unsigned dim
= a
->den
.NumCols();
1457 order_sign sign
= order_eq
;
1458 EDomain
*D
= ind
->D
;
1459 unsigned MaxRays
= ind
->options
->barvinok
->MaxRays
;
1460 bool rational
= a
->sign
== 0 || b
->sign
== 0;
1462 order_sign cached_sign
= order_eq
;
1463 for (int k
= 0; k
< dim
; ++k
) {
1464 cached_sign
= evalue_diff_constant_sign(a
->vertex
[k
], b
->vertex
[k
]);
1465 if (cached_sign
!= order_eq
)
1468 if (cached_sign
!= order_undefined
)
1471 order_cache_el cache_el
;
1472 cached_sign
= order_undefined
;
1474 cached_sign
= cache
.check(a
, b
, cache_el
);
1475 if (cached_sign
!= order_undefined
) {
1480 if (rational
&& POL_ISSET(ind
->options
->barvinok
->MaxRays
, POL_INTEGER
)) {
1481 ind
->options
->barvinok
->MaxRays
&= ~POL_INTEGER
;
1482 if (ind
->options
->barvinok
->MaxRays
)
1483 ind
->options
->barvinok
->MaxRays
|= POL_HIGH_BIT
;
1488 vector
<indicator_term
*> term
;
1490 for (int k
= 0; k
< dim
; ++k
) {
1491 /* compute a->vertex[k] - b->vertex[k] */
1493 if (cache_el
.e
.size() <= k
) {
1494 diff
= ediff(a
->vertex
[k
], b
->vertex
[k
]);
1495 cache_el
.e
.push_back(diff
);
1497 diff
= cache_el
.e
[k
];
1500 tdiff
= diff
= ediff(term
[0]->vertex
[k
], term
[1]->vertex
[k
]);
1501 order_sign diff_sign
;
1503 diff_sign
= order_undefined
;
1504 else if (eequal(a
->vertex
[k
], b
->vertex
[k
]))
1505 diff_sign
= order_eq
;
1507 diff_sign
= evalue_sign(diff
, D
, ind
->options
);
1509 if (diff_sign
== order_undefined
) {
1510 assert(sign
== order_le
|| sign
== order_ge
);
1511 if (sign
== order_le
)
1517 if (diff_sign
== order_lt
) {
1518 if (sign
== order_eq
|| sign
== order_le
)
1521 sign
= order_unknown
;
1524 if (diff_sign
== order_gt
) {
1525 if (sign
== order_eq
|| sign
== order_ge
)
1528 sign
= order_unknown
;
1531 if (diff_sign
== order_eq
) {
1532 if (D
== ind
->D
&& term
.size() == 0 && !rational
&&
1533 !EVALUE_IS_ZERO(*diff
))
1534 ind
->add_substitution(diff
);
1537 if ((diff_sign
== order_unknown
) ||
1538 ((diff_sign
== order_lt
|| diff_sign
== order_le
) && sign
== order_ge
) ||
1539 ((diff_sign
== order_gt
|| diff_sign
== order_ge
) && sign
== order_le
)) {
1540 sign
= order_unknown
;
1547 term
.push_back(new indicator_term(*a
));
1548 term
.push_back(new indicator_term(*b
));
1550 substitute(term
, diff
);
1554 cache
.add(cache_el
, sign
);
1558 if (D
&& D
!= ind
->D
)
1566 ind
->options
->barvinok
->MaxRays
= MaxRays
;
1570 bool partial_order::compared(const indicator_term
* a
, const indicator_term
* b
)
1572 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator j
;
1575 if (j
!= lt
.end() && find(lt
[a
].begin(), lt
[a
].end(), b
) != lt
[a
].end())
1579 if (j
!= le
.end() && find(le
[a
].begin(), le
[a
].end(), b
) != le
[a
].end())
1585 void partial_order::add(const indicator_term
* it
,
1586 std::set
<const indicator_term
*> *filter
)
1588 if (eq
.find(it
) != eq
.end() && eq
[it
].size() == 1)
1591 typeof(head
) head_copy(head
);
1596 std::set
<const indicator_term
*>::iterator i
;
1597 for (i
= head_copy
.begin(); i
!= head_copy
.end(); ++i
) {
1600 if (eq
.find(*i
) != eq
.end() && eq
[*i
].size() == 1)
1603 if (filter
->find(*i
) == filter
->end())
1605 if (compared(*i
, it
))
1608 order_sign sign
= compare(it
, *i
);
1609 if (sign
== order_lt
) {
1610 lt
[it
].push_back(*i
);
1612 } else if (sign
== order_le
) {
1613 le
[it
].push_back(*i
);
1615 } else if (sign
== order_eq
) {
1618 } else if (sign
== order_gt
) {
1619 pending
[*i
].push_back(it
);
1620 lt
[*i
].push_back(it
);
1622 } else if (sign
== order_ge
) {
1623 pending
[*i
].push_back(it
);
1624 le
[*i
].push_back(it
);
1630 void partial_order::remove(const indicator_term
* it
)
1632 std::set
<const indicator_term
*> filter
;
1633 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1635 assert(head
.find(it
) != head
.end());
1638 if (i
!= eq
.end()) {
1639 assert(eq
[it
].size() >= 1);
1640 const indicator_term
*base
;
1641 if (eq
[it
].size() == 1) {
1645 vector
<const indicator_term
* >::iterator j
;
1646 j
= find(eq
[base
].begin(), eq
[base
].end(), it
);
1647 assert(j
!= eq
[base
].end());
1650 /* "it" may no longer be the smallest, since the order
1651 * structure may have been copied from another one.
1653 sort(eq
[it
].begin()+1, eq
[it
].end(), pred
.key_comp());
1654 assert(eq
[it
][0] == it
);
1655 eq
[it
].erase(eq
[it
].begin());
1660 for (int j
= 1; j
< eq
[base
].size(); ++j
)
1661 eq
[eq
[base
][j
]][0] = base
;
1664 if (i
!= lt
.end()) {
1670 if (i
!= le
.end()) {
1675 i
= pending
.find(it
);
1676 if (i
!= pending
.end()) {
1677 pending
[base
] = pending
[it
];
1682 if (eq
[base
].size() == 1)
1691 if (i
!= lt
.end()) {
1692 for (int j
= 0; j
< lt
[it
].size(); ++j
) {
1693 filter
.insert(lt
[it
][j
]);
1694 dec_pred(lt
[it
][j
]);
1700 if (i
!= le
.end()) {
1701 for (int j
= 0; j
< le
[it
].size(); ++j
) {
1702 filter
.insert(le
[it
][j
]);
1703 dec_pred(le
[it
][j
]);
1710 i
= pending
.find(it
);
1711 if (i
!= pending
.end()) {
1712 vector
<const indicator_term
*> it_pending
= pending
[it
];
1714 for (int j
= 0; j
< it_pending
.size(); ++j
) {
1715 filter
.erase(it_pending
[j
]);
1716 add(it_pending
[j
], &filter
);
1721 void partial_order::print(ostream
& os
, char **p
)
1723 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1724 map
<const indicator_term
*, int >::iterator j
;
1725 std::set
<const indicator_term
*>::iterator k
;
1726 for (k
= head
.begin(); k
!= head
.end(); ++k
) {
1730 for (j
= pred
.begin(); j
!= pred
.end(); ++j
) {
1731 (*j
).first
->print(os
, p
);
1732 os
<< ": " << (*j
).second
<< endl
;
1734 for (i
= lt
.begin(); i
!= lt
.end(); ++i
) {
1735 (*i
).first
->print(os
, p
);
1736 assert(head
.find((*i
).first
) != head
.end() ||
1737 pred
.find((*i
).first
) != pred
.end());
1738 if (pred
.find((*i
).first
) != pred
.end())
1739 os
<< "(" << pred
[(*i
).first
] << ")";
1741 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1744 (*i
).second
[j
]->print(os
, p
);
1745 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1746 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1750 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1751 (*i
).first
->print(os
, p
);
1752 assert(head
.find((*i
).first
) != head
.end() ||
1753 pred
.find((*i
).first
) != pred
.end());
1754 if (pred
.find((*i
).first
) != pred
.end())
1755 os
<< "(" << pred
[(*i
).first
] << ")";
1757 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1760 (*i
).second
[j
]->print(os
, p
);
1761 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1762 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1766 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1767 if ((*i
).second
.size() <= 1)
1769 (*i
).first
->print(os
, p
);
1770 assert(head
.find((*i
).first
) != head
.end() ||
1771 pred
.find((*i
).first
) != pred
.end());
1772 if (pred
.find((*i
).first
) != pred
.end())
1773 os
<< "(" << pred
[(*i
).first
] << ")";
1774 for (int j
= 1; j
< (*i
).second
.size(); ++j
) {
1777 (*i
).second
[j
]->print(os
, p
);
1778 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1779 pred
.find((*i
).second
[j
]) != pred
.end());
1780 if (pred
.find((*i
).second
[j
]) != pred
.end())
1781 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1785 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1786 os
<< "pending on ";
1787 (*i
).first
->print(os
, p
);
1788 assert(head
.find((*i
).first
) != head
.end() ||
1789 pred
.find((*i
).first
) != pred
.end());
1790 if (pred
.find((*i
).first
) != pred
.end())
1791 os
<< "(" << pred
[(*i
).first
] << ")";
1793 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1796 (*i
).second
[j
]->print(os
, p
);
1797 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1798 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1804 void indicator::add(const indicator_term
* it
)
1806 indicator_term
*nt
= new indicator_term(*it
);
1807 if (options
->reduce
)
1808 nt
->reduce_in_domain(P
? P
: D
->D
);
1810 order
.add(nt
, NULL
);
1811 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1814 void indicator::remove(const indicator_term
* it
)
1816 vector
<indicator_term
*>::iterator i
;
1817 i
= find(term
.begin(), term
.end(), it
);
1818 assert(i
!= term
.end());
1821 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1825 void indicator::expand_rational_vertex(const indicator_term
*initial
)
1827 int pos
= initial
->pos
;
1829 if (ic
.terms
[pos
].size() == 0) {
1831 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, ic
.PP
) // _i is internal counter
1833 ic
.decompose_at_vertex(V
, pos
, options
->barvinok
);
1836 END_FORALL_PVertex_in_ParamPolyhedron
;
1838 for (int j
= 0; j
< ic
.terms
[pos
].size(); ++j
)
1839 add(ic
.terms
[pos
][j
]);
1842 void indicator::remove_initial_rational_vertices()
1845 const indicator_term
*initial
= NULL
;
1846 std::set
<const indicator_term
*>::iterator i
;
1847 for (i
= order
.head
.begin(); i
!= order
.head
.end(); ++i
) {
1848 if ((*i
)->sign
!= 0)
1850 if (order
.eq
.find(*i
) != order
.eq
.end() && order
.eq
[*i
].size() <= 1)
1857 expand_rational_vertex(initial
);
1861 void indicator::reduce_in_domain(Polyhedron
*D
)
1863 for (int i
= 0; i
< term
.size(); ++i
)
1864 term
[i
]->reduce_in_domain(D
);
1867 void indicator::print(ostream
& os
, char **p
)
1869 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1870 for (int i
= 0; i
< term
.size(); ++i
) {
1871 term
[i
]->print(os
, p
);
1873 os
<< ": " << term
[i
]->eval(D
->sample
->p
);
1880 /* Remove pairs of opposite terms */
1881 void indicator::simplify()
1883 for (int i
= 0; i
< term
.size(); ++i
) {
1884 for (int j
= i
+1; j
< term
.size(); ++j
) {
1885 if (term
[i
]->sign
+ term
[j
]->sign
!= 0)
1887 if (term
[i
]->den
!= term
[j
]->den
)
1890 for (k
= 0; k
< term
[i
]->den
.NumCols(); ++k
)
1891 if (!eequal(term
[i
]->vertex
[k
], term
[j
]->vertex
[k
]))
1893 if (k
< term
[i
]->den
.NumCols())
1897 term
.erase(term
.begin()+j
);
1898 term
.erase(term
.begin()+i
);
1905 void indicator::peel(int i
, int j
)
1913 int dim
= term
[i
]->den
.NumCols();
1918 int n_common
= 0, n_i
= 0, n_j
= 0;
1920 common
.SetDims(min(term
[i
]->den
.NumRows(), term
[j
]->den
.NumRows()), dim
);
1921 rest_i
.SetDims(term
[i
]->den
.NumRows(), dim
);
1922 rest_j
.SetDims(term
[j
]->den
.NumRows(), dim
);
1925 for (k
= 0, l
= 0; k
< term
[i
]->den
.NumRows() && l
< term
[j
]->den
.NumRows(); ) {
1926 int s
= lex_cmp(term
[i
]->den
[k
], term
[j
]->den
[l
]);
1928 common
[n_common
++] = term
[i
]->den
[k
];
1932 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1934 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1936 while (k
< term
[i
]->den
.NumRows())
1937 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1938 while (l
< term
[j
]->den
.NumRows())
1939 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1940 common
.SetDims(n_common
, dim
);
1941 rest_i
.SetDims(n_i
, dim
);
1942 rest_j
.SetDims(n_j
, dim
);
1944 for (k
= 0; k
<= n_i
; ++k
) {
1945 indicator_term
*it
= new indicator_term(*term
[i
]);
1946 it
->den
.SetDims(n_common
+ k
, dim
);
1947 for (l
= 0; l
< n_common
; ++l
)
1948 it
->den
[l
] = common
[l
];
1949 for (l
= 0; l
< k
; ++l
)
1950 it
->den
[n_common
+l
] = rest_i
[l
];
1951 lex_order_rows(it
->den
);
1953 for (l
= 0; l
< dim
; ++l
)
1954 evalue_add_constant(it
->vertex
[l
], rest_i
[k
-1][l
]);
1958 for (k
= 0; k
<= n_j
; ++k
) {
1959 indicator_term
*it
= new indicator_term(*term
[j
]);
1960 it
->den
.SetDims(n_common
+ k
, dim
);
1961 for (l
= 0; l
< n_common
; ++l
)
1962 it
->den
[l
] = common
[l
];
1963 for (l
= 0; l
< k
; ++l
)
1964 it
->den
[n_common
+l
] = rest_j
[l
];
1965 lex_order_rows(it
->den
);
1967 for (l
= 0; l
< dim
; ++l
)
1968 evalue_add_constant(it
->vertex
[l
], rest_j
[k
-1][l
]);
1971 term
.erase(term
.begin()+j
);
1972 term
.erase(term
.begin()+i
);
1975 void indicator::combine(const indicator_term
*a
, const indicator_term
*b
)
1977 int dim
= a
->den
.NumCols();
1980 mat_ZZ rest_i
; /* factors in a, but not in b */
1981 mat_ZZ rest_j
; /* factors in b, but not in a */
1982 int n_common
= 0, n_i
= 0, n_j
= 0;
1984 common
.SetDims(min(a
->den
.NumRows(), b
->den
.NumRows()), dim
);
1985 rest_i
.SetDims(a
->den
.NumRows(), dim
);
1986 rest_j
.SetDims(b
->den
.NumRows(), dim
);
1989 for (k
= 0, l
= 0; k
< a
->den
.NumRows() && l
< b
->den
.NumRows(); ) {
1990 int s
= lex_cmp(a
->den
[k
], b
->den
[l
]);
1992 common
[n_common
++] = a
->den
[k
];
1996 rest_i
[n_i
++] = a
->den
[k
++];
1998 rest_j
[n_j
++] = b
->den
[l
++];
2000 while (k
< a
->den
.NumRows())
2001 rest_i
[n_i
++] = a
->den
[k
++];
2002 while (l
< b
->den
.NumRows())
2003 rest_j
[n_j
++] = b
->den
[l
++];
2004 common
.SetDims(n_common
, dim
);
2005 rest_i
.SetDims(n_i
, dim
);
2006 rest_j
.SetDims(n_j
, dim
);
2008 assert(order
.eq
[a
].size() > 1);
2009 indicator_term
*prev
;
2012 for (int k
= n_i
-1; k
>= 0; --k
) {
2013 indicator_term
*it
= new indicator_term(*b
);
2014 it
->den
.SetDims(n_common
+ n_j
+ n_i
-k
, dim
);
2015 for (int l
= k
; l
< n_i
; ++l
)
2016 it
->den
[n_common
+n_j
+l
-k
] = rest_i
[l
];
2017 lex_order_rows(it
->den
);
2018 for (int m
= 0; m
< dim
; ++m
)
2019 evalue_add_constant(it
->vertex
[m
], rest_i
[k
][m
]);
2020 it
->sign
= -it
->sign
;
2022 order
.pending
[it
].push_back(prev
);
2023 order
.lt
[it
].push_back(prev
);
2024 order
.inc_pred(prev
);
2027 order
.head
.insert(it
);
2031 indicator_term
*it
= new indicator_term(*b
);
2032 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
2033 for (l
= 0; l
< n_i
; ++l
)
2034 it
->den
[n_common
+n_j
+l
] = rest_i
[l
];
2035 lex_order_rows(it
->den
);
2037 order
.pending
[a
].push_back(prev
);
2038 order
.lt
[a
].push_back(prev
);
2039 order
.inc_pred(prev
);
2040 order
.replace(b
, it
);
2045 for (int k
= n_j
-1; k
>= 0; --k
) {
2046 indicator_term
*it
= new indicator_term(*a
);
2047 it
->den
.SetDims(n_common
+ n_i
+ n_j
-k
, dim
);
2048 for (int l
= k
; l
< n_j
; ++l
)
2049 it
->den
[n_common
+n_i
+l
-k
] = rest_j
[l
];
2050 lex_order_rows(it
->den
);
2051 for (int m
= 0; m
< dim
; ++m
)
2052 evalue_add_constant(it
->vertex
[m
], rest_j
[k
][m
]);
2053 it
->sign
= -it
->sign
;
2055 order
.pending
[it
].push_back(prev
);
2056 order
.lt
[it
].push_back(prev
);
2057 order
.inc_pred(prev
);
2060 order
.head
.insert(it
);
2064 indicator_term
*it
= new indicator_term(*a
);
2065 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
2066 for (l
= 0; l
< n_j
; ++l
)
2067 it
->den
[n_common
+n_i
+l
] = rest_j
[l
];
2068 lex_order_rows(it
->den
);
2070 order
.pending
[a
].push_back(prev
);
2071 order
.lt
[a
].push_back(prev
);
2072 order
.inc_pred(prev
);
2073 order
.replace(a
, it
);
2077 assert(term
.size() == order
.head
.size() + order
.pred
.size());
2080 bool indicator::handle_equal_numerators(const indicator_term
*base
)
2082 for (int i
= 0; i
< order
.eq
[base
].size(); ++i
) {
2083 for (int j
= i
+1; j
< order
.eq
[base
].size(); ++j
) {
2084 if (order
.eq
[base
][i
]->is_opposite(order
.eq
[base
][j
])) {
2085 remove(order
.eq
[base
][j
]);
2086 remove(i
? order
.eq
[base
][i
] : base
);
2091 for (int j
= 1; j
< order
.eq
[base
].size(); ++j
)
2092 if (order
.eq
[base
][j
]->sign
!= base
->sign
) {
2093 combine(base
, order
.eq
[base
][j
]);
2099 void indicator::substitute(evalue
*equation
)
2101 ::substitute(term
, equation
);
2104 void indicator::add_substitution(evalue
*equation
)
2106 for (int i
= 0; i
< substitutions
.size(); ++i
)
2107 if (eequal(substitutions
[i
], equation
))
2109 evalue
*copy
= new evalue();
2110 value_init(copy
->d
);
2111 evalue_copy(copy
, equation
);
2112 substitutions
.push_back(copy
);
2115 void indicator::perform_pending_substitutions()
2117 if (substitutions
.size() == 0)
2120 for (int i
= 0; i
< substitutions
.size(); ++i
) {
2121 substitute(substitutions
[i
]);
2122 free_evalue_refs(substitutions
[i
]);
2123 delete substitutions
[i
];
2125 substitutions
.clear();
2129 static void print_varlist(ostream
& os
, int n
, char **names
)
2133 for (i
= 0; i
< n
; ++i
) {
2141 void max_term::print(ostream
& os
, char **p
, barvinok_options
*options
) const
2144 print_varlist(os
, domain
->dimension(), p
);
2147 for (int i
= 0; i
< max
.size(); ++i
) {
2150 evalue_print(os
, max
[i
], p
);
2154 domain
->print_constraints(os
, p
, options
);
2158 Matrix
*left_inverse(Matrix
*M
, Matrix
**Eq
)
2161 Matrix
*L
, *H
, *Q
, *U
, *ratH
, *invH
, *Ut
, *inv
;
2166 L
= Matrix_Alloc(M
->NbRows
-1, M
->NbColumns
-1);
2167 for (i
= 0; i
< L
->NbRows
; ++i
)
2168 Vector_Copy(M
->p
[i
], L
->p
[i
], L
->NbColumns
);
2169 right_hermite(L
, &H
, &U
, &Q
);
2172 t
= Vector_Alloc(U
->NbColumns
);
2173 for (i
= 0; i
< U
->NbColumns
; ++i
)
2174 value_oppose(t
->p
[i
], M
->p
[i
][M
->NbColumns
-1]);
2176 *Eq
= Matrix_Alloc(H
->NbRows
- H
->NbColumns
, 2 + U
->NbColumns
);
2177 for (i
= 0; i
< H
->NbRows
- H
->NbColumns
; ++i
) {
2178 Vector_Copy(U
->p
[H
->NbColumns
+i
], (*Eq
)->p
[i
]+1, U
->NbColumns
);
2179 Inner_Product(U
->p
[H
->NbColumns
+i
], t
->p
, U
->NbColumns
,
2180 (*Eq
)->p
[i
]+1+U
->NbColumns
);
2183 ratH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
2184 invH
= Matrix_Alloc(H
->NbColumns
+1, H
->NbColumns
+1);
2185 for (i
= 0; i
< H
->NbColumns
; ++i
)
2186 Vector_Copy(H
->p
[i
], ratH
->p
[i
], H
->NbColumns
);
2187 value_set_si(ratH
->p
[ratH
->NbRows
-1][ratH
->NbColumns
-1], 1);
2189 ok
= Matrix_Inverse(ratH
, invH
);
2192 Ut
= Matrix_Alloc(invH
->NbRows
, U
->NbColumns
+1);
2193 for (i
= 0; i
< Ut
->NbRows
-1; ++i
) {
2194 Vector_Copy(U
->p
[i
], Ut
->p
[i
], U
->NbColumns
);
2195 Inner_Product(U
->p
[i
], t
->p
, U
->NbColumns
, &Ut
->p
[i
][Ut
->NbColumns
-1]);
2199 value_set_si(Ut
->p
[Ut
->NbRows
-1][Ut
->NbColumns
-1], 1);
2200 inv
= Matrix_Alloc(invH
->NbRows
, Ut
->NbColumns
);
2201 Matrix_Product(invH
, Ut
, inv
);
2207 /* T maps the compressed parameters to the original parameters,
2208 * while this max_term is based on the compressed parameters
2209 * and we want get the original parameters back.
2211 void max_term::substitute(Matrix
*T
, barvinok_options
*options
)
2213 assert(domain
->dimension() == T
->NbColumns
-1);
2214 int nexist
= domain
->D
->Dimension
- (T
->NbColumns
-1);
2216 Matrix
*inv
= left_inverse(T
, &Eq
);
2219 value_init(denom
.d
);
2220 value_init(denom
.x
.n
);
2221 value_set_si(denom
.x
.n
, 1);
2222 value_assign(denom
.d
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
2225 v
.SetLength(inv
->NbColumns
);
2226 evalue
* subs
[inv
->NbRows
-1];
2227 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2228 values2zz(inv
->p
[i
], v
, v
.length());
2229 subs
[i
] = multi_monom(v
);
2230 emul(&denom
, subs
[i
]);
2232 free_evalue_refs(&denom
);
2234 domain
->substitute(subs
, inv
, Eq
, options
->MaxRays
);
2237 for (int i
= 0; i
< max
.size(); ++i
) {
2238 evalue_substitute(max
[i
], subs
);
2239 reduce_evalue(max
[i
]);
2242 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2243 free_evalue_refs(subs
[i
]);
2249 int Last_Non_Zero(Value
*p
, unsigned len
)
2251 for (int i
= len
-1; i
>= 0; --i
)
2252 if (value_notzero_p(p
[i
]))
2257 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
2259 for (int r
= 0; r
< n
; ++r
)
2260 value_swap(V
[r
][i
], V
[r
][j
]);
2263 static void SwapColumns(Polyhedron
*P
, int i
, int j
)
2265 SwapColumns(P
->Constraint
, P
->NbConstraints
, i
, j
);
2266 SwapColumns(P
->Ray
, P
->NbRays
, i
, j
);
2269 Vector
*max_term::eval(Value
*val
, unsigned MaxRays
) const
2271 if (!domain
->contains(val
, domain
->dimension()))
2273 Vector
*res
= Vector_Alloc(max
.size());
2274 for (int i
= 0; i
< max
.size(); ++i
) {
2275 compute_evalue(max
[i
], val
, &res
->p
[i
]);
2282 enum sign
{ le
, ge
, lge
} sign
;
2284 split (evalue
*c
, enum sign s
) : constraint(c
), sign(s
) {}
2287 static void split_on(const split
& sp
, EDomain
*D
,
2288 EDomain
**Dlt
, EDomain
**Deq
, EDomain
**Dgt
,
2289 lexmin_options
*options
)
2295 ge_constraint
*ge
= D
->compute_ge_constraint(sp
.constraint
);
2296 if (sp
.sign
== split::lge
|| sp
.sign
== split::ge
)
2297 ED
[2] = EDomain::new_from_ge_constraint(ge
, 1, options
->barvinok
);
2300 if (sp
.sign
== split::lge
|| sp
.sign
== split::le
)
2301 ED
[0] = EDomain::new_from_ge_constraint(ge
, -1, options
->barvinok
);
2305 assert(sp
.sign
== split::lge
|| sp
.sign
== split::ge
|| sp
.sign
== split::le
);
2306 ED
[1] = EDomain::new_from_ge_constraint(ge
, 0, options
->barvinok
);
2310 for (int i
= 0; i
< 3; ++i
) {
2313 if (D
->sample
&& ED
[i
]->contains(D
->sample
->p
, D
->sample
->Size
-1)) {
2314 ED
[i
]->sample
= Vector_Alloc(D
->sample
->Size
);
2315 Vector_Copy(D
->sample
->p
, ED
[i
]->sample
->p
, D
->sample
->Size
);
2316 } else if (emptyQ2(ED
[i
]->D
) ||
2317 (options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2318 !(ED
[i
]->not_empty(options
)))) {
2328 ostream
& operator<< (ostream
& os
, const vector
<int> & v
)
2331 for (int i
= 0; i
< v
.size(); ++i
) {
2340 static bool isTranslation(Matrix
*M
)
2343 if (M
->NbRows
!= M
->NbColumns
)
2346 for (i
= 0;i
< M
->NbRows
; i
++)
2347 for (j
= 0; j
< M
->NbColumns
-1; j
++)
2349 if(value_notone_p(M
->p
[i
][j
]))
2352 if(value_notzero_p(M
->p
[i
][j
]))
2355 return value_one_p(M
->p
[M
->NbRows
-1][M
->NbColumns
-1]);
2358 static Matrix
*compress_parameters(Polyhedron
**P
, Polyhedron
**C
,
2359 unsigned nparam
, unsigned MaxRays
)
2363 /* compress_parms doesn't like equalities that only involve parameters */
2364 for (int i
= 0; i
< (*P
)->NbEq
; ++i
)
2365 assert(First_Non_Zero((*P
)->Constraint
[i
]+1, (*P
)->Dimension
-nparam
) != -1);
2367 M
= Matrix_Alloc((*P
)->NbEq
, (*P
)->Dimension
+2);
2368 Vector_Copy((*P
)->Constraint
[0], M
->p
[0], (*P
)->NbEq
* ((*P
)->Dimension
+2));
2369 CP
= compress_parms(M
, nparam
);
2372 if (isTranslation(CP
)) {
2377 T
= align_matrix(CP
, (*P
)->Dimension
+1);
2378 *P
= Polyhedron_Preimage(*P
, T
, MaxRays
);
2381 *C
= Polyhedron_Preimage(*C
, CP
, MaxRays
);
2386 void construct_rational_vertices(Param_Polyhedron
*PP
, Matrix
*T
, unsigned dim
,
2387 int nparam
, vector
<indicator_term
*>& vertices
)
2396 v
.SetLength(nparam
+1);
2399 value_init(factor
.d
);
2400 value_init(factor
.x
.n
);
2401 value_set_si(factor
.x
.n
, 1);
2402 value_set_si(factor
.d
, 1);
2404 for (i
= 0, PV
= PP
->V
; PV
; ++i
, PV
= PV
->next
) {
2405 indicator_term
*term
= new indicator_term(dim
, i
);
2406 vertices
.push_back(term
);
2407 Matrix
*M
= Matrix_Alloc(PV
->Vertex
->NbRows
+nparam
+1, nparam
+1);
2408 value_set_si(lcm
, 1);
2409 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
)
2410 value_lcm(lcm
, PV
->Vertex
->p
[j
][nparam
+1], &lcm
);
2411 value_assign(M
->p
[M
->NbRows
-1][M
->NbColumns
-1], lcm
);
2412 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
) {
2413 value_division(tmp
, lcm
, PV
->Vertex
->p
[j
][nparam
+1]);
2414 Vector_Scale(PV
->Vertex
->p
[j
], M
->p
[j
], tmp
, nparam
+1);
2416 for (int j
= 0; j
< nparam
; ++j
)
2417 value_assign(M
->p
[PV
->Vertex
->NbRows
+j
][j
], lcm
);
2419 Matrix
*M2
= Matrix_Alloc(T
->NbRows
, M
->NbColumns
);
2420 Matrix_Product(T
, M
, M2
);
2424 for (int j
= 0; j
< dim
; ++j
) {
2425 values2zz(M
->p
[j
], v
, nparam
+1);
2426 term
->vertex
[j
] = multi_monom(v
);
2427 value_assign(factor
.d
, lcm
);
2428 emul(&factor
, term
->vertex
[j
]);
2432 assert(i
== PP
->nbV
);
2433 free_evalue_refs(&factor
);
2438 static vector
<max_term
*> lexmin(indicator
& ind
, unsigned nparam
,
2441 vector
<max_term
*> maxima
;
2442 std::set
<const indicator_term
*>::iterator i
;
2443 vector
<int> best_score
;
2444 vector
<int> second_score
;
2445 vector
<int> neg_score
;
2448 ind
.perform_pending_substitutions();
2449 const indicator_term
*best
= NULL
, *second
= NULL
, *neg
= NULL
,
2450 *neg_eq
= NULL
, *neg_le
= NULL
;
2451 for (i
= ind
.order
.head
.begin(); i
!= ind
.order
.head
.end(); ++i
) {
2453 const indicator_term
*term
= *i
;
2454 if (term
->sign
== 0) {
2455 ind
.expand_rational_vertex(term
);
2459 if (ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2461 if (ind
.order
.eq
[term
].size() <= 1)
2463 for (j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2464 if (ind
.order
.pred
.find(ind
.order
.eq
[term
][j
]) !=
2465 ind
.order
.pred
.end())
2467 if (j
< ind
.order
.eq
[term
].size())
2469 score
.push_back(ind
.order
.eq
[term
].size());
2472 if (ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2473 score
.push_back(ind
.order
.le
[term
].size());
2476 if (ind
.order
.lt
.find(term
) != ind
.order
.lt
.end())
2477 score
.push_back(-ind
.order
.lt
[term
].size());
2481 if (term
->sign
> 0) {
2482 if (!best
|| score
< best_score
) {
2484 second_score
= best_score
;
2487 } else if (!second
|| score
< second_score
) {
2489 second_score
= score
;
2492 if (!neg_eq
&& ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2493 for (int j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2494 if (ind
.order
.eq
[term
][j
]->sign
!= term
->sign
) {
2499 if (!neg_le
&& ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2501 if (!neg
|| score
< neg_score
) {
2507 if (i
!= ind
.order
.head
.end())
2510 if (!best
&& neg_eq
) {
2511 assert(ind
.order
.eq
[neg_eq
].size() != 0);
2512 bool handled
= ind
.handle_equal_numerators(neg_eq
);
2517 if (!best
&& neg_le
) {
2518 /* The smallest term is negative and <= some positive term */
2524 /* apparently there can be negative initial term on empty domains */
2525 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2526 ind
.options
->polysign
== BV_LEXMIN_POLYSIGN_POLYLIB
)
2531 if (!second
&& !neg
) {
2532 const indicator_term
*rat
= NULL
;
2534 if (ind
.order
.le
.find(best
) == ind
.order
.le
.end()) {
2535 if (ind
.order
.eq
.find(best
) != ind
.order
.eq
.end()) {
2536 bool handled
= ind
.handle_equal_numerators(best
);
2537 if (ind
.options
->emptiness_check
!=
2538 BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2539 ind
.options
->polysign
== BV_LEXMIN_POLYSIGN_POLYLIB
)
2541 /* If !handled then the leading coefficient is bigger than one;
2542 * must be an empty domain
2549 maxima
.push_back(ind
.create_max_term(best
));
2552 for (int j
= 0; j
< ind
.order
.le
[best
].size(); ++j
) {
2553 if (ind
.order
.le
[best
][j
]->sign
== 0) {
2554 if (!rat
&& ind
.order
.pred
[ind
.order
.le
[best
][j
]] == 1)
2555 rat
= ind
.order
.le
[best
][j
];
2556 } else if (ind
.order
.le
[best
][j
]->sign
> 0) {
2557 second
= ind
.order
.le
[best
][j
];
2560 neg
= ind
.order
.le
[best
][j
];
2563 if (!second
&& !neg
) {
2565 ind
.order
.unset_le(best
, rat
);
2566 ind
.expand_rational_vertex(rat
);
2573 ind
.order
.unset_le(best
, second
);
2579 unsigned dim
= best
->den
.NumCols();
2582 for (int k
= 0; k
< dim
; ++k
) {
2583 diff
= ediff(best
->vertex
[k
], second
->vertex
[k
]);
2584 sign
= evalue_sign(diff
, ind
.D
, ind
.options
);
2586 /* neg can never be smaller than best, unless it may still cancel.
2587 * This can happen if positive terms have been determined to
2588 * be equal or less than or equal to some negative term.
2590 if (second
== neg
&& !neg_eq
&& !neg_le
) {
2591 if (sign
== order_ge
)
2593 if (sign
== order_unknown
)
2597 if (sign
!= order_eq
)
2599 if (!EVALUE_IS_ZERO(*diff
)) {
2600 ind
.add_substitution(diff
);
2601 ind
.perform_pending_substitutions();
2604 if (sign
== order_eq
) {
2605 ind
.order
.set_equal(best
, second
);
2608 if (sign
== order_lt
) {
2609 ind
.order
.lt
[best
].push_back(second
);
2610 ind
.order
.inc_pred(second
);
2613 if (sign
== order_gt
) {
2614 ind
.order
.lt
[second
].push_back(best
);
2615 ind
.order
.inc_pred(best
);
2619 split
sp(diff
, sign
== order_le
? split::le
:
2620 sign
== order_ge
? split::ge
: split::lge
);
2622 EDomain
*Dlt
, *Deq
, *Dgt
;
2623 split_on(sp
, ind
.D
, &Dlt
, &Deq
, &Dgt
, ind
.options
);
2624 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
)
2625 assert(Dlt
|| Deq
|| Dgt
);
2626 else if (!(Dlt
|| Deq
|| Dgt
))
2627 /* Must have been empty all along */
2630 if (Deq
&& (Dlt
|| Dgt
)) {
2631 int locsize
= loc
.size();
2633 indicator
indeq(ind
, Deq
);
2635 indeq
.add_substitution(diff
);
2636 indeq
.perform_pending_substitutions();
2637 vector
<max_term
*> maxeq
= lexmin(indeq
, nparam
, loc
);
2638 maxima
.insert(maxima
.end(), maxeq
.begin(), maxeq
.end());
2639 loc
.resize(locsize
);
2642 int locsize
= loc
.size();
2644 indicator
indgt(ind
, Dgt
);
2646 /* we don't know the new location of these terms in indgt */
2648 indgt.order.lt[second].push_back(best);
2649 indgt.order.inc_pred(best);
2651 vector
<max_term
*> maxgt
= lexmin(indgt
, nparam
, loc
);
2652 maxima
.insert(maxima
.end(), maxgt
.begin(), maxgt
.end());
2653 loc
.resize(locsize
);
2658 ind
.set_domain(Deq
);
2659 ind
.add_substitution(diff
);
2660 ind
.perform_pending_substitutions();
2664 ind
.set_domain(Dlt
);
2665 ind
.order
.lt
[best
].push_back(second
);
2666 ind
.order
.inc_pred(second
);
2670 ind
.set_domain(Dgt
);
2671 ind
.order
.lt
[second
].push_back(best
);
2672 ind
.order
.inc_pred(best
);
2679 static vector
<max_term
*> lexmin(Polyhedron
*P
, Polyhedron
*C
,
2680 lexmin_options
*options
)
2682 unsigned nparam
= C
->Dimension
;
2683 Param_Polyhedron
*PP
= NULL
;
2684 Polyhedron
*CEq
= NULL
, *pVD
;
2686 Matrix
*T
= NULL
, *CP
= NULL
;
2687 Param_Domain
*D
, *next
;
2689 Polyhedron
*Porig
= P
;
2690 Polyhedron
*Corig
= C
;
2691 vector
<max_term
*> all_max
;
2693 unsigned P2PSD_MaxRays
;
2698 POL_ENSURE_VERTICES(P
);
2703 assert(P
->NbBid
== 0);
2706 remove_all_equalities(&P
, &C
, &CP
, &T
, nparam
, options
->barvinok
->MaxRays
);
2708 nparam
= CP
->NbColumns
-1;
2716 if (options
->barvinok
->MaxRays
& POL_NO_DUAL
)
2719 P2PSD_MaxRays
= options
->barvinok
->MaxRays
;
2722 PP
= Polyhedron2Param_SimplifiedDomain(&P
, C
, P2PSD_MaxRays
, &CEq
, &CT
);
2723 if (P
!= Q
&& Q
!= Porig
)
2727 if (isIdentity(CT
)) {
2731 nparam
= CT
->NbRows
- 1;
2735 unsigned dim
= P
->Dimension
- nparam
;
2738 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
2739 Polyhedron
**fVD
= new Polyhedron
*[nd
];
2741 indicator_constructor
ic(P
, dim
, PP
, T
);
2743 vector
<indicator_term
*> all_vertices
;
2744 construct_rational_vertices(PP
, T
, T
? T
->NbRows
-nparam
-1 : dim
,
2745 nparam
, all_vertices
);
2747 for (nd
= 0, D
=PP
->D
; D
; D
=next
) {
2750 Polyhedron
*rVD
= reduce_domain(D
->Domain
, CT
, CEq
,
2751 fVD
, nd
, options
->barvinok
->MaxRays
);
2755 pVD
= CT
? DomainImage(rVD
, CT
, options
->barvinok
->MaxRays
) : rVD
;
2757 EDomain
*epVD
= new EDomain(pVD
);
2758 indicator
ind(ic
, D
, epVD
, options
);
2760 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
2761 ind
.add(all_vertices
[_i
]);
2762 END_FORALL_PVertex_in_ParamPolyhedron
;
2764 ind
.remove_initial_rational_vertices();
2767 vector
<max_term
*> maxima
= lexmin(ind
, nparam
, loc
);
2769 for (int j
= 0; j
< maxima
.size(); ++j
)
2770 maxima
[j
]->substitute(CP
, options
->barvinok
);
2771 all_max
.insert(all_max
.end(), maxima
.begin(), maxima
.end());
2778 for (int i
= 0; i
< all_vertices
.size(); ++i
)
2779 delete all_vertices
[i
];
2784 Param_Polyhedron_Free(PP
);
2786 Polyhedron_Free(CEq
);
2787 for (--nd
; nd
>= 0; --nd
) {
2788 Domain_Free(fVD
[nd
]);
2799 static void verify_results(Polyhedron
*A
, Polyhedron
*C
,
2800 vector
<max_term
*>& maxima
, int m
, int M
,
2801 int print_all
, unsigned MaxRays
);
2803 int main(int argc
, char **argv
)
2808 char **iter_names
, **param_names
;
2809 int print_solution
= 1;
2810 struct lexmin_options options
;
2811 static struct argp_child argp_children
[] = {
2812 { &barvinok_argp
, 0, 0, 0 },
2813 { &verify_argp
, 0, "verification", 1 },
2816 static struct argp argp
= { argp_options
, parse_opt
, 0, 0, argp_children
};
2817 struct barvinok_options
*bv_options
;
2819 bv_options
= barvinok_options_new_with_defaults();
2820 bv_options
->lookup_table
= 0;
2822 options
.barvinok
= bv_options
;
2823 argp_parse(&argp
, argc
, argv
, 0, 0, &options
);
2826 C
= Constraints2Polyhedron(MA
, bv_options
->MaxRays
);
2828 fscanf(stdin
, " %d", &bignum
);
2829 assert(bignum
== -1);
2831 A
= Constraints2Polyhedron(MA
, bv_options
->MaxRays
);
2834 verify_options_set_range(&options
.verify
, A
);
2836 if (options
.verify
.verify
)
2839 iter_names
= util_generate_names(A
->Dimension
- C
->Dimension
, "i");
2840 param_names
= util_generate_names(C
->Dimension
, "p");
2841 if (print_solution
) {
2842 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
2843 Polyhedron_Print(stdout
, P_VALUE_FMT
, C
);
2845 vector
<max_term
*> maxima
= lexmin(A
, C
, &options
);
2847 for (int i
= 0; i
< maxima
.size(); ++i
)
2848 maxima
[i
]->print(cout
, param_names
, options
.barvinok
);
2850 if (options
.verify
.verify
)
2851 verify_results(A
, C
, maxima
, options
.verify
.m
, options
.verify
.M
,
2852 options
.verify
.print_all
, bv_options
->MaxRays
);
2854 for (int i
= 0; i
< maxima
.size(); ++i
)
2857 util_free_names(A
->Dimension
- C
->Dimension
, iter_names
);
2858 util_free_names(C
->Dimension
, param_names
);
2867 static bool lexmin(int pos
, Polyhedron
*P
, Value
*context
)
2876 value_init(LB
); value_init(UB
); value_init(k
);
2879 lu_flags
= lower_upper_bounds(pos
,P
,context
,&LB
,&UB
);
2880 assert(!(lu_flags
& LB_INFINITY
));
2882 value_set_si(context
[pos
],0);
2883 if (!lu_flags
&& value_lt(UB
,LB
)) {
2884 value_clear(LB
); value_clear(UB
); value_clear(k
);
2888 value_assign(context
[pos
], LB
);
2889 value_clear(LB
); value_clear(UB
); value_clear(k
);
2892 for (value_assign(k
,LB
); lu_flags
|| value_le(k
,UB
); value_increment(k
,k
)) {
2893 value_assign(context
[pos
],k
);
2894 if ((found
= lexmin(pos
+1, P
->next
, context
)))
2898 value_set_si(context
[pos
],0);
2899 value_clear(LB
); value_clear(UB
); value_clear(k
);
2903 static void print_list(FILE *out
, Value
*z
, char* brackets
, int len
)
2905 fprintf(out
, "%c", brackets
[0]);
2906 value_print(out
, VALUE_FMT
,z
[0]);
2907 for (int k
= 1; k
< len
; ++k
) {
2909 value_print(out
, VALUE_FMT
,z
[k
]);
2911 fprintf(out
, "%c", brackets
[1]);
2914 static int check_poly(Polyhedron
*S
, Polyhedron
*CS
, vector
<max_term
*>& maxima
,
2915 int nparam
, int pos
, Value
*z
, int print_all
, int st
,
2918 if (pos
== nparam
) {
2920 bool found
= lexmin(1, S
, z
);
2924 print_list(stdout
, z
+S
->Dimension
-nparam
+1, "()", nparam
);
2927 print_list(stdout
, z
+1, "[]", S
->Dimension
-nparam
);
2932 for (int i
= 0; i
< maxima
.size(); ++i
)
2933 if ((min
= maxima
[i
]->eval(z
+S
->Dimension
-nparam
+1, MaxRays
)))
2936 int ok
= !(found
^ !!min
);
2938 for (int i
= 0; i
< S
->Dimension
-nparam
; ++i
)
2939 if (value_ne(z
[1+i
], min
->p
[i
])) {
2946 fprintf(stderr
, "Error !\n");
2947 fprintf(stderr
, "lexmin");
2948 print_list(stderr
, z
+S
->Dimension
-nparam
+1, "()", nparam
);
2949 fprintf(stderr
, " should be ");
2951 print_list(stderr
, z
+1, "[]", S
->Dimension
-nparam
);
2952 fprintf(stderr
, " while digging gives ");
2954 print_list(stderr
, min
->p
, "[]", S
->Dimension
-nparam
);
2955 fprintf(stderr
, ".\n");
2957 } else if (print_all
)
2962 for (k
= 1; k
<= S
->Dimension
-nparam
; ++k
)
2963 value_set_si(z
[k
], 0);
2971 !(lower_upper_bounds(1+pos
, CS
, &z
[S
->Dimension
-nparam
], &LB
, &UB
));
2972 for (value_assign(tmp
,LB
); value_le(tmp
,UB
); value_increment(tmp
,tmp
)) {
2974 int k
= VALUE_TO_INT(tmp
);
2975 if (!pos
&& !(k
%st
)) {
2980 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
2981 if (!check_poly(S
, CS
->next
, maxima
, nparam
, pos
+1, z
, print_all
, st
,
2989 value_set_si(z
[pos
+S
->Dimension
-nparam
+1],0);
2997 void verify_results(Polyhedron
*A
, Polyhedron
*C
, vector
<max_term
*>& maxima
,
2998 int m
, int M
, int print_all
, unsigned MaxRays
)
3000 Polyhedron
*CC
, *CC2
, *CS
, *S
;
3001 unsigned nparam
= C
->Dimension
;
3006 CC
= Polyhedron_Project(A
, nparam
);
3007 CC2
= DomainIntersection(C
, CC
, MaxRays
);
3011 /* Intersect context with range */
3016 MM
= Matrix_Alloc(2*C
->Dimension
, C
->Dimension
+2);
3017 for (int i
= 0; i
< C
->Dimension
; ++i
) {
3018 value_set_si(MM
->p
[2*i
][0], 1);
3019 value_set_si(MM
->p
[2*i
][1+i
], 1);
3020 value_set_si(MM
->p
[2*i
][1+C
->Dimension
], -m
);
3021 value_set_si(MM
->p
[2*i
+1][0], 1);
3022 value_set_si(MM
->p
[2*i
+1][1+i
], -1);
3023 value_set_si(MM
->p
[2*i
+1][1+C
->Dimension
], M
);
3025 CC2
= AddConstraints(MM
->p
[0], 2*CC
->Dimension
, CC
, MaxRays
);
3026 U
= Universe_Polyhedron(0);
3027 CS
= Polyhedron_Scan(CC2
, U
, MaxRays
& POL_NO_DUAL
? 0 : MaxRays
);
3029 Polyhedron_Free(CC2
);
3034 p
= ALLOCN(Value
, A
->Dimension
+2);
3035 for (i
=0; i
<= A
->Dimension
; i
++) {
3037 value_set_si(p
[i
],0);
3040 value_set_si(p
[i
], 1);
3042 S
= Polyhedron_Scan(A
, C
, MaxRays
& POL_NO_DUAL
? 0 : MaxRays
);
3044 if (!print_all
&& C
->Dimension
> 0) {
3049 for (int i
= m
; i
<= M
; i
+= st
)
3056 if (!(CS
&& emptyQ2(CS
)))
3057 check_poly(S
, CS
, maxima
, nparam
, 0, p
, print_all
, st
, MaxRays
);
3064 for (i
=0; i
<= (A
->Dimension
+1); i
++)
3069 Polyhedron_Free(CC
);