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>
14 #include "conversion.h"
15 #include "decomposer.h"
16 #include "lattice_point.h"
17 #include "reduce_domain.h"
21 #include "evalue_util.h"
22 #include "remove_equalities.h"
27 #undef CS /* for Solaris 10 */
40 #define EMPTINESS_CHECK (BV_OPT_LAST+1)
41 #define NO_REDUCTION (BV_OPT_LAST+2)
43 struct argp_option argp_options
[] = {
44 { "emptiness-check", EMPTINESS_CHECK
, "[none|count]", 0 },
45 { "no-reduction", NO_REDUCTION
, 0, 0 },
49 static error_t
parse_opt(int key
, char *arg
, struct argp_state
*state
)
51 struct lexmin_options
*options
= (struct lexmin_options
*)(state
->input
);
52 struct barvinok_options
*bv_options
= options
->verify
.barvinok
;
56 state
->child_inputs
[0] = options
->verify
.barvinok
;
57 state
->child_inputs
[1] = &options
->verify
;
58 options
->emptiness_check
= BV_LEXMIN_EMPTINESS_CHECK_SAMPLE
;
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 return ARGP_ERR_UNKNOWN
;
78 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
80 static int type_offset(enode
*p
)
82 return p
->type
== fractional
? 1 :
83 p
->type
== flooring
? 1 : 0;
86 void compute_evalue(evalue
*e
, Value
*val
, Value
*res
)
88 double d
= compute_evalue(e
, val
);
93 value_set_double(*res
, d
);
96 struct indicator_term
{
98 int pos
; /* number of rational vertex */
99 int n
; /* number of cone associated to given rational vertex */
103 indicator_term(unsigned dim
, int pos
) {
105 vertex
= new evalue
* [dim
];
110 indicator_term(unsigned dim
, int pos
, int n
) {
111 den
.SetDims(dim
, dim
);
112 vertex
= new evalue
* [dim
];
116 indicator_term(const indicator_term
& src
) {
121 unsigned dim
= den
.NumCols();
122 vertex
= new evalue
* [dim
];
123 for (int i
= 0; i
< dim
; ++i
) {
124 vertex
[i
] = new evalue();
125 value_init(vertex
[i
]->d
);
126 evalue_copy(vertex
[i
], src
.vertex
[i
]);
129 void swap(indicator_term
*other
) {
131 tmp
= sign
; sign
= other
->sign
; other
->sign
= tmp
;
132 tmp
= pos
; pos
= other
->pos
; other
->pos
= tmp
;
133 tmp
= n
; n
= other
->n
; other
->n
= tmp
;
134 mat_ZZ tmp_den
= den
; den
= other
->den
; other
->den
= tmp_den
;
135 unsigned dim
= den
.NumCols();
136 for (int i
= 0; i
< dim
; ++i
) {
137 evalue
*tmp
= vertex
[i
];
138 vertex
[i
] = other
->vertex
[i
];
139 other
->vertex
[i
] = tmp
;
143 unsigned dim
= den
.NumCols();
144 for (int i
= 0; i
< dim
; ++i
) {
145 free_evalue_refs(vertex
[i
]);
150 void print(ostream
& os
, char **p
) const;
151 void substitute(Matrix
*T
);
153 void substitute(evalue
*fract
, evalue
*val
);
154 void substitute(int pos
, evalue
*val
);
155 void reduce_in_domain(Polyhedron
*D
);
156 bool is_opposite(const indicator_term
*neg
) const;
157 vec_ZZ
eval(Value
*val
) const {
159 unsigned dim
= den
.NumCols();
163 for (int i
= 0; i
< dim
; ++i
) {
164 compute_evalue(vertex
[i
], val
, &tmp
);
172 static int evalue_rational_cmp(const evalue
*e1
, const evalue
*e2
)
180 assert(value_notzero_p(e1
->d
));
181 assert(value_notzero_p(e2
->d
));
182 value_multiply(m
, e1
->x
.n
, e2
->d
);
183 value_multiply(m2
, e2
->x
.n
, e1
->d
);
186 else if (value_gt(m
, m2
))
196 static int evalue_cmp(const evalue
*e1
, const evalue
*e2
)
198 if (value_notzero_p(e1
->d
)) {
199 if (value_zero_p(e2
->d
))
201 return evalue_rational_cmp(e1
, e2
);
203 if (value_notzero_p(e2
->d
))
205 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
206 return e1
->x
.p
->type
- e2
->x
.p
->type
;
207 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
208 return e1
->x
.p
->size
- e2
->x
.p
->size
;
209 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
210 return e1
->x
.p
->pos
- e2
->x
.p
->pos
;
211 assert(e1
->x
.p
->type
== polynomial
||
212 e1
->x
.p
->type
== fractional
||
213 e1
->x
.p
->type
== flooring
);
214 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
) {
215 int s
= evalue_cmp(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]);
222 void evalue_length(evalue
*e
, int len
[2])
227 while (value_zero_p(e
->d
)) {
228 assert(e
->x
.p
->type
== polynomial
||
229 e
->x
.p
->type
== fractional
||
230 e
->x
.p
->type
== flooring
);
231 if (e
->x
.p
->type
== polynomial
)
235 int offset
= type_offset(e
->x
.p
);
236 assert(e
->x
.p
->size
== offset
+2);
237 e
= &e
->x
.p
->arr
[offset
];
241 static bool it_smaller(const indicator_term
* it1
, const indicator_term
* it2
)
245 int len1
[2], len2
[2];
246 unsigned dim
= it1
->den
.NumCols();
247 for (int i
= 0; i
< dim
; ++i
) {
248 evalue_length(it1
->vertex
[i
], len1
);
249 evalue_length(it2
->vertex
[i
], len2
);
250 if (len1
[0] != len2
[0])
251 return len1
[0] < len2
[0];
252 if (len1
[1] != len2
[1])
253 return len1
[1] < len2
[1];
255 if (it1
->pos
!= it2
->pos
)
256 return it1
->pos
< it2
->pos
;
257 if (it1
->n
!= it2
->n
)
258 return it1
->n
< it2
->n
;
259 int s
= lex_cmp(it1
->den
, it2
->den
);
262 for (int i
= 0; i
< dim
; ++i
) {
263 s
= evalue_cmp(it1
->vertex
[i
], it2
->vertex
[i
]);
267 assert(it1
->sign
!= 0);
268 assert(it2
->sign
!= 0);
269 if (it1
->sign
!= it2
->sign
)
270 return it1
->sign
> 0;
275 static const int requires_resort
;
276 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
277 return it_smaller(it1
, it2
);
280 const int smaller_it::requires_resort
= 1;
282 struct smaller_it_p
{
283 static const int requires_resort
;
284 bool operator()(const indicator_term
* it1
, const indicator_term
* it2
) const {
288 const int smaller_it_p::requires_resort
= 0;
290 /* Returns true if this and neg are opposite using the knowledge
291 * that they have the same numerator.
292 * In particular, we check that the signs are different and that
293 * the denominator is the same.
295 bool indicator_term::is_opposite(const indicator_term
*neg
) const
297 if (sign
+ neg
->sign
!= 0)
304 void indicator_term::reduce_in_domain(Polyhedron
*D
)
306 for (int k
= 0; k
< den
.NumCols(); ++k
) {
307 reduce_evalue_in_domain(vertex
[k
], D
);
308 if (evalue_range_reduction_in_domain(vertex
[k
], D
))
309 reduce_evalue(vertex
[k
]);
313 void indicator_term::print(ostream
& os
, char **p
) const
315 unsigned dim
= den
.NumCols();
316 unsigned factors
= den
.NumRows();
324 for (int i
= 0; i
< dim
; ++i
) {
327 evalue_print(os
, vertex
[i
], p
);
330 for (int i
= 0; i
< factors
; ++i
) {
331 os
<< " + t" << i
<< "*[";
332 for (int j
= 0; j
< dim
; ++j
) {
339 os
<< " ((" << pos
<< ", " << n
<< ", " << (void*)this << "))";
342 /* Perform the substitution specified by T on the variables.
343 * T has dimension (newdim+nparam+1) x (olddim + nparam + 1).
344 * The substitution is performed as in gen_fun::substitute
346 void indicator_term::substitute(Matrix
*T
)
348 unsigned dim
= den
.NumCols();
349 unsigned nparam
= T
->NbColumns
- dim
- 1;
350 unsigned newdim
= T
->NbRows
- nparam
- 1;
353 matrix2zz(T
, trans
, newdim
, dim
);
354 trans
= transpose(trans
);
356 newvertex
= new evalue
* [newdim
];
359 v
.SetLength(nparam
+1);
362 value_init(factor
.d
);
363 value_set_si(factor
.d
, 1);
364 value_init(factor
.x
.n
);
365 for (int i
= 0; i
< newdim
; ++i
) {
366 values2zz(T
->p
[i
]+dim
, v
, nparam
+1);
367 newvertex
[i
] = multi_monom(v
);
369 for (int j
= 0; j
< dim
; ++j
) {
370 if (value_zero_p(T
->p
[i
][j
]))
374 evalue_copy(&term
, vertex
[j
]);
375 value_assign(factor
.x
.n
, T
->p
[i
][j
]);
376 emul(&factor
, &term
);
377 eadd(&term
, newvertex
[i
]);
378 free_evalue_refs(&term
);
381 free_evalue_refs(&factor
);
382 for (int i
= 0; i
< dim
; ++i
) {
383 free_evalue_refs(vertex
[i
]);
390 static void evalue_add_constant(evalue
*e
, ZZ v
)
395 /* go down to constant term */
396 while (value_zero_p(e
->d
))
397 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)];
400 value_multiply(tmp
, tmp
, e
->d
);
401 value_addto(e
->x
.n
, e
->x
.n
, tmp
);
406 /* Make all powers in denominator lexico-positive */
407 void indicator_term::normalize()
410 extra_vertex
.SetLength(den
.NumCols());
411 for (int r
= 0; r
< den
.NumRows(); ++r
) {
412 for (int k
= 0; k
< den
.NumCols(); ++k
) {
419 extra_vertex
+= den
[r
];
423 for (int k
= 0; k
< extra_vertex
.length(); ++k
)
424 if (extra_vertex
[k
] != 0)
425 evalue_add_constant(vertex
[k
], extra_vertex
[k
]);
428 static void substitute(evalue
*e
, evalue
*fract
, evalue
*val
)
432 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
433 if (t
->x
.p
->type
== fractional
&& eequal(&t
->x
.p
->arr
[0], fract
))
436 if (value_notzero_p(t
->d
))
439 free_evalue_refs(&t
->x
.p
->arr
[0]);
440 evalue
*term
= &t
->x
.p
->arr
[2];
447 free_evalue_refs(term
);
453 void indicator_term::substitute(evalue
*fract
, evalue
*val
)
455 unsigned dim
= den
.NumCols();
456 for (int i
= 0; i
< dim
; ++i
) {
457 ::substitute(vertex
[i
], fract
, val
);
461 static void substitute(evalue
*e
, int pos
, evalue
*val
)
465 for (t
= e
; value_zero_p(t
->d
); t
= &t
->x
.p
->arr
[type_offset(t
->x
.p
)]) {
466 if (t
->x
.p
->type
== polynomial
&& t
->x
.p
->pos
== pos
)
469 if (value_notzero_p(t
->d
))
472 evalue
*term
= &t
->x
.p
->arr
[1];
479 free_evalue_refs(term
);
485 void indicator_term::substitute(int pos
, evalue
*val
)
487 unsigned dim
= den
.NumCols();
488 for (int i
= 0; i
< dim
; ++i
) {
489 ::substitute(vertex
[i
], pos
, val
);
493 struct indicator_constructor
: public signed_cone_consumer
,
494 public vertex_decomposer
{
496 vector
<indicator_term
*> *terms
;
497 Matrix
*T
; /* Transformation to original space */
498 Param_Polyhedron
*PP
;
502 indicator_constructor(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
,
504 vertex_decomposer(P
, PP
->nbV
, *this), T(T
), PP(PP
) {
505 vertex
.SetLength(dim
);
506 terms
= new vector
<indicator_term
*>[nbV
];
508 ~indicator_constructor() {
509 for (int i
= 0; i
< nbV
; ++i
)
510 for (int j
= 0; j
< terms
[i
].size(); ++j
)
514 void substitute(Matrix
*T
);
516 void print(ostream
& os
, char **p
);
518 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
519 void decompose_at_vertex(Param_Vertices
*V
, int _i
,
520 barvinok_options
*options
) {
523 vertex_decomposer::decompose_at_vertex(V
, _i
, options
);
527 void indicator_constructor::handle(const signed_cone
& sc
, barvinok_options
*options
)
531 unsigned dim
= vertex
.length();
533 assert(sc
.rays
.NumRows() == dim
);
535 indicator_term
*term
= new indicator_term(dim
, pos
, n
++);
536 term
->sign
= sc
.sign
;
537 terms
[vert
].push_back(term
);
539 lattice_point(V
, sc
.rays
, vertex
, term
->vertex
, options
);
542 for (int r
= 0; r
< dim
; ++r
) {
543 for (int j
= 0; j
< dim
; ++j
) {
544 if (term
->den
[r
][j
] == 0)
546 if (term
->den
[r
][j
] > 0)
548 term
->sign
= -term
->sign
;
549 term
->den
[r
] = -term
->den
[r
];
550 vertex
+= term
->den
[r
];
555 for (int i
= 0; i
< dim
; ++i
) {
556 if (!term
->vertex
[i
]) {
557 term
->vertex
[i
] = new evalue();
558 value_init(term
->vertex
[i
]->d
);
559 value_init(term
->vertex
[i
]->x
.n
);
560 zz2value(vertex
[i
], term
->vertex
[i
]->x
.n
);
561 value_set_si(term
->vertex
[i
]->d
, 1);
566 evalue_add_constant(term
->vertex
[i
], vertex
[i
]);
574 lex_order_rows(term
->den
);
577 void indicator_constructor::print(ostream
& os
, char **p
)
579 for (int i
= 0; i
< nbV
; ++i
)
580 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
581 os
<< "i: " << i
<< ", j: " << j
<< endl
;
582 terms
[i
][j
]->print(os
, p
);
587 void indicator_constructor::normalize()
589 for (int i
= 0; i
< nbV
; ++i
)
590 for (int j
= 0; j
< terms
[i
].size(); ++j
) {
592 vertex
.SetLength(terms
[i
][j
]->den
.NumCols());
593 for (int r
= 0; r
< terms
[i
][j
]->den
.NumRows(); ++r
) {
594 for (int k
= 0; k
< terms
[i
][j
]->den
.NumCols(); ++k
) {
595 if (terms
[i
][j
]->den
[r
][k
] == 0)
597 if (terms
[i
][j
]->den
[r
][k
] > 0)
599 terms
[i
][j
]->sign
= -terms
[i
][j
]->sign
;
600 terms
[i
][j
]->den
[r
] = -terms
[i
][j
]->den
[r
];
601 vertex
+= terms
[i
][j
]->den
[r
];
605 lex_order_rows(terms
[i
][j
]->den
);
606 for (int k
= 0; k
< vertex
.length(); ++k
)
608 evalue_add_constant(terms
[i
][j
]->vertex
[k
], vertex
[k
]);
612 struct order_cache_el
{
614 order_cache_el
copy() const {
616 for (int i
= 0; i
< e
.size(); ++i
) {
617 evalue
*c
= new evalue
;
619 evalue_copy(c
, e
[i
]);
625 for (int i
= 0; i
< e
.size(); ++i
) {
626 free_evalue_refs(e
[i
]);
633 evalue_set_si(&mone
, -1, 1);
634 for (int i
= 0; i
< e
.size(); ++i
)
636 free_evalue_refs(&mone
);
638 void print(ostream
& os
, char **p
);
641 void order_cache_el::print(ostream
& os
, char **p
)
644 for (int i
= 0; i
< e
.size(); ++i
) {
647 evalue_print(os
, e
[i
], p
);
653 vector
<order_cache_el
> lt
;
654 vector
<order_cache_el
> le
;
655 vector
<order_cache_el
> unknown
;
657 void clear_transients() {
658 for (int i
= 0; i
< le
.size(); ++i
)
660 for (int i
= 0; i
< unknown
.size(); ++i
)
667 for (int i
= 0; i
< lt
.size(); ++i
)
671 void add(order_cache_el
& cache_el
, order_sign sign
);
672 order_sign
check_lt(vector
<order_cache_el
>* list
,
673 const indicator_term
*a
, const indicator_term
*b
,
674 order_cache_el
& cache_el
);
675 order_sign
check_lt(const indicator_term
*a
, const indicator_term
*b
,
676 order_cache_el
& cache_el
);
677 order_sign
check_direct(const indicator_term
*a
, const indicator_term
*b
,
678 order_cache_el
& cache_el
);
679 order_sign
check(const indicator_term
*a
, const indicator_term
*b
,
680 order_cache_el
& cache_el
);
681 void copy(const order_cache
& cache
);
682 void print(ostream
& os
, char **p
);
685 void order_cache::copy(const order_cache
& cache
)
687 for (int i
= 0; i
< cache
.lt
.size(); ++i
) {
688 order_cache_el n
= cache
.lt
[i
].copy();
693 void order_cache::add(order_cache_el
& cache_el
, order_sign sign
)
695 if (sign
== order_lt
) {
696 lt
.push_back(cache_el
);
697 } else if (sign
== order_gt
) {
699 lt
.push_back(cache_el
);
700 } else if (sign
== order_le
) {
701 le
.push_back(cache_el
);
702 } else if (sign
== order_ge
) {
704 le
.push_back(cache_el
);
705 } else if (sign
== order_unknown
) {
706 unknown
.push_back(cache_el
);
708 assert(sign
== order_eq
);
715 static evalue
*ediff(const evalue
*a
, const evalue
*b
)
719 evalue_set_si(&mone
, -1, 1);
720 evalue
*diff
= new evalue
;
722 evalue_copy(diff
, b
);
726 free_evalue_refs(&mone
);
730 static bool evalue_first_difference(const evalue
*e1
, const evalue
*e2
,
731 const evalue
**d1
, const evalue
**d2
)
736 if (value_ne(e1
->d
, e2
->d
))
739 if (value_notzero_p(e1
->d
)) {
740 if (value_eq(e1
->x
.n
, e2
->x
.n
))
744 if (e1
->x
.p
->type
!= e2
->x
.p
->type
)
746 if (e1
->x
.p
->size
!= e2
->x
.p
->size
)
748 if (e1
->x
.p
->pos
!= e2
->x
.p
->pos
)
751 assert(e1
->x
.p
->type
== polynomial
||
752 e1
->x
.p
->type
== fractional
||
753 e1
->x
.p
->type
== flooring
);
754 int offset
= type_offset(e1
->x
.p
);
755 assert(e1
->x
.p
->size
== offset
+2);
756 for (int i
= 0; i
< e1
->x
.p
->size
; ++i
)
757 if (i
!= type_offset(e1
->x
.p
) &&
758 !eequal(&e1
->x
.p
->arr
[i
], &e2
->x
.p
->arr
[i
]))
761 return evalue_first_difference(&e1
->x
.p
->arr
[offset
],
762 &e2
->x
.p
->arr
[offset
], d1
, d2
);
765 static order_sign
evalue_diff_constant_sign(const evalue
*e1
, const evalue
*e2
)
767 if (!evalue_first_difference(e1
, e2
, &e1
, &e2
))
769 if (value_zero_p(e1
->d
) || value_zero_p(e2
->d
))
770 return order_undefined
;
771 int s
= evalue_rational_cmp(e1
, e2
);
780 order_sign
order_cache::check_lt(vector
<order_cache_el
>* list
,
781 const indicator_term
*a
, const indicator_term
*b
,
782 order_cache_el
& cache_el
)
784 order_sign sign
= order_undefined
;
785 for (int i
= 0; i
< list
->size(); ++i
) {
787 for (j
= cache_el
.e
.size(); j
< (*list
)[i
].e
.size(); ++j
)
788 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
789 for (j
= 0; j
< (*list
)[i
].e
.size(); ++j
) {
790 order_sign diff_sign
;
791 diff_sign
= evalue_diff_constant_sign((*list
)[i
].e
[j
], cache_el
.e
[j
]);
792 if (diff_sign
== order_gt
) {
795 } else if (diff_sign
== order_lt
)
797 else if (diff_sign
== order_undefined
)
800 assert(diff_sign
== order_eq
);
802 if (j
== (*list
)[i
].e
.size())
803 sign
= list
== <
? order_lt
: order_le
;
804 if (sign
!= order_undefined
)
810 order_sign
order_cache::check_direct(const indicator_term
*a
,
811 const indicator_term
*b
,
812 order_cache_el
& cache_el
)
814 order_sign sign
= check_lt(<
, a
, b
, cache_el
);
815 if (sign
!= order_undefined
)
817 sign
= check_lt(&le
, a
, b
, cache_el
);
818 if (sign
!= order_undefined
)
821 for (int i
= 0; i
< unknown
.size(); ++i
) {
823 for (j
= cache_el
.e
.size(); j
< unknown
[i
].e
.size(); ++j
)
824 cache_el
.e
.push_back(ediff(a
->vertex
[j
], b
->vertex
[j
]));
825 for (j
= 0; j
< unknown
[i
].e
.size(); ++j
) {
826 if (!eequal(unknown
[i
].e
[j
], cache_el
.e
[j
]))
829 if (j
== unknown
[i
].e
.size()) {
830 sign
= order_unknown
;
837 order_sign
order_cache::check(const indicator_term
*a
, const indicator_term
*b
,
838 order_cache_el
& cache_el
)
840 order_sign sign
= check_direct(a
, b
, cache_el
);
841 if (sign
!= order_undefined
)
843 int size
= cache_el
.e
.size();
845 sign
= check_direct(a
, b
, cache_el
);
847 assert(cache_el
.e
.size() == size
);
848 if (sign
== order_undefined
)
850 if (sign
== order_lt
)
852 else if (sign
== order_le
)
855 assert(sign
== order_unknown
);
861 struct partial_order
{
864 std::set
<const indicator_term
*, smaller_it
> head
;
865 map
<const indicator_term
*, int, smaller_it
> pred
;
866 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> lt
;
867 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> le
;
868 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> eq
;
870 map
<const indicator_term
*, vector
<const indicator_term
* >, smaller_it
> pending
;
874 partial_order(indicator
*ind
) : ind(ind
) {}
875 void copy(const partial_order
& order
,
876 map
< const indicator_term
*, indicator_term
* > old2new
);
878 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
879 map
<const indicator_term
*, int >::iterator j
;
880 std::set
<const indicator_term
*>::iterator k
;
882 if (head
.key_comp().requires_resort
) {
883 typeof(head
) new_head
;
884 for (k
= head
.begin(); k
!= head
.end(); ++k
)
890 if (pred
.key_comp().requires_resort
) {
891 typeof(pred
) new_pred
;
892 for (j
= pred
.begin(); j
!= pred
.end(); ++j
)
893 new_pred
[(*j
).first
] = (*j
).second
;
898 if (lt
.key_comp().requires_resort
) {
900 for (i
= lt
.begin(); i
!= lt
.end(); ++i
)
901 m
[(*i
).first
] = (*i
).second
;
906 if (le
.key_comp().requires_resort
) {
908 for (i
= le
.begin(); i
!= le
.end(); ++i
)
909 m
[(*i
).first
] = (*i
).second
;
914 if (eq
.key_comp().requires_resort
) {
916 for (i
= eq
.begin(); i
!= eq
.end(); ++i
)
917 m
[(*i
).first
] = (*i
).second
;
922 if (pending
.key_comp().requires_resort
) {
924 for (i
= pending
.begin(); i
!= pending
.end(); ++i
)
925 m
[(*i
).first
] = (*i
).second
;
931 order_sign
compare(const indicator_term
*a
, const indicator_term
*b
);
932 void set_equal(const indicator_term
*a
, const indicator_term
*b
);
933 void unset_le(const indicator_term
*a
, const indicator_term
*b
);
934 void dec_pred(const indicator_term
*it
) {
935 if (--pred
[it
] == 0) {
940 void inc_pred(const indicator_term
*it
) {
941 if (head
.find(it
) != head
.end())
946 bool compared(const indicator_term
* a
, const indicator_term
* b
);
947 void add(const indicator_term
* it
, std::set
<const indicator_term
*> *filter
);
948 void remove(const indicator_term
* it
);
950 void print(ostream
& os
, char **p
);
952 /* replace references to orig to references to replacement */
953 void replace(const indicator_term
* orig
, indicator_term
* replacement
);
954 void sanity_check() const;
957 /* We actually replace the contents of orig by that of replacement,
958 * but we have to be careful since replacing the content changes
959 * the order of orig in the maps.
961 void partial_order::replace(const indicator_term
* orig
, indicator_term
* replacement
)
963 std::set
<const indicator_term
*>::iterator k
;
965 bool is_head
= k
!= head
.end();
970 orig_pred
= pred
[orig
];
973 vector
<const indicator_term
* > orig_lt
;
974 vector
<const indicator_term
* > orig_le
;
975 vector
<const indicator_term
* > orig_eq
;
976 vector
<const indicator_term
* > orig_pending
;
977 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
978 bool in_lt
= ((i
= lt
.find(orig
)) != lt
.end());
980 orig_lt
= (*i
).second
;
983 bool in_le
= ((i
= le
.find(orig
)) != le
.end());
985 orig_le
= (*i
).second
;
988 bool in_eq
= ((i
= eq
.find(orig
)) != eq
.end());
990 orig_eq
= (*i
).second
;
993 bool in_pending
= ((i
= pending
.find(orig
)) != pending
.end());
995 orig_pending
= (*i
).second
;
998 indicator_term
*old
= const_cast<indicator_term
*>(orig
);
999 old
->swap(replacement
);
1003 pred
[old
] = orig_pred
;
1011 pending
[old
] = orig_pending
;
1014 void partial_order::unset_le(const indicator_term
*a
, const indicator_term
*b
)
1016 vector
<const indicator_term
*>::iterator i
;
1017 i
= find(le
[a
].begin(), le
[a
].end(), b
);
1019 if (le
[a
].size() == 0)
1022 i
= find(pending
[a
].begin(), pending
[a
].end(), b
);
1023 if (i
!= pending
[a
].end())
1024 pending
[a
].erase(i
);
1027 void partial_order::set_equal(const indicator_term
*a
, const indicator_term
*b
)
1029 if (eq
[a
].size() == 0)
1031 if (eq
[b
].size() == 0)
1036 if (pred
.key_comp()(b
, a
)) {
1037 const indicator_term
*c
= a
;
1042 const indicator_term
*base
= a
;
1044 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1046 for (int j
= 0; j
< eq
[b
].size(); ++j
) {
1047 eq
[base
].push_back(eq
[b
][j
]);
1048 eq
[eq
[b
][j
]][0] = base
;
1053 if (i
!= lt
.end()) {
1054 for (int j
= 0; j
< lt
[b
].size(); ++j
) {
1055 if (find(eq
[base
].begin(), eq
[base
].end(), lt
[b
][j
]) != eq
[base
].end())
1057 else if (find(lt
[base
].begin(), lt
[base
].end(), lt
[b
][j
])
1061 lt
[base
].push_back(lt
[b
][j
]);
1067 if (i
!= le
.end()) {
1068 for (int j
= 0; j
< le
[b
].size(); ++j
) {
1069 if (find(eq
[base
].begin(), eq
[base
].end(), le
[b
][j
]) != eq
[base
].end())
1071 else if (find(le
[base
].begin(), le
[base
].end(), le
[b
][j
])
1075 le
[base
].push_back(le
[b
][j
]);
1080 i
= pending
.find(base
);
1081 if (i
!= pending
.end()) {
1082 vector
<const indicator_term
* > old
= pending
[base
];
1083 pending
[base
].clear();
1084 for (int j
= 0; j
< old
.size(); ++j
) {
1085 if (find(eq
[base
].begin(), eq
[base
].end(), old
[j
]) == eq
[base
].end())
1086 pending
[base
].push_back(old
[j
]);
1090 i
= pending
.find(b
);
1091 if (i
!= pending
.end()) {
1092 for (int j
= 0; j
< pending
[b
].size(); ++j
) {
1093 if (find(eq
[base
].begin(), eq
[base
].end(), pending
[b
][j
]) == eq
[base
].end())
1094 pending
[base
].push_back(pending
[b
][j
]);
1100 void partial_order::copy(const partial_order
& order
,
1101 map
< const indicator_term
*, indicator_term
* > old2new
)
1103 cache
.copy(order
.cache
);
1105 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator i
;
1106 map
<const indicator_term
*, int >::const_iterator j
;
1107 std::set
<const indicator_term
*>::const_iterator k
;
1109 for (k
= order
.head
.begin(); k
!= order
.head
.end(); ++k
)
1110 head
.insert(old2new
[*k
]);
1112 for (j
= order
.pred
.begin(); j
!= order
.pred
.end(); ++j
)
1113 pred
[old2new
[(*j
).first
]] = (*j
).second
;
1115 for (i
= order
.lt
.begin(); i
!= order
.lt
.end(); ++i
) {
1116 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1117 lt
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1119 for (i
= order
.le
.begin(); i
!= order
.le
.end(); ++i
) {
1120 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1121 le
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1123 for (i
= order
.eq
.begin(); i
!= order
.eq
.end(); ++i
) {
1124 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1125 eq
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1127 for (i
= order
.pending
.begin(); i
!= order
.pending
.end(); ++i
) {
1128 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1129 pending
[old2new
[(*i
).first
]].push_back(old2new
[(*i
).second
[j
]]);
1135 vector
<evalue
*> max
;
1137 void print(ostream
& os
, char **p
, barvinok_options
*options
) const;
1138 void substitute(Matrix
*T
, barvinok_options
*options
);
1139 Vector
*eval(Value
*val
, unsigned MaxRays
) const;
1142 for (int i
= 0; i
< max
.size(); ++i
) {
1143 free_evalue_refs(max
[i
]);
1151 * Project on first dim dimensions
1153 Polyhedron
* Polyhedron_Project_Initial(Polyhedron
*P
, int dim
)
1159 if (P
->Dimension
== dim
)
1160 return Polyhedron_Copy(P
);
1162 T
= Matrix_Alloc(dim
+1, P
->Dimension
+1);
1163 for (i
= 0; i
< dim
; ++i
)
1164 value_set_si(T
->p
[i
][i
], 1);
1165 value_set_si(T
->p
[dim
][P
->Dimension
], 1);
1166 I
= Polyhedron_Image(P
, T
, P
->NbConstraints
);
1172 vector
<indicator_term
*> term
;
1173 indicator_constructor
& ic
;
1174 partial_order order
;
1178 lexmin_options
*options
;
1179 vector
<evalue
*> substitutions
;
1181 indicator(indicator_constructor
& ic
, Param_Domain
*PD
, EDomain
*D
,
1182 lexmin_options
*options
) :
1183 ic(ic
), PD(PD
), D(D
), order(this), options(options
), P(NULL
) {}
1184 indicator(const indicator
& ind
, EDomain
*D
) :
1185 ic(ind
.ic
), PD(ind
.PD
), D(NULL
), order(this), options(ind
.options
),
1186 P(Polyhedron_Copy(ind
.P
)) {
1187 map
< const indicator_term
*, indicator_term
* > old2new
;
1188 for (int i
= 0; i
< ind
.term
.size(); ++i
) {
1189 indicator_term
*it
= new indicator_term(*ind
.term
[i
]);
1190 old2new
[ind
.term
[i
]] = it
;
1193 order
.copy(ind
.order
, old2new
);
1197 for (int i
= 0; i
< term
.size(); ++i
)
1205 void set_domain(EDomain
*D
) {
1206 order
.cache
.clear_transients();
1210 int nparam
= ic
.P
->Dimension
- ic
.vertex
.length();
1211 if (options
->reduce
) {
1212 Polyhedron
*Q
= Polyhedron_Project_Initial(D
->D
, nparam
);
1213 Q
= DomainConstraintSimplify(Q
, options
->verify
.barvinok
->MaxRays
);
1214 if (!P
|| !PolyhedronIncludes(Q
, P
))
1215 reduce_in_domain(Q
);
1223 void add(const indicator_term
* it
);
1224 void remove(const indicator_term
* it
);
1225 void remove_initial_rational_vertices();
1226 void expand_rational_vertex(const indicator_term
*initial
);
1228 void print(ostream
& os
, char **p
);
1230 void peel(int i
, int j
);
1231 void combine(const indicator_term
*a
, const indicator_term
*b
);
1232 void add_substitution(evalue
*equation
);
1233 void perform_pending_substitutions();
1234 void reduce_in_domain(Polyhedron
*D
);
1235 bool handle_equal_numerators(const indicator_term
*base
);
1237 max_term
* create_max_term(const indicator_term
*it
);
1239 void substitute(evalue
*equation
);
1242 void partial_order::sanity_check() const
1244 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator i
;
1245 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator prev
;
1246 map
<const indicator_term
*, vector
<const indicator_term
* > >::const_iterator l
;
1247 map
<const indicator_term
*, int >::const_iterator k
, prev_k
;
1249 for (k
= pred
.begin(); k
!= pred
.end(); prev_k
= k
, ++k
)
1250 if (k
!= pred
.begin())
1251 assert(pred
.key_comp()((*prev_k
).first
, (*k
).first
));
1252 for (i
= lt
.begin(); i
!= lt
.end(); prev
= i
, ++i
) {
1255 i_v
= (*i
).first
->eval(ind
->D
->sample
->p
);
1256 if (i
!= lt
.begin())
1257 assert(lt
.key_comp()((*prev
).first
, (*i
).first
));
1258 l
= eq
.find((*i
).first
);
1260 assert((*l
).second
.size() > 1);
1261 assert(head
.find((*i
).first
) != head
.end() ||
1262 pred
.find((*i
).first
) != pred
.end());
1263 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1264 k
= pred
.find((*i
).second
[j
]);
1265 assert(k
!= pred
.end());
1266 assert((*k
).second
!= 0);
1267 if ((*i
).first
->sign
!= 0 &&
1268 (*i
).second
[j
]->sign
!= 0 && ind
->D
->sample
) {
1269 vec_ZZ j_v
= (*i
).second
[j
]->eval(ind
->D
->sample
->p
);
1270 assert(lex_cmp(i_v
, j_v
) < 0);
1274 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1275 assert((*i
).second
.size() > 0);
1276 assert(head
.find((*i
).first
) != head
.end() ||
1277 pred
.find((*i
).first
) != pred
.end());
1278 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1279 k
= pred
.find((*i
).second
[j
]);
1280 assert(k
!= pred
.end());
1281 assert((*k
).second
!= 0);
1284 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1285 assert(head
.find((*i
).first
) != head
.end() ||
1286 pred
.find((*i
).first
) != pred
.end());
1287 assert((*i
).second
.size() >= 1);
1289 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1290 assert(head
.find((*i
).first
) != head
.end() ||
1291 pred
.find((*i
).first
) != pred
.end());
1292 for (int j
= 0; j
< (*i
).second
.size(); ++j
)
1293 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1294 pred
.find((*i
).second
[j
]) != pred
.end());
1298 max_term
* indicator::create_max_term(const indicator_term
*it
)
1300 int dim
= it
->den
.NumCols();
1301 int nparam
= ic
.P
->Dimension
- ic
.vertex
.length();
1302 max_term
*maximum
= new max_term
;
1303 maximum
->domain
= new EDomain(D
);
1304 for (int j
= 0; j
< dim
; ++j
) {
1305 evalue
*E
= new evalue
;
1307 evalue_copy(E
, it
->vertex
[j
]);
1308 if (evalue_frac2floor_in_domain(E
, D
->D
))
1310 maximum
->max
.push_back(E
);
1315 static order_sign
evalue_sign(evalue
*diff
, EDomain
*D
, barvinok_options
*options
)
1317 order_sign sign
= order_eq
;
1320 evalue_set_si(&mone
, -1, 1);
1321 int len
= 1 + D
->D
->Dimension
+ 1;
1322 Vector
*c
= Vector_Alloc(len
);
1323 Matrix
*T
= Matrix_Alloc(2, len
-1);
1325 int fract
= evalue2constraint(D
, diff
, c
->p
, len
);
1326 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1327 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1329 order_sign upper_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1330 if (upper_sign
== order_lt
|| !fract
)
1334 evalue2constraint(D
, diff
, c
->p
, len
);
1336 Vector_Copy(c
->p
+1, T
->p
[0], len
-1);
1337 value_assign(T
->p
[1][len
-2], c
->p
[0]);
1339 order_sign neg_lower_sign
= polyhedron_affine_sign(D
->D
, T
, options
);
1341 if (neg_lower_sign
== order_lt
)
1343 else if (neg_lower_sign
== order_eq
|| neg_lower_sign
== order_le
) {
1344 if (upper_sign
== order_eq
|| upper_sign
== order_le
)
1349 if (upper_sign
== order_lt
|| upper_sign
== order_le
||
1350 upper_sign
== order_eq
)
1353 sign
= order_unknown
;
1359 free_evalue_refs(&mone
);
1364 /* An auxiliary class that keeps a reference to an evalue
1365 * and frees it when it goes out of scope.
1367 struct temp_evalue
{
1369 temp_evalue() : E(NULL
) {}
1370 temp_evalue(evalue
*e
) : E(e
) {}
1371 operator evalue
* () const { return E
; }
1372 evalue
*operator=(evalue
*e
) {
1374 free_evalue_refs(E
);
1382 free_evalue_refs(E
);
1388 static void substitute(vector
<indicator_term
*>& term
, evalue
*equation
)
1390 evalue
*fract
= NULL
;
1391 evalue
*val
= new evalue
;
1393 evalue_copy(val
, equation
);
1396 value_init(factor
.d
);
1397 value_init(factor
.x
.n
);
1400 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= fractional
;
1401 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1404 if (value_zero_p(e
->d
) && e
->x
.p
->type
== fractional
)
1405 fract
= &e
->x
.p
->arr
[0];
1407 for (e
= val
; value_zero_p(e
->d
) && e
->x
.p
->type
!= polynomial
;
1408 e
= &e
->x
.p
->arr
[type_offset(e
->x
.p
)])
1410 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== polynomial
);
1413 int offset
= type_offset(e
->x
.p
);
1415 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].d
));
1416 assert(value_notzero_p(e
->x
.p
->arr
[offset
+1].x
.n
));
1417 if (value_neg_p(e
->x
.p
->arr
[offset
+1].x
.n
)) {
1418 value_oppose(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1419 value_assign(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1421 value_assign(factor
.d
, e
->x
.p
->arr
[offset
+1].x
.n
);
1422 value_oppose(factor
.x
.n
, e
->x
.p
->arr
[offset
+1].d
);
1425 free_evalue_refs(&e
->x
.p
->arr
[offset
+1]);
1428 *e
= e
->x
.p
->arr
[offset
];
1433 for (int i
= 0; i
< term
.size(); ++i
)
1434 term
[i
]->substitute(fract
, val
);
1436 free_evalue_refs(&p
->arr
[0]);
1438 for (int i
= 0; i
< term
.size(); ++i
)
1439 term
[i
]->substitute(p
->pos
, val
);
1442 free_evalue_refs(&factor
);
1443 free_evalue_refs(val
);
1449 order_sign
partial_order::compare(const indicator_term
*a
, const indicator_term
*b
)
1451 unsigned dim
= a
->den
.NumCols();
1452 order_sign sign
= order_eq
;
1453 EDomain
*D
= ind
->D
;
1454 unsigned MaxRays
= ind
->options
->verify
.barvinok
->MaxRays
;
1455 bool rational
= a
->sign
== 0 || b
->sign
== 0;
1457 order_sign cached_sign
= order_eq
;
1458 for (int k
= 0; k
< dim
; ++k
) {
1459 cached_sign
= evalue_diff_constant_sign(a
->vertex
[k
], b
->vertex
[k
]);
1460 if (cached_sign
!= order_eq
)
1463 if (cached_sign
!= order_undefined
)
1466 order_cache_el cache_el
;
1467 cached_sign
= order_undefined
;
1469 cached_sign
= cache
.check(a
, b
, cache_el
);
1470 if (cached_sign
!= order_undefined
) {
1475 if (rational
&& POL_ISSET(MaxRays
, POL_INTEGER
)) {
1476 ind
->options
->verify
.barvinok
->MaxRays
&= ~POL_INTEGER
;
1477 if (ind
->options
->verify
.barvinok
->MaxRays
)
1478 ind
->options
->verify
.barvinok
->MaxRays
|= POL_HIGH_BIT
;
1483 vector
<indicator_term
*> term
;
1485 for (int k
= 0; k
< dim
; ++k
) {
1486 /* compute a->vertex[k] - b->vertex[k] */
1488 if (cache_el
.e
.size() <= k
) {
1489 diff
= ediff(a
->vertex
[k
], b
->vertex
[k
]);
1490 cache_el
.e
.push_back(diff
);
1492 diff
= cache_el
.e
[k
];
1495 tdiff
= diff
= ediff(term
[0]->vertex
[k
], term
[1]->vertex
[k
]);
1496 order_sign diff_sign
;
1498 diff_sign
= order_undefined
;
1499 else if (eequal(a
->vertex
[k
], b
->vertex
[k
]))
1500 diff_sign
= order_eq
;
1502 diff_sign
= evalue_sign(diff
, D
, ind
->options
->verify
.barvinok
);
1504 if (diff_sign
== order_undefined
) {
1505 assert(sign
== order_le
|| sign
== order_ge
);
1506 if (sign
== order_le
)
1512 if (diff_sign
== order_lt
) {
1513 if (sign
== order_eq
|| sign
== order_le
)
1516 sign
= order_unknown
;
1519 if (diff_sign
== order_gt
) {
1520 if (sign
== order_eq
|| sign
== order_ge
)
1523 sign
= order_unknown
;
1526 if (diff_sign
== order_eq
) {
1527 if (D
== ind
->D
&& term
.size() == 0 && !rational
&&
1528 !EVALUE_IS_ZERO(*diff
))
1529 ind
->add_substitution(diff
);
1532 if ((diff_sign
== order_unknown
) ||
1533 ((diff_sign
== order_lt
|| diff_sign
== order_le
) && sign
== order_ge
) ||
1534 ((diff_sign
== order_gt
|| diff_sign
== order_ge
) && sign
== order_le
)) {
1535 sign
= order_unknown
;
1542 term
.push_back(new indicator_term(*a
));
1543 term
.push_back(new indicator_term(*b
));
1545 substitute(term
, diff
);
1549 cache
.add(cache_el
, sign
);
1553 if (D
&& D
!= ind
->D
)
1561 ind
->options
->verify
.barvinok
->MaxRays
= MaxRays
;
1565 bool partial_order::compared(const indicator_term
* a
, const indicator_term
* b
)
1567 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator j
;
1570 if (j
!= lt
.end() && find(lt
[a
].begin(), lt
[a
].end(), b
) != lt
[a
].end())
1574 if (j
!= le
.end() && find(le
[a
].begin(), le
[a
].end(), b
) != le
[a
].end())
1580 void partial_order::add(const indicator_term
* it
,
1581 std::set
<const indicator_term
*> *filter
)
1583 if (eq
.find(it
) != eq
.end() && eq
[it
].size() == 1)
1586 typeof(head
) head_copy(head
);
1591 std::set
<const indicator_term
*>::iterator i
;
1592 for (i
= head_copy
.begin(); i
!= head_copy
.end(); ++i
) {
1595 if (eq
.find(*i
) != eq
.end() && eq
[*i
].size() == 1)
1598 if (filter
->find(*i
) == filter
->end())
1600 if (compared(*i
, it
))
1603 order_sign sign
= compare(it
, *i
);
1604 if (sign
== order_lt
) {
1605 lt
[it
].push_back(*i
);
1607 } else if (sign
== order_le
) {
1608 le
[it
].push_back(*i
);
1610 } else if (sign
== order_eq
) {
1613 } else if (sign
== order_gt
) {
1614 pending
[*i
].push_back(it
);
1615 lt
[*i
].push_back(it
);
1617 } else if (sign
== order_ge
) {
1618 pending
[*i
].push_back(it
);
1619 le
[*i
].push_back(it
);
1625 void partial_order::remove(const indicator_term
* it
)
1627 std::set
<const indicator_term
*> filter
;
1628 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1630 assert(head
.find(it
) != head
.end());
1633 if (i
!= eq
.end()) {
1634 assert(eq
[it
].size() >= 1);
1635 const indicator_term
*base
;
1636 if (eq
[it
].size() == 1) {
1640 vector
<const indicator_term
* >::iterator j
;
1641 j
= find(eq
[base
].begin(), eq
[base
].end(), it
);
1642 assert(j
!= eq
[base
].end());
1645 /* "it" may no longer be the smallest, since the order
1646 * structure may have been copied from another one.
1648 sort(eq
[it
].begin()+1, eq
[it
].end(), pred
.key_comp());
1649 assert(eq
[it
][0] == it
);
1650 eq
[it
].erase(eq
[it
].begin());
1655 for (int j
= 1; j
< eq
[base
].size(); ++j
)
1656 eq
[eq
[base
][j
]][0] = base
;
1659 if (i
!= lt
.end()) {
1665 if (i
!= le
.end()) {
1670 i
= pending
.find(it
);
1671 if (i
!= pending
.end()) {
1672 pending
[base
] = pending
[it
];
1677 if (eq
[base
].size() == 1)
1686 if (i
!= lt
.end()) {
1687 for (int j
= 0; j
< lt
[it
].size(); ++j
) {
1688 filter
.insert(lt
[it
][j
]);
1689 dec_pred(lt
[it
][j
]);
1695 if (i
!= le
.end()) {
1696 for (int j
= 0; j
< le
[it
].size(); ++j
) {
1697 filter
.insert(le
[it
][j
]);
1698 dec_pred(le
[it
][j
]);
1705 i
= pending
.find(it
);
1706 if (i
!= pending
.end()) {
1707 vector
<const indicator_term
*> it_pending
= pending
[it
];
1709 for (int j
= 0; j
< it_pending
.size(); ++j
) {
1710 filter
.erase(it_pending
[j
]);
1711 add(it_pending
[j
], &filter
);
1716 void partial_order::print(ostream
& os
, char **p
)
1718 map
<const indicator_term
*, vector
<const indicator_term
* > >::iterator i
;
1719 map
<const indicator_term
*, int >::iterator j
;
1720 std::set
<const indicator_term
*>::iterator k
;
1721 for (k
= head
.begin(); k
!= head
.end(); ++k
) {
1725 for (j
= pred
.begin(); j
!= pred
.end(); ++j
) {
1726 (*j
).first
->print(os
, p
);
1727 os
<< ": " << (*j
).second
<< endl
;
1729 for (i
= lt
.begin(); i
!= lt
.end(); ++i
) {
1730 (*i
).first
->print(os
, p
);
1731 assert(head
.find((*i
).first
) != head
.end() ||
1732 pred
.find((*i
).first
) != pred
.end());
1733 if (pred
.find((*i
).first
) != pred
.end())
1734 os
<< "(" << pred
[(*i
).first
] << ")";
1736 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1739 (*i
).second
[j
]->print(os
, p
);
1740 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1741 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1745 for (i
= le
.begin(); i
!= le
.end(); ++i
) {
1746 (*i
).first
->print(os
, p
);
1747 assert(head
.find((*i
).first
) != head
.end() ||
1748 pred
.find((*i
).first
) != pred
.end());
1749 if (pred
.find((*i
).first
) != pred
.end())
1750 os
<< "(" << pred
[(*i
).first
] << ")";
1752 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1755 (*i
).second
[j
]->print(os
, p
);
1756 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1757 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1761 for (i
= eq
.begin(); i
!= eq
.end(); ++i
) {
1762 if ((*i
).second
.size() <= 1)
1764 (*i
).first
->print(os
, p
);
1765 assert(head
.find((*i
).first
) != head
.end() ||
1766 pred
.find((*i
).first
) != pred
.end());
1767 if (pred
.find((*i
).first
) != pred
.end())
1768 os
<< "(" << pred
[(*i
).first
] << ")";
1769 for (int j
= 1; j
< (*i
).second
.size(); ++j
) {
1772 (*i
).second
[j
]->print(os
, p
);
1773 assert(head
.find((*i
).second
[j
]) != head
.end() ||
1774 pred
.find((*i
).second
[j
]) != pred
.end());
1775 if (pred
.find((*i
).second
[j
]) != pred
.end())
1776 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1780 for (i
= pending
.begin(); i
!= pending
.end(); ++i
) {
1781 os
<< "pending on ";
1782 (*i
).first
->print(os
, p
);
1783 assert(head
.find((*i
).first
) != head
.end() ||
1784 pred
.find((*i
).first
) != pred
.end());
1785 if (pred
.find((*i
).first
) != pred
.end())
1786 os
<< "(" << pred
[(*i
).first
] << ")";
1788 for (int j
= 0; j
< (*i
).second
.size(); ++j
) {
1791 (*i
).second
[j
]->print(os
, p
);
1792 assert(pred
.find((*i
).second
[j
]) != pred
.end());
1793 os
<< "(" << pred
[(*i
).second
[j
]] << ")";
1799 void indicator::add(const indicator_term
* it
)
1801 indicator_term
*nt
= new indicator_term(*it
);
1802 if (options
->reduce
)
1803 nt
->reduce_in_domain(P
? P
: D
->D
);
1805 order
.add(nt
, NULL
);
1806 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1809 void indicator::remove(const indicator_term
* it
)
1811 vector
<indicator_term
*>::iterator i
;
1812 i
= find(term
.begin(), term
.end(), it
);
1813 assert(i
!= term
.end());
1816 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1820 void indicator::expand_rational_vertex(const indicator_term
*initial
)
1822 int pos
= initial
->pos
;
1824 if (ic
.terms
[pos
].size() == 0) {
1826 FORALL_PVertex_in_ParamPolyhedron(V
, PD
, ic
.PP
) // _i is internal counter
1828 ic
.decompose_at_vertex(V
, pos
, options
->verify
.barvinok
);
1831 END_FORALL_PVertex_in_ParamPolyhedron
;
1833 for (int j
= 0; j
< ic
.terms
[pos
].size(); ++j
)
1834 add(ic
.terms
[pos
][j
]);
1837 void indicator::remove_initial_rational_vertices()
1840 const indicator_term
*initial
= NULL
;
1841 std::set
<const indicator_term
*>::iterator i
;
1842 for (i
= order
.head
.begin(); i
!= order
.head
.end(); ++i
) {
1843 if ((*i
)->sign
!= 0)
1845 if (order
.eq
.find(*i
) != order
.eq
.end() && order
.eq
[*i
].size() <= 1)
1852 expand_rational_vertex(initial
);
1856 void indicator::reduce_in_domain(Polyhedron
*D
)
1858 for (int i
= 0; i
< term
.size(); ++i
)
1859 term
[i
]->reduce_in_domain(D
);
1862 void indicator::print(ostream
& os
, char **p
)
1864 assert(term
.size() == order
.head
.size() + order
.pred
.size());
1865 for (int i
= 0; i
< term
.size(); ++i
) {
1866 term
[i
]->print(os
, p
);
1868 os
<< ": " << term
[i
]->eval(D
->sample
->p
);
1875 /* Remove pairs of opposite terms */
1876 void indicator::simplify()
1878 for (int i
= 0; i
< term
.size(); ++i
) {
1879 for (int j
= i
+1; j
< term
.size(); ++j
) {
1880 if (term
[i
]->sign
+ term
[j
]->sign
!= 0)
1882 if (term
[i
]->den
!= term
[j
]->den
)
1885 for (k
= 0; k
< term
[i
]->den
.NumCols(); ++k
)
1886 if (!eequal(term
[i
]->vertex
[k
], term
[j
]->vertex
[k
]))
1888 if (k
< term
[i
]->den
.NumCols())
1892 term
.erase(term
.begin()+j
);
1893 term
.erase(term
.begin()+i
);
1900 void indicator::peel(int i
, int j
)
1908 int dim
= term
[i
]->den
.NumCols();
1913 int n_common
= 0, n_i
= 0, n_j
= 0;
1915 common
.SetDims(min(term
[i
]->den
.NumRows(), term
[j
]->den
.NumRows()), dim
);
1916 rest_i
.SetDims(term
[i
]->den
.NumRows(), dim
);
1917 rest_j
.SetDims(term
[j
]->den
.NumRows(), dim
);
1920 for (k
= 0, l
= 0; k
< term
[i
]->den
.NumRows() && l
< term
[j
]->den
.NumRows(); ) {
1921 int s
= lex_cmp(term
[i
]->den
[k
], term
[j
]->den
[l
]);
1923 common
[n_common
++] = term
[i
]->den
[k
];
1927 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1929 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1931 while (k
< term
[i
]->den
.NumRows())
1932 rest_i
[n_i
++] = term
[i
]->den
[k
++];
1933 while (l
< term
[j
]->den
.NumRows())
1934 rest_j
[n_j
++] = term
[j
]->den
[l
++];
1935 common
.SetDims(n_common
, dim
);
1936 rest_i
.SetDims(n_i
, dim
);
1937 rest_j
.SetDims(n_j
, dim
);
1939 for (k
= 0; k
<= n_i
; ++k
) {
1940 indicator_term
*it
= new indicator_term(*term
[i
]);
1941 it
->den
.SetDims(n_common
+ k
, dim
);
1942 for (l
= 0; l
< n_common
; ++l
)
1943 it
->den
[l
] = common
[l
];
1944 for (l
= 0; l
< k
; ++l
)
1945 it
->den
[n_common
+l
] = rest_i
[l
];
1946 lex_order_rows(it
->den
);
1948 for (l
= 0; l
< dim
; ++l
)
1949 evalue_add_constant(it
->vertex
[l
], rest_i
[k
-1][l
]);
1953 for (k
= 0; k
<= n_j
; ++k
) {
1954 indicator_term
*it
= new indicator_term(*term
[j
]);
1955 it
->den
.SetDims(n_common
+ k
, dim
);
1956 for (l
= 0; l
< n_common
; ++l
)
1957 it
->den
[l
] = common
[l
];
1958 for (l
= 0; l
< k
; ++l
)
1959 it
->den
[n_common
+l
] = rest_j
[l
];
1960 lex_order_rows(it
->den
);
1962 for (l
= 0; l
< dim
; ++l
)
1963 evalue_add_constant(it
->vertex
[l
], rest_j
[k
-1][l
]);
1966 term
.erase(term
.begin()+j
);
1967 term
.erase(term
.begin()+i
);
1970 void indicator::combine(const indicator_term
*a
, const indicator_term
*b
)
1972 int dim
= a
->den
.NumCols();
1975 mat_ZZ rest_i
; /* factors in a, but not in b */
1976 mat_ZZ rest_j
; /* factors in b, but not in a */
1977 int n_common
= 0, n_i
= 0, n_j
= 0;
1979 common
.SetDims(min(a
->den
.NumRows(), b
->den
.NumRows()), dim
);
1980 rest_i
.SetDims(a
->den
.NumRows(), dim
);
1981 rest_j
.SetDims(b
->den
.NumRows(), dim
);
1984 for (k
= 0, l
= 0; k
< a
->den
.NumRows() && l
< b
->den
.NumRows(); ) {
1985 int s
= lex_cmp(a
->den
[k
], b
->den
[l
]);
1987 common
[n_common
++] = a
->den
[k
];
1991 rest_i
[n_i
++] = a
->den
[k
++];
1993 rest_j
[n_j
++] = b
->den
[l
++];
1995 while (k
< a
->den
.NumRows())
1996 rest_i
[n_i
++] = a
->den
[k
++];
1997 while (l
< b
->den
.NumRows())
1998 rest_j
[n_j
++] = b
->den
[l
++];
1999 common
.SetDims(n_common
, dim
);
2000 rest_i
.SetDims(n_i
, dim
);
2001 rest_j
.SetDims(n_j
, dim
);
2003 assert(order
.eq
[a
].size() > 1);
2004 indicator_term
*prev
;
2007 for (int k
= n_i
-1; k
>= 0; --k
) {
2008 indicator_term
*it
= new indicator_term(*b
);
2009 it
->den
.SetDims(n_common
+ n_j
+ n_i
-k
, dim
);
2010 for (int l
= k
; l
< n_i
; ++l
)
2011 it
->den
[n_common
+n_j
+l
-k
] = rest_i
[l
];
2012 lex_order_rows(it
->den
);
2013 for (int m
= 0; m
< dim
; ++m
)
2014 evalue_add_constant(it
->vertex
[m
], rest_i
[k
][m
]);
2015 it
->sign
= -it
->sign
;
2017 order
.pending
[it
].push_back(prev
);
2018 order
.lt
[it
].push_back(prev
);
2019 order
.inc_pred(prev
);
2022 order
.head
.insert(it
);
2026 indicator_term
*it
= new indicator_term(*b
);
2027 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
2028 for (l
= 0; l
< n_i
; ++l
)
2029 it
->den
[n_common
+n_j
+l
] = rest_i
[l
];
2030 lex_order_rows(it
->den
);
2032 order
.pending
[a
].push_back(prev
);
2033 order
.lt
[a
].push_back(prev
);
2034 order
.inc_pred(prev
);
2035 order
.replace(b
, it
);
2040 for (int k
= n_j
-1; k
>= 0; --k
) {
2041 indicator_term
*it
= new indicator_term(*a
);
2042 it
->den
.SetDims(n_common
+ n_i
+ n_j
-k
, dim
);
2043 for (int l
= k
; l
< n_j
; ++l
)
2044 it
->den
[n_common
+n_i
+l
-k
] = rest_j
[l
];
2045 lex_order_rows(it
->den
);
2046 for (int m
= 0; m
< dim
; ++m
)
2047 evalue_add_constant(it
->vertex
[m
], rest_j
[k
][m
]);
2048 it
->sign
= -it
->sign
;
2050 order
.pending
[it
].push_back(prev
);
2051 order
.lt
[it
].push_back(prev
);
2052 order
.inc_pred(prev
);
2055 order
.head
.insert(it
);
2059 indicator_term
*it
= new indicator_term(*a
);
2060 it
->den
.SetDims(n_common
+ n_i
+ n_j
, dim
);
2061 for (l
= 0; l
< n_j
; ++l
)
2062 it
->den
[n_common
+n_i
+l
] = rest_j
[l
];
2063 lex_order_rows(it
->den
);
2065 order
.pending
[a
].push_back(prev
);
2066 order
.lt
[a
].push_back(prev
);
2067 order
.inc_pred(prev
);
2068 order
.replace(a
, it
);
2072 assert(term
.size() == order
.head
.size() + order
.pred
.size());
2075 bool indicator::handle_equal_numerators(const indicator_term
*base
)
2077 for (int i
= 0; i
< order
.eq
[base
].size(); ++i
) {
2078 for (int j
= i
+1; j
< order
.eq
[base
].size(); ++j
) {
2079 if (order
.eq
[base
][i
]->is_opposite(order
.eq
[base
][j
])) {
2080 remove(order
.eq
[base
][j
]);
2081 remove(i
? order
.eq
[base
][i
] : base
);
2086 for (int j
= 1; j
< order
.eq
[base
].size(); ++j
)
2087 if (order
.eq
[base
][j
]->sign
!= base
->sign
) {
2088 combine(base
, order
.eq
[base
][j
]);
2094 void indicator::substitute(evalue
*equation
)
2096 ::substitute(term
, equation
);
2099 void indicator::add_substitution(evalue
*equation
)
2101 for (int i
= 0; i
< substitutions
.size(); ++i
)
2102 if (eequal(substitutions
[i
], equation
))
2104 evalue
*copy
= new evalue();
2105 value_init(copy
->d
);
2106 evalue_copy(copy
, equation
);
2107 substitutions
.push_back(copy
);
2110 void indicator::perform_pending_substitutions()
2112 if (substitutions
.size() == 0)
2115 for (int i
= 0; i
< substitutions
.size(); ++i
) {
2116 substitute(substitutions
[i
]);
2117 free_evalue_refs(substitutions
[i
]);
2118 delete substitutions
[i
];
2120 substitutions
.clear();
2124 static void print_varlist(ostream
& os
, int n
, char **names
)
2128 for (i
= 0; i
< n
; ++i
) {
2136 void max_term::print(ostream
& os
, char **p
, barvinok_options
*options
) const
2139 print_varlist(os
, domain
->dimension(), p
);
2142 for (int i
= 0; i
< max
.size(); ++i
) {
2145 evalue_print(os
, max
[i
], p
);
2149 domain
->print_constraints(os
, p
, options
);
2153 /* T maps the compressed parameters to the original parameters,
2154 * while this max_term is based on the compressed parameters
2155 * and we want get the original parameters back.
2157 void max_term::substitute(Matrix
*T
, barvinok_options
*options
)
2159 assert(domain
->dimension() == T
->NbColumns
-1);
2160 int nexist
= domain
->D
->Dimension
- (T
->NbColumns
-1);
2162 Matrix
*inv
= left_inverse(T
, &Eq
);
2165 value_init(denom
.d
);
2166 value_init(denom
.x
.n
);
2167 value_set_si(denom
.x
.n
, 1);
2168 value_assign(denom
.d
, inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
2171 v
.SetLength(inv
->NbColumns
);
2172 evalue
* subs
[inv
->NbRows
-1];
2173 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2174 values2zz(inv
->p
[i
], v
, v
.length());
2175 subs
[i
] = multi_monom(v
);
2176 emul(&denom
, subs
[i
]);
2178 free_evalue_refs(&denom
);
2180 domain
->substitute(subs
, inv
, Eq
, options
->MaxRays
);
2183 for (int i
= 0; i
< max
.size(); ++i
) {
2184 evalue_substitute(max
[i
], subs
);
2185 reduce_evalue(max
[i
]);
2188 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
2189 free_evalue_refs(subs
[i
]);
2195 int Last_Non_Zero(Value
*p
, unsigned len
)
2197 for (int i
= len
-1; i
>= 0; --i
)
2198 if (value_notzero_p(p
[i
]))
2203 static void SwapColumns(Value
**V
, int n
, int i
, int j
)
2205 for (int r
= 0; r
< n
; ++r
)
2206 value_swap(V
[r
][i
], V
[r
][j
]);
2209 static void SwapColumns(Polyhedron
*P
, int i
, int j
)
2211 SwapColumns(P
->Constraint
, P
->NbConstraints
, i
, j
);
2212 SwapColumns(P
->Ray
, P
->NbRays
, i
, j
);
2215 Vector
*max_term::eval(Value
*val
, unsigned MaxRays
) const
2217 if (!domain
->contains(val
, domain
->dimension()))
2219 Vector
*res
= Vector_Alloc(max
.size());
2220 for (int i
= 0; i
< max
.size(); ++i
) {
2221 compute_evalue(max
[i
], val
, &res
->p
[i
]);
2228 enum sign
{ le
, ge
, lge
} sign
;
2230 split (evalue
*c
, enum sign s
) : constraint(c
), sign(s
) {}
2233 static void split_on(const split
& sp
, EDomain
*D
,
2234 EDomain
**Dlt
, EDomain
**Deq
, EDomain
**Dgt
,
2235 lexmin_options
*options
)
2241 ge_constraint
*ge
= D
->compute_ge_constraint(sp
.constraint
);
2242 if (sp
.sign
== split::lge
|| sp
.sign
== split::ge
)
2243 ED
[2] = EDomain::new_from_ge_constraint(ge
, 1, options
->verify
.barvinok
);
2246 if (sp
.sign
== split::lge
|| sp
.sign
== split::le
)
2247 ED
[0] = EDomain::new_from_ge_constraint(ge
, -1, options
->verify
.barvinok
);
2251 assert(sp
.sign
== split::lge
|| sp
.sign
== split::ge
|| sp
.sign
== split::le
);
2252 ED
[1] = EDomain::new_from_ge_constraint(ge
, 0, options
->verify
.barvinok
);
2256 for (int i
= 0; i
< 3; ++i
) {
2259 if (D
->sample
&& ED
[i
]->contains(D
->sample
->p
, D
->sample
->Size
-1)) {
2260 ED
[i
]->sample
= Vector_Alloc(D
->sample
->Size
);
2261 Vector_Copy(D
->sample
->p
, ED
[i
]->sample
->p
, D
->sample
->Size
);
2262 } else if (emptyQ2(ED
[i
]->D
) ||
2263 (options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2264 !(ED
[i
]->not_empty(options
)))) {
2274 ostream
& operator<< (ostream
& os
, const vector
<int> & v
)
2277 for (int i
= 0; i
< v
.size(); ++i
) {
2286 static bool isTranslation(Matrix
*M
)
2289 if (M
->NbRows
!= M
->NbColumns
)
2292 for (i
= 0;i
< M
->NbRows
; i
++)
2293 for (j
= 0; j
< M
->NbColumns
-1; j
++)
2295 if(value_notone_p(M
->p
[i
][j
]))
2298 if(value_notzero_p(M
->p
[i
][j
]))
2301 return value_one_p(M
->p
[M
->NbRows
-1][M
->NbColumns
-1]);
2304 static Matrix
*compress_parameters(Polyhedron
**P
, Polyhedron
**C
,
2305 unsigned nparam
, unsigned MaxRays
)
2309 /* compress_parms doesn't like equalities that only involve parameters */
2310 for (int i
= 0; i
< (*P
)->NbEq
; ++i
)
2311 assert(First_Non_Zero((*P
)->Constraint
[i
]+1, (*P
)->Dimension
-nparam
) != -1);
2313 M
= Matrix_Alloc((*P
)->NbEq
, (*P
)->Dimension
+2);
2314 Vector_Copy((*P
)->Constraint
[0], M
->p
[0], (*P
)->NbEq
* ((*P
)->Dimension
+2));
2315 CP
= compress_parms(M
, nparam
);
2318 if (isTranslation(CP
)) {
2323 T
= align_matrix(CP
, (*P
)->Dimension
+1);
2324 *P
= Polyhedron_Preimage(*P
, T
, MaxRays
);
2327 *C
= Polyhedron_Preimage(*C
, CP
, MaxRays
);
2332 void construct_rational_vertices(Param_Polyhedron
*PP
, Matrix
*T
, unsigned dim
,
2333 int nparam
, vector
<indicator_term
*>& vertices
)
2342 v
.SetLength(nparam
+1);
2345 value_init(factor
.d
);
2346 value_init(factor
.x
.n
);
2347 value_set_si(factor
.x
.n
, 1);
2348 value_set_si(factor
.d
, 1);
2350 for (i
= 0, PV
= PP
->V
; PV
; ++i
, PV
= PV
->next
) {
2351 indicator_term
*term
= new indicator_term(dim
, i
);
2352 vertices
.push_back(term
);
2353 Matrix
*M
= Matrix_Alloc(PV
->Vertex
->NbRows
+nparam
+1, nparam
+1);
2354 value_set_si(lcm
, 1);
2355 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
)
2356 value_lcm(lcm
, PV
->Vertex
->p
[j
][nparam
+1], &lcm
);
2357 value_assign(M
->p
[M
->NbRows
-1][M
->NbColumns
-1], lcm
);
2358 for (int j
= 0; j
< PV
->Vertex
->NbRows
; ++j
) {
2359 value_division(tmp
, lcm
, PV
->Vertex
->p
[j
][nparam
+1]);
2360 Vector_Scale(PV
->Vertex
->p
[j
], M
->p
[j
], tmp
, nparam
+1);
2362 for (int j
= 0; j
< nparam
; ++j
)
2363 value_assign(M
->p
[PV
->Vertex
->NbRows
+j
][j
], lcm
);
2365 Matrix
*M2
= Matrix_Alloc(T
->NbRows
, M
->NbColumns
);
2366 Matrix_Product(T
, M
, M2
);
2370 for (int j
= 0; j
< dim
; ++j
) {
2371 values2zz(M
->p
[j
], v
, nparam
+1);
2372 term
->vertex
[j
] = multi_monom(v
);
2373 value_assign(factor
.d
, lcm
);
2374 emul(&factor
, term
->vertex
[j
]);
2378 assert(i
== PP
->nbV
);
2379 free_evalue_refs(&factor
);
2384 static vector
<max_term
*> lexmin(indicator
& ind
, unsigned nparam
,
2387 vector
<max_term
*> maxima
;
2388 std::set
<const indicator_term
*>::iterator i
;
2389 vector
<int> best_score
;
2390 vector
<int> second_score
;
2391 vector
<int> neg_score
;
2394 ind
.perform_pending_substitutions();
2395 const indicator_term
*best
= NULL
, *second
= NULL
, *neg
= NULL
,
2396 *neg_eq
= NULL
, *neg_le
= NULL
;
2397 for (i
= ind
.order
.head
.begin(); i
!= ind
.order
.head
.end(); ++i
) {
2399 const indicator_term
*term
= *i
;
2400 if (term
->sign
== 0) {
2401 ind
.expand_rational_vertex(term
);
2405 if (ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2407 if (ind
.order
.eq
[term
].size() <= 1)
2409 for (j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2410 if (ind
.order
.pred
.find(ind
.order
.eq
[term
][j
]) !=
2411 ind
.order
.pred
.end())
2413 if (j
< ind
.order
.eq
[term
].size())
2415 score
.push_back(ind
.order
.eq
[term
].size());
2418 if (ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2419 score
.push_back(ind
.order
.le
[term
].size());
2422 if (ind
.order
.lt
.find(term
) != ind
.order
.lt
.end())
2423 score
.push_back(-ind
.order
.lt
[term
].size());
2427 if (term
->sign
> 0) {
2428 if (!best
|| score
< best_score
) {
2430 second_score
= best_score
;
2433 } else if (!second
|| score
< second_score
) {
2435 second_score
= score
;
2438 if (!neg_eq
&& ind
.order
.eq
.find(term
) != ind
.order
.eq
.end()) {
2439 for (int j
= 1; j
< ind
.order
.eq
[term
].size(); ++j
)
2440 if (ind
.order
.eq
[term
][j
]->sign
!= term
->sign
) {
2445 if (!neg_le
&& ind
.order
.le
.find(term
) != ind
.order
.le
.end())
2447 if (!neg
|| score
< neg_score
) {
2453 if (i
!= ind
.order
.head
.end())
2456 if (!best
&& neg_eq
) {
2457 assert(ind
.order
.eq
[neg_eq
].size() != 0);
2458 bool handled
= ind
.handle_equal_numerators(neg_eq
);
2463 if (!best
&& neg_le
) {
2464 /* The smallest term is negative and <= some positive term */
2470 /* apparently there can be negative initial term on empty domains */
2471 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2472 ind
.options
->verify
.barvinok
->lp_solver
== BV_LP_POLYLIB
)
2477 if (!second
&& !neg
) {
2478 const indicator_term
*rat
= NULL
;
2480 if (ind
.order
.le
.find(best
) == ind
.order
.le
.end()) {
2481 if (ind
.order
.eq
.find(best
) != ind
.order
.eq
.end()) {
2482 bool handled
= ind
.handle_equal_numerators(best
);
2483 if (ind
.options
->emptiness_check
!=
2484 BV_LEXMIN_EMPTINESS_CHECK_NONE
&&
2485 ind
.options
->verify
.barvinok
->lp_solver
== BV_LP_POLYLIB
)
2487 /* If !handled then the leading coefficient is bigger than one;
2488 * must be an empty domain
2495 maxima
.push_back(ind
.create_max_term(best
));
2498 for (int j
= 0; j
< ind
.order
.le
[best
].size(); ++j
) {
2499 if (ind
.order
.le
[best
][j
]->sign
== 0) {
2500 if (!rat
&& ind
.order
.pred
[ind
.order
.le
[best
][j
]] == 1)
2501 rat
= ind
.order
.le
[best
][j
];
2502 } else if (ind
.order
.le
[best
][j
]->sign
> 0) {
2503 second
= ind
.order
.le
[best
][j
];
2506 neg
= ind
.order
.le
[best
][j
];
2509 if (!second
&& !neg
) {
2511 ind
.order
.unset_le(best
, rat
);
2512 ind
.expand_rational_vertex(rat
);
2519 ind
.order
.unset_le(best
, second
);
2525 unsigned dim
= best
->den
.NumCols();
2528 for (int k
= 0; k
< dim
; ++k
) {
2529 diff
= ediff(best
->vertex
[k
], second
->vertex
[k
]);
2530 sign
= evalue_sign(diff
, ind
.D
, ind
.options
->verify
.barvinok
);
2532 /* neg can never be smaller than best, unless it may still cancel.
2533 * This can happen if positive terms have been determined to
2534 * be equal or less than or equal to some negative term.
2536 if (second
== neg
&& !neg_eq
&& !neg_le
) {
2537 if (sign
== order_ge
)
2539 if (sign
== order_unknown
)
2543 if (sign
!= order_eq
)
2545 if (!EVALUE_IS_ZERO(*diff
)) {
2546 ind
.add_substitution(diff
);
2547 ind
.perform_pending_substitutions();
2550 if (sign
== order_eq
) {
2551 ind
.order
.set_equal(best
, second
);
2554 if (sign
== order_lt
) {
2555 ind
.order
.lt
[best
].push_back(second
);
2556 ind
.order
.inc_pred(second
);
2559 if (sign
== order_gt
) {
2560 ind
.order
.lt
[second
].push_back(best
);
2561 ind
.order
.inc_pred(best
);
2565 split
sp(diff
, sign
== order_le
? split::le
:
2566 sign
== order_ge
? split::ge
: split::lge
);
2568 EDomain
*Dlt
, *Deq
, *Dgt
;
2569 split_on(sp
, ind
.D
, &Dlt
, &Deq
, &Dgt
, ind
.options
);
2570 if (ind
.options
->emptiness_check
!= BV_LEXMIN_EMPTINESS_CHECK_NONE
)
2571 assert(Dlt
|| Deq
|| Dgt
);
2572 else if (!(Dlt
|| Deq
|| Dgt
))
2573 /* Must have been empty all along */
2576 if (Deq
&& (Dlt
|| Dgt
)) {
2577 int locsize
= loc
.size();
2579 indicator
indeq(ind
, Deq
);
2581 indeq
.add_substitution(diff
);
2582 indeq
.perform_pending_substitutions();
2583 vector
<max_term
*> maxeq
= lexmin(indeq
, nparam
, loc
);
2584 maxima
.insert(maxima
.end(), maxeq
.begin(), maxeq
.end());
2585 loc
.resize(locsize
);
2588 int locsize
= loc
.size();
2590 indicator
indgt(ind
, Dgt
);
2592 /* we don't know the new location of these terms in indgt */
2594 indgt.order.lt[second].push_back(best);
2595 indgt.order.inc_pred(best);
2597 vector
<max_term
*> maxgt
= lexmin(indgt
, nparam
, loc
);
2598 maxima
.insert(maxima
.end(), maxgt
.begin(), maxgt
.end());
2599 loc
.resize(locsize
);
2604 ind
.set_domain(Deq
);
2605 ind
.add_substitution(diff
);
2606 ind
.perform_pending_substitutions();
2610 ind
.set_domain(Dlt
);
2611 ind
.order
.lt
[best
].push_back(second
);
2612 ind
.order
.inc_pred(second
);
2616 ind
.set_domain(Dgt
);
2617 ind
.order
.lt
[second
].push_back(best
);
2618 ind
.order
.inc_pred(best
);
2625 static vector
<max_term
*> lexmin(Polyhedron
*P
, Polyhedron
*C
,
2626 lexmin_options
*options
)
2628 unsigned nparam
= C
->Dimension
;
2629 Param_Polyhedron
*PP
= NULL
;
2630 Matrix
*T
= NULL
, *CP
= NULL
;
2631 Polyhedron
*Porig
= P
;
2632 Polyhedron
*Corig
= C
;
2633 vector
<max_term
*> all_max
;
2635 unsigned P2PSD_MaxRays
;
2640 POL_ENSURE_VERTICES(P
);
2645 assert(P
->NbBid
== 0);
2648 remove_all_equalities(&P
, &C
, &CP
, &T
, nparam
,
2649 options
->verify
.barvinok
->MaxRays
);
2651 nparam
= CP
->NbColumns
-1;
2659 if (options
->verify
.barvinok
->MaxRays
& POL_NO_DUAL
)
2662 P2PSD_MaxRays
= options
->verify
.barvinok
->MaxRays
;
2664 PP
= Polyhedron2Param_Domain(P
, C
, P2PSD_MaxRays
);
2666 unsigned dim
= P
->Dimension
- nparam
;
2670 indicator_constructor
ic(P
, dim
, PP
, T
);
2672 vector
<indicator_term
*> all_vertices
;
2673 construct_rational_vertices(PP
, T
, T
? T
->NbRows
-nparam
-1 : dim
,
2674 nparam
, all_vertices
);
2676 Polyhedron
*TC
= true_context(P
, C
, options
->verify
.barvinok
->MaxRays
);
2677 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
->verify
.barvinok
, i
, D
, rVD
)
2680 EDomain
*epVD
= new EDomain(rVD
);
2681 indicator
ind(ic
, D
, epVD
, options
);
2683 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
2684 ind
.add(all_vertices
[_i
]);
2685 END_FORALL_PVertex_in_ParamPolyhedron
;
2687 ind
.remove_initial_rational_vertices();
2690 vector
<max_term
*> maxima
= lexmin(ind
, nparam
, loc
);
2692 for (int j
= 0; j
< maxima
.size(); ++j
)
2693 maxima
[j
]->substitute(CP
, options
->verify
.barvinok
);
2694 all_max
.insert(all_max
.end(), maxima
.begin(), maxima
.end());
2697 END_FORALL_REDUCED_DOMAIN
2698 Polyhedron_Free(TC
);
2699 for (int i
= 0; i
< all_vertices
.size(); ++i
)
2700 delete all_vertices
[i
];
2705 Param_Polyhedron_Free(PP
);
2714 static void verify_results(Polyhedron
*A
, Polyhedron
*C
,
2715 vector
<max_term
*>& maxima
,
2716 struct verify_options
*options
);
2718 int main(int argc
, char **argv
)
2723 char **iter_names
, **param_names
;
2724 int print_solution
= 1;
2725 struct lexmin_options options
;
2726 static struct argp_child argp_children
[] = {
2727 { &barvinok_argp
, 0, 0, 0 },
2728 { &verify_argp
, 0, "verification", 1 },
2731 static struct argp argp
= { argp_options
, parse_opt
, 0, 0, argp_children
};
2732 struct barvinok_options
*bv_options
;
2734 bv_options
= barvinok_options_new_with_defaults();
2735 bv_options
->lookup_table
= 0;
2737 options
.verify
.barvinok
= bv_options
;
2738 set_program_name(argv
[0]);
2739 argp_parse(&argp
, argc
, argv
, 0, 0, &options
);
2742 C
= Constraints2Polyhedron(MA
, bv_options
->MaxRays
);
2744 fscanf(stdin
, " %d", &bignum
);
2745 assert(bignum
== -1);
2747 A
= Constraints2Polyhedron(MA
, bv_options
->MaxRays
);
2750 verify_options_set_range(&options
.verify
, A
->Dimension
);
2752 if (options
.verify
.verify
)
2755 iter_names
= util_generate_names(A
->Dimension
- C
->Dimension
, "i");
2756 param_names
= util_generate_names(C
->Dimension
, "p");
2757 if (print_solution
) {
2758 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
2759 Polyhedron_Print(stdout
, P_VALUE_FMT
, C
);
2761 vector
<max_term
*> maxima
= lexmin(A
, C
, &options
);
2763 for (int i
= 0; i
< maxima
.size(); ++i
)
2764 maxima
[i
]->print(cout
, param_names
, options
.verify
.barvinok
);
2766 if (options
.verify
.verify
)
2767 verify_results(A
, C
, maxima
, &options
.verify
);
2769 for (int i
= 0; i
< maxima
.size(); ++i
)
2772 util_free_names(A
->Dimension
- C
->Dimension
, iter_names
);
2773 util_free_names(C
->Dimension
, param_names
);
2777 barvinok_options_free(bv_options
);
2782 static bool lexmin(int pos
, Polyhedron
*P
, Value
*context
)
2791 value_init(LB
); value_init(UB
); value_init(k
);
2794 lu_flags
= lower_upper_bounds(pos
,P
,context
,&LB
,&UB
);
2795 assert(!(lu_flags
& LB_INFINITY
));
2797 value_set_si(context
[pos
],0);
2798 if (!lu_flags
&& value_lt(UB
,LB
)) {
2799 value_clear(LB
); value_clear(UB
); value_clear(k
);
2803 value_assign(context
[pos
], LB
);
2804 value_clear(LB
); value_clear(UB
); value_clear(k
);
2807 for (value_assign(k
,LB
); lu_flags
|| value_le(k
,UB
); value_increment(k
,k
)) {
2808 value_assign(context
[pos
],k
);
2809 if ((found
= lexmin(pos
+1, P
->next
, context
)))
2813 value_set_si(context
[pos
],0);
2814 value_clear(LB
); value_clear(UB
); value_clear(k
);
2818 static void print_list(FILE *out
, Value
*z
, char* brackets
, int len
)
2820 fprintf(out
, "%c", brackets
[0]);
2821 value_print(out
, VALUE_FMT
,z
[0]);
2822 for (int k
= 1; k
< len
; ++k
) {
2824 value_print(out
, VALUE_FMT
,z
[k
]);
2826 fprintf(out
, "%c", brackets
[1]);
2829 static int check_poly_lexmin(const struct check_poly_data
*data
,
2830 int nparam
, Value
*z
,
2831 const struct verify_options
*options
);
2833 struct check_poly_lexmin_data
: public check_poly_data
{
2835 vector
<max_term
*>& maxima
;
2837 check_poly_lexmin_data(Polyhedron
*S
, Value
*z
,
2838 vector
<max_term
*>& maxima
) : S(S
), maxima(maxima
) {
2840 this->check
= check_poly_lexmin
;
2844 static int check_poly_lexmin(const struct check_poly_data
*data
,
2845 int nparam
, Value
*z
,
2846 const struct verify_options
*options
)
2848 const check_poly_lexmin_data
*lexmin_data
;
2849 lexmin_data
= static_cast<const check_poly_lexmin_data
*>(data
);
2850 Polyhedron
*S
= lexmin_data
->S
;
2851 vector
<max_term
*>& maxima
= lexmin_data
->maxima
;
2853 bool found
= lexmin(1, S
, lexmin_data
->z
);
2855 if (options
->print_all
) {
2857 print_list(stdout
, z
, "()", nparam
);
2860 print_list(stdout
, lexmin_data
->z
+1, "[]", S
->Dimension
-nparam
);
2865 for (int i
= 0; i
< maxima
.size(); ++i
)
2866 if ((min
= maxima
[i
]->eval(z
, options
->barvinok
->MaxRays
)))
2869 int ok
= !(found
^ !!min
);
2871 for (int i
= 0; i
< S
->Dimension
-nparam
; ++i
)
2872 if (value_ne(lexmin_data
->z
[1+i
], min
->p
[i
])) {
2879 fprintf(stderr
, "Error !\n");
2880 fprintf(stderr
, "lexmin");
2881 print_list(stderr
, z
, "()", nparam
);
2882 fprintf(stderr
, " should be ");
2884 print_list(stderr
, lexmin_data
->z
+1, "[]", S
->Dimension
-nparam
);
2885 fprintf(stderr
, " while digging gives ");
2887 print_list(stderr
, min
->p
, "[]", S
->Dimension
-nparam
);
2888 fprintf(stderr
, ".\n");
2890 } else if (options
->print_all
)
2895 for (k
= 1; k
<= S
->Dimension
-nparam
; ++k
)
2896 value_set_si(lexmin_data
->z
[k
], 0);
2899 void verify_results(Polyhedron
*A
, Polyhedron
*C
, vector
<max_term
*>& maxima
,
2900 struct verify_options
*options
)
2903 unsigned nparam
= C
->Dimension
;
2904 unsigned MaxRays
= options
->barvinok
->MaxRays
;
2909 CS
= check_poly_context_scan(A
, &C
, nparam
, options
);
2911 p
= Vector_Alloc(A
->Dimension
+2);
2912 value_set_si(p
->p
[A
->Dimension
+1], 1);
2914 S
= Polyhedron_Scan(A
, C
, MaxRays
& POL_NO_DUAL
? 0 : MaxRays
);
2916 check_poly_init(C
, options
);
2919 if (!(CS
&& emptyQ2(CS
))) {
2920 check_poly_lexmin_data
data(S
, p
->p
, maxima
);
2921 check_poly(CS
, &data
, nparam
, 0, p
->p
+S
->Dimension
-nparam
+1, options
);
2926 if (!options
->print_all
)