9 #include <NTL/vec_ZZ.h>
10 #include <NTL/mat_ZZ.h>
12 #include <isl/space.h>
15 #include <isl/union_set.h>
16 #include <isl/union_map.h>
17 #include <isl/polynomial.h>
18 #include <isl_set_polylib.h>
19 #include <barvinok/polylib.h>
20 #include <barvinok/util.h>
21 #include <barvinok/evalue.h>
23 #include <barvinok/barvinok.h>
24 #include <barvinok/genfun.h>
25 #include <barvinok/options.h>
26 #include <barvinok/sample.h>
27 #include "bfcounter.h"
28 #include "conversion.h"
30 #include "decomposer.h"
32 #include "lattice_point.h"
34 #include "reduce_domain.h"
35 #include "remove_equalities.h"
38 #include "bernoulli.h"
39 #include "param_util.h"
49 using std::ostringstream
;
51 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
64 coeff
= Matrix_Alloc(d
+1, d
+1+1);
65 value_set_si(coeff
->p
[0][0], 1);
66 value_set_si(coeff
->p
[0][d
+1], 1);
67 for (int i
= 1; i
<= d
; ++i
) {
68 value_multiply(coeff
->p
[i
][0], coeff
->p
[i
-1][0], d0
);
69 Vector_Combine(coeff
->p
[i
-1], coeff
->p
[i
-1]+1, coeff
->p
[i
]+1,
71 value_set_si(coeff
->p
[i
][d
+1], i
);
72 value_multiply(coeff
->p
[i
][d
+1], coeff
->p
[i
][d
+1], coeff
->p
[i
-1][d
+1]);
73 value_decrement(d0
, d0
);
78 void div(dpoly
& d
, Vector
*count
, int sign
) {
79 int len
= coeff
->NbRows
;
80 Matrix
* c
= Matrix_Alloc(coeff
->NbRows
, coeff
->NbColumns
);
83 for (int i
= 0; i
< len
; ++i
) {
84 Vector_Copy(coeff
->p
[i
], c
->p
[i
], len
+1);
85 for (int j
= 1; j
<= i
; ++j
) {
86 value_multiply(tmp
, d
.coeff
->p
[j
], c
->p
[i
][len
]);
87 value_oppose(tmp
, tmp
);
88 Vector_Combine(c
->p
[i
], c
->p
[i
-j
], c
->p
[i
],
89 c
->p
[i
-j
][len
], tmp
, len
);
90 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], c
->p
[i
-j
][len
]);
92 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], d
.coeff
->p
[0]);
95 value_set_si(tmp
, -1);
96 Vector_Scale(c
->p
[len
-1], count
->p
, tmp
, len
);
97 value_assign(count
->p
[len
], c
->p
[len
-1][len
]);
99 Vector_Copy(c
->p
[len
-1], count
->p
, len
+1);
100 Vector_Normalize(count
->p
, len
+1);
106 struct bfe_term
: public bfc_term_base
{
107 vector
<evalue
*> factors
;
109 bfe_term(int len
) : bfc_term_base(len
) {
113 for (int i
= 0; i
< factors
.size(); ++i
) {
116 free_evalue_refs(factors
[i
]);
122 static void print_int_vector(int *v
, int len
, const char *name
)
124 cerr
<< name
<< endl
;
125 for (int j
= 0; j
< len
; ++j
) {
131 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
132 __attribute__((unused
));
133 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
136 cerr
<< "factors" << endl
;
137 cerr
<< factors
<< endl
;
138 for (int i
= 0; i
< v
.size(); ++i
) {
139 cerr
<< "term: " << i
<< endl
;
140 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
141 cerr
<< "terms" << endl
;
142 cerr
<< v
[i
]->terms
<< endl
;
143 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
144 cerr
<< bfct
->c
<< endl
;
148 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
149 __attribute__((unused
));
150 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
153 cerr
<< "factors" << endl
;
154 cerr
<< factors
<< endl
;
155 for (int i
= 0; i
< v
.size(); ++i
) {
156 cerr
<< "term: " << i
<< endl
;
157 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
158 cerr
<< "terms" << endl
;
159 cerr
<< v
[i
]->terms
<< endl
;
160 bfe_term
* bfet
= static_cast<bfe_term
*>(v
[i
]);
161 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
162 const char * test
[] = {"a", "b"};
163 print_evalue(stderr
, bfet
->factors
[j
], test
);
164 fprintf(stderr
, "\n");
169 struct bfcounter
: public bfcounter_base
{
173 bfcounter(unsigned dim
) : bfcounter_base(dim
) {
182 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
183 virtual void get_count(Value
*result
) {
184 assert(value_one_p(&count
[0]._mp_den
));
185 value_assign(*result
, &count
[0]._mp_num
);
189 void bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
191 unsigned nf
= factors
.NumRows();
193 for (int i
= 0; i
< v
.size(); ++i
) {
194 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
196 // factor is always positive, so we always
198 for (int k
= 0; k
< nf
; ++k
)
199 total_power
+= v
[i
]->powers
[k
];
202 for (j
= 0; j
< nf
; ++j
)
203 if (v
[i
]->powers
[j
] > 0)
206 zz2value(factors
[j
][0], tz
);
207 dpoly
D(total_power
, tz
, 1);
208 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
209 zz2value(factors
[j
][0], tz
);
210 dpoly
fact(total_power
, tz
, 1);
214 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
215 zz2value(factors
[j
][0], tz
);
216 dpoly
fact(total_power
, tz
, 1);
220 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
221 zz2value(v
[i
]->terms
[k
][0], tz
);
222 dpoly
n(total_power
, tz
);
223 mpq_set_si(tcount
, 0, 1);
226 bfct
->c
[k
].n
= -bfct
->c
[k
].n
;
227 zz2value(bfct
->c
[k
].n
, tn
);
228 zz2value(bfct
->c
[k
].d
, td
);
230 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
231 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
232 mpq_canonicalize(tcount
);
233 mpq_add(count
, count
, tcount
);
240 /* Check whether the polyhedron is unbounded and if so,
241 * check whether it has any (and therefore an infinite number of)
243 * If one of the vertices is integer, then we are done.
244 * Otherwise, transform the polyhedron such that one of the rays
245 * is the first unit vector and cut it off at a height that ensures
246 * that if the whole polyhedron has any points, then the remaining part
247 * has integer points. In particular we add the largest coefficient
248 * of a ray to the highest vertex (rounded up).
250 static bool Polyhedron_is_infinite(Polyhedron
*P
, Value
* result
,
251 barvinok_options
*options
)
263 for (; r
< P
->NbRays
; ++r
)
264 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
266 if (P
->NbBid
== 0 && r
== P
->NbRays
)
269 if (options
->count_sample_infinite
) {
272 sample
= Polyhedron_Sample(P
, options
);
274 value_set_si(*result
, 0);
276 value_set_si(*result
, -1);
282 for (int i
= 0; i
< P
->NbRays
; ++i
)
283 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
284 value_set_si(*result
, -1);
289 M
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
290 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
291 Vector_AntiScale(P
->Ray
[r
]+1, M
->p
[0], g
, P
->Dimension
+1);
292 int ok
= unimodular_complete(M
, 1);
294 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
297 P
= Polyhedron_Preimage(P
, M2
, 0);
305 value_set_si(size
, 0);
307 for (int i
= 0; i
< P
->NbBid
; ++i
) {
308 value_absolute(tmp
, P
->Ray
[i
][1]);
309 if (value_gt(tmp
, size
))
310 value_assign(size
, tmp
);
312 for (int i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
313 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
314 if (value_gt(P
->Ray
[i
][1], size
))
315 value_assign(size
, P
->Ray
[i
][1]);
318 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
319 if (first
|| value_gt(tmp
, offset
)) {
320 value_assign(offset
, tmp
);
324 value_addto(offset
, offset
, size
);
328 v
= Vector_Alloc(P
->Dimension
+2);
329 value_set_si(v
->p
[0], 1);
330 value_set_si(v
->p
[1], -1);
331 value_assign(v
->p
[1+P
->Dimension
], offset
);
332 R
= AddConstraints(v
->p
, 1, P
, options
->MaxRays
);
340 barvinok_count_with_options(P
, &c
, options
);
343 value_set_si(*result
, 0);
345 value_set_si(*result
, -1);
351 static void evalue2value(evalue
*e
, Value
*v
)
353 if (EVALUE_IS_ZERO(*e
)) {
358 if (value_notzero_p(e
->d
)) {
359 assert(value_one_p(e
->d
));
360 value_assign(*v
, e
->x
.n
);
364 assert(e
->x
.p
->type
== partition
);
365 assert(e
->x
.p
->size
== 2);
366 assert(EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
== 0);
367 evalue2value(&e
->x
.p
->arr
[1], v
);
370 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
371 barvinok_options
*options
);
373 void barvinok_count_with_options(Polyhedron
*P
, Value
* result
,
374 struct barvinok_options
*options
)
378 bool infinite
= false;
382 "barvinok_count: input is a union; only first polyhedron is counted\n");
385 value_set_si(*result
, 0);
391 P
= remove_equalities(P
, options
->MaxRays
);
393 P
= DomainConstraintSimplify(P
, options
->MaxRays
);
397 } while (P
&& !emptyQ(P
) && P
->NbEq
!= 0);
398 if (!P
|| emptyQ(P
)) {
400 value_set_si(*result
, 0);
405 if (Polyhedron_is_infinite(P
, result
, options
)) {
410 if (P
->Dimension
== 0) {
411 /* Test whether the constraints are satisfied */
412 POL_ENSURE_VERTICES(P
);
413 value_set_si(*result
, !emptyQ(P
));
418 if (options
->summation
== BV_SUM_BERNOULLI
) {
419 Polyhedron
*C
= Universe_Polyhedron(0);
420 evalue
*sum
= barvinok_summate_unweighted(P
, C
, options
);
422 evalue2value(sum
, result
);
426 Q
= Polyhedron_Factor(P
, 0, NULL
, options
->MaxRays
);
434 barvinok_count_f(P
, result
, options
);
435 if (value_neg_p(*result
))
437 if (Q
&& P
->next
&& value_notzero_p(*result
)) {
441 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
442 barvinok_count_f(Q
, &factor
, options
);
443 if (value_neg_p(factor
)) {
446 } else if (Q
->next
&& value_zero_p(factor
)) {
447 value_set_si(*result
, 0);
450 value_multiply(*result
, *result
, factor
);
459 value_set_si(*result
, -1);
462 void barvinok_count(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
464 barvinok_options
*options
= barvinok_options_new_with_defaults();
465 options
->MaxRays
= NbMaxCons
;
466 barvinok_count_with_options(P
, result
, options
);
467 barvinok_options_free(options
);
470 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
471 barvinok_options
*options
)
474 value_set_si(*result
, 0);
478 if (P
->Dimension
== 1)
479 return Line_Length(P
, result
);
481 int c
= P
->NbConstraints
;
482 POL_ENSURE_FACETS(P
);
483 if (c
!= P
->NbConstraints
|| P
->NbEq
!= 0) {
484 Polyhedron
*next
= P
->next
;
486 barvinok_count_with_options(P
, result
, options
);
491 POL_ENSURE_VERTICES(P
);
493 if (Polyhedron_is_infinite(P
, result
, options
))
497 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
498 cnt
= new bfcounter(P
->Dimension
);
499 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
500 cnt
= new icounter(P
->Dimension
);
501 else if (options
->incremental_specialization
== BV_SPECIALIZATION_TODD
)
502 cnt
= new tcounter(P
->Dimension
, options
->max_index
);
504 cnt
= new counter(P
->Dimension
, options
->max_index
);
505 cnt
->start(P
, options
);
507 cnt
->get_count(result
);
511 typedef evalue
* evalue_p
;
513 struct enumerator_base
{
517 vertex_decomposer
*vpd
;
519 enumerator_base(unsigned dim
, vertex_decomposer
*vpd
)
524 vE
= new evalue_p
[vpd
->PP
->nbV
];
525 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
529 evalue_set_si(&mone
, -1, 1);
532 void decompose_at(Param_Vertices
*V
, int _i
, barvinok_options
*options
) {
536 value_init(vE
[_i
]->d
);
537 evalue_set_si(vE
[_i
], 0, 1);
539 vpd
->decompose_at_vertex(V
, _i
, options
);
542 virtual ~enumerator_base() {
543 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
545 free_evalue_refs(vE
[j
]);
550 free_evalue_refs(&mone
);
553 static enumerator_base
*create(Polyhedron
*P
, unsigned dim
,
554 Param_Polyhedron
*PP
,
555 barvinok_options
*options
);
558 struct enumerator
: public signed_cone_consumer
, public vertex_decomposer
,
559 public enumerator_base
{
567 enumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
568 vertex_decomposer(PP
, *this), enumerator_base(dim
, this) {
569 randomvector(P
, lambda
, dim
, 0);
571 c
= Vector_Alloc(dim
+2);
583 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
586 void enumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
589 assert(sc
.rays
.NumRows() == dim
);
590 for (int k
= 0; k
< dim
; ++k
) {
591 if (lambda
* sc
.rays
[k
] == 0)
595 lattice_point(V
, sc
.rays
, lambda
, &num
, sc
.det
, options
);
596 den
= sc
.rays
* lambda
;
601 zz2value(den
[0], tz
);
603 for (int k
= 1; k
< dim
; ++k
) {
604 zz2value(den
[k
], tz
);
605 dpoly
fact(dim
, tz
, 1);
611 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
612 evalue
*EV
= evalue_polynomial(c
, num
.E
[i
]);
615 evalue_free(num
.E
[i
]);
619 mpq_set_si(count
, 0, 1);
620 if (num
.constant
.length() == 1) {
621 zz2value(num
.constant
[0], tz
);
623 d
.div(n
, count
, sign
);
630 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
631 value_assign(acc
, c
->p
[dim
]);
632 zz2value(num
.constant
[i
], x
);
633 for (int j
= dim
-1; j
>= 0; --j
) {
634 value_multiply(acc
, acc
, x
);
635 value_addto(acc
, acc
, c
->p
[j
]);
637 value_addto(mpq_numref(count
), mpq_numref(count
), acc
);
639 mpz_set(mpq_denref(count
), c
->p
[dim
+1]);
645 evalue_set(&EV
, &count
[0]._mp_num
, &count
[0]._mp_den
);
647 free_evalue_refs(&EV
);
651 struct ienumerator_base
: enumerator_base
{
654 ienumerator_base(unsigned dim
, vertex_decomposer
*vpd
) :
655 enumerator_base(dim
,vpd
) {
656 E_vertex
= new evalue_p
[dim
];
659 virtual ~ienumerator_base() {
663 evalue
*E_num(int i
, int d
) {
664 return E_vertex
[i
+ (dim
-d
)];
673 cumulator(evalue
*factor
, evalue
*v
, dpoly_r
*r
) :
674 factor(factor
), v(v
), r(r
) {}
676 void cumulate(barvinok_options
*options
);
678 virtual void add_term(const vector
<int>& powers
, evalue
*f2
) = 0;
679 virtual ~cumulator() {}
682 void cumulator::cumulate(barvinok_options
*options
)
684 evalue cum
; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
686 evalue t
; // E_num[0] - (m-1)
690 if (options
->lookup_table
) {
692 evalue_set_si(&mone
, -1, 1);
696 evalue_copy(&cum
, factor
);
699 value_set_si(f
.d
, 1);
700 value_set_si(f
.x
.n
, 1);
704 if (!options
->lookup_table
) {
705 for (cst
= &t
; value_zero_p(cst
->d
); ) {
706 if (cst
->x
.p
->type
== fractional
)
707 cst
= &cst
->x
.p
->arr
[1];
709 cst
= &cst
->x
.p
->arr
[0];
713 for (int m
= 0; m
< r
->len
; ++m
) {
716 value_set_si(f
.d
, m
);
718 if (!options
->lookup_table
)
719 value_subtract(cst
->x
.n
, cst
->x
.n
, cst
->d
);
725 dpoly_r_term_list
& current
= r
->c
[r
->len
-1-m
];
726 dpoly_r_term_list::iterator j
;
727 for (j
= current
.begin(); j
!= current
.end(); ++j
) {
728 if ((*j
)->coeff
== 0)
730 evalue
*f2
= new evalue
;
733 zz2value((*j
)->coeff
, f2
->x
.n
);
734 zz2value(r
->denom
, f2
->d
);
737 add_term((*j
)->powers
, f2
);
740 free_evalue_refs(&f
);
741 free_evalue_refs(&t
);
742 free_evalue_refs(&cum
);
743 if (options
->lookup_table
)
744 free_evalue_refs(&mone
);
752 struct ie_cum
: public cumulator
{
753 vector
<E_poly_term
*> terms
;
755 ie_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
) : cumulator(factor
, v
, r
) {}
757 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
760 void ie_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
763 for (k
= 0; k
< terms
.size(); ++k
) {
764 if (terms
[k
]->powers
== powers
) {
765 eadd(f2
, terms
[k
]->E
);
766 free_evalue_refs(f2
);
771 if (k
>= terms
.size()) {
772 E_poly_term
*ET
= new E_poly_term
;
779 struct ienumerator
: public signed_cone_consumer
, public vertex_decomposer
,
780 public ienumerator_base
{
786 ienumerator(unsigned dim
, Param_Polyhedron
*PP
) :
787 vertex_decomposer(PP
, *this), ienumerator_base(dim
, this) {
788 vertex
.SetDims(1, dim
);
790 den
.SetDims(dim
, dim
);
800 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
801 void reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
802 barvinok_options
*options
);
805 void ienumerator::reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
806 barvinok_options
*options
)
808 unsigned len
= den_f
.NumRows(); // number of factors in den
809 unsigned dim
= num
.NumCols();
810 assert(num
.NumRows() == 1);
813 eadd(factor
, vE
[vert
]);
822 split_one(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
825 den_p
.SetLength(len
);
829 normalize(one
, num_s
, num_p
, den_s
, den_p
, den_r
);
835 for (int k
= 0; k
< len
; ++k
) {
838 else if (den_s
[k
] == 0)
842 reduce(factor
, num_p
, den_r
, options
);
846 pden
.SetDims(only_param
, dim
-1);
848 for (k
= 0, l
= 0; k
< len
; ++k
)
850 pden
[l
++] = den_r
[k
];
852 for (k
= 0; k
< len
; ++k
)
856 zz2value(num_s
[0], tz
);
857 dpoly
n(no_param
, tz
);
858 zz2value(den_s
[k
], tz
);
859 dpoly
D(no_param
, tz
, 1);
862 zz2value(den_s
[k
], tz
);
863 dpoly
fact(no_param
, tz
, 1);
868 // if no_param + only_param == len then all powers
869 // below will be all zero
870 if (no_param
+ only_param
== len
) {
871 if (E_num(0, dim
) != 0)
872 r
= new dpoly_r(n
, len
);
874 mpq_set_si(tcount
, 0, 1);
878 if (value_notzero_p(mpq_numref(tcount
))) {
882 value_assign(f
.x
.n
, mpq_numref(tcount
));
883 value_assign(f
.d
, mpq_denref(tcount
));
885 reduce(factor
, num_p
, pden
, options
);
886 free_evalue_refs(&f
);
891 for (k
= 0; k
< len
; ++k
) {
892 if (den_s
[k
] == 0 || den_p
[k
] == 0)
895 zz2value(den_s
[k
], tz
);
896 dpoly
pd(no_param
-1, tz
, 1);
899 for (l
= 0; l
< k
; ++l
)
900 if (den_r
[l
] == den_r
[k
])
904 r
= new dpoly_r(n
, pd
, l
, len
);
906 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
912 dpoly_r
*rc
= r
->div(D
);
915 if (E_num(0, dim
) == 0) {
916 int common
= pden
.NumRows();
917 dpoly_r_term_list
& final
= r
->c
[r
->len
-1];
923 zz2value(r
->denom
, f
.d
);
924 dpoly_r_term_list::iterator j
;
925 for (j
= final
.begin(); j
!= final
.end(); ++j
) {
926 if ((*j
)->coeff
== 0)
929 for (int k
= 0; k
< r
->dim
; ++k
) {
930 int n
= (*j
)->powers
[k
];
933 pden
.SetDims(rows
+n
, pden
.NumCols());
934 for (int l
= 0; l
< n
; ++l
)
935 pden
[rows
+l
] = den_r
[k
];
939 evalue_copy(&t
, factor
);
940 zz2value((*j
)->coeff
, f
.x
.n
);
942 reduce(&t
, num_p
, pden
, options
);
943 free_evalue_refs(&t
);
945 free_evalue_refs(&f
);
947 ie_cum
cum(factor
, E_num(0, dim
), r
);
948 cum
.cumulate(options
);
950 int common
= pden
.NumRows();
952 for (int j
= 0; j
< cum
.terms
.size(); ++j
) {
954 pden
.SetDims(rows
, pden
.NumCols());
955 for (int k
= 0; k
< r
->dim
; ++k
) {
956 int n
= cum
.terms
[j
]->powers
[k
];
959 pden
.SetDims(rows
+n
, pden
.NumCols());
960 for (int l
= 0; l
< n
; ++l
)
961 pden
[rows
+l
] = den_r
[k
];
964 reduce(cum
.terms
[j
]->E
, num_p
, pden
, options
);
965 free_evalue_refs(cum
.terms
[j
]->E
);
966 delete cum
.terms
[j
]->E
;
974 void ienumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
977 assert(sc
.rays
.NumRows() == dim
);
979 lattice_point(V
, sc
.rays
, vertex
[0], E_vertex
, options
);
985 evalue_set_si(&one
, sc
.sign
, 1);
986 reduce(&one
, vertex
, den
, options
);
987 free_evalue_refs(&one
);
989 for (int i
= 0; i
< dim
; ++i
)
991 evalue_free(E_vertex
[i
]);
994 struct bfenumerator
: public vertex_decomposer
, public bf_base
,
995 public ienumerator_base
{
998 bfenumerator(unsigned dim
, Param_Polyhedron
*PP
) :
999 vertex_decomposer(PP
, *this),
1000 bf_base(dim
), ienumerator_base(dim
, this) {
1008 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
1009 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1011 bfc_term_base
* new_bf_term(int len
) {
1012 bfe_term
* t
= new bfe_term(len
);
1016 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
1017 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1018 factor
= bfet
->factors
[k
];
1019 assert(factor
!= NULL
);
1020 bfet
->factors
[k
] = NULL
;
1022 emul(&mone
, factor
);
1025 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&q
, int change
) {
1026 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1027 factor
= bfet
->factors
[k
];
1028 assert(factor
!= NULL
);
1029 bfet
->factors
[k
] = NULL
;
1035 value_oppose(f
.x
.n
, mpq_numref(q
));
1037 value_assign(f
.x
.n
, mpq_numref(q
));
1038 value_assign(f
.d
, mpq_denref(q
));
1040 free_evalue_refs(&f
);
1043 virtual void set_factor(bfc_term_base
*t
, int k
, const QQ
& c
, int change
) {
1044 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1046 factor
= new evalue
;
1051 zz2value(c
.n
, f
.x
.n
);
1053 value_oppose(f
.x
.n
, f
.x
.n
);
1056 value_init(factor
->d
);
1057 evalue_copy(factor
, bfet
->factors
[k
]);
1059 free_evalue_refs(&f
);
1062 void set_factor(evalue
*f
, int change
) {
1068 virtual void insert_term(bfc_term_base
*t
, int i
) {
1069 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1070 int len
= t
->terms
.NumRows()-1; // already increased by one
1072 bfet
->factors
.resize(len
+1);
1073 for (int j
= len
; j
> i
; --j
) {
1074 bfet
->factors
[j
] = bfet
->factors
[j
-1];
1075 t
->terms
[j
] = t
->terms
[j
-1];
1077 bfet
->factors
[i
] = factor
;
1081 virtual void update_term(bfc_term_base
*t
, int i
) {
1082 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1084 eadd(factor
, bfet
->factors
[i
]);
1085 free_evalue_refs(factor
);
1089 virtual bool constant_vertex(int dim
) { return E_num(0, dim
) == 0; }
1091 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
,
1092 barvinok_options
*options
);
1095 enumerator_base
*enumerator_base::create(Polyhedron
*P
, unsigned dim
,
1096 Param_Polyhedron
*PP
,
1097 barvinok_options
*options
)
1099 enumerator_base
*eb
;
1101 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
1102 eb
= new bfenumerator(dim
, PP
);
1103 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
1104 eb
= new ienumerator(dim
, PP
);
1106 eb
= new enumerator(P
, dim
, PP
);
1111 struct bfe_cum
: public cumulator
{
1113 bfc_term_base
*told
;
1117 bfe_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
, bf_reducer
*bfr
,
1118 bfc_term_base
*t
, int k
, bfenumerator
*e
) :
1119 cumulator(factor
, v
, r
), bfe(e
), told(t
), k(k
),
1123 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
1126 void bfe_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
1128 bfr
->update_powers(powers
);
1130 bfc_term_base
* t
= bfe
->find_bfc_term(bfr
->vn
, bfr
->npowers
, bfr
->nnf
);
1131 bfe
->set_factor(f2
, bfr
->l_changes
% 2);
1132 bfe
->add_term(t
, told
->terms
[k
], bfr
->l_extra_num
);
1135 void bfenumerator::cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
1136 dpoly_r
*r
, barvinok_options
*options
)
1138 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1139 bfe_cum
cum(bfet
->factors
[k
], E_num(0, bfr
->d
), r
, bfr
, t
, k
, this);
1140 cum
.cumulate(options
);
1143 void bfenumerator::base(mat_ZZ
& factors
, bfc_vec
& v
)
1145 for (int i
= 0; i
< v
.size(); ++i
) {
1146 assert(v
[i
]->terms
.NumRows() == 1);
1147 evalue
*factor
= static_cast<bfe_term
*>(v
[i
])->factors
[0];
1148 eadd(factor
, vE
[vert
]);
1153 void bfenumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1155 assert(sc
.det
== 1);
1156 assert(sc
.rays
.NumRows() == enumerator_base::dim
);
1158 bfe_term
* t
= new bfe_term(enumerator_base::dim
);
1159 vector
< bfc_term_base
* > v
;
1162 t
->factors
.resize(1);
1164 t
->terms
.SetDims(1, enumerator_base::dim
);
1165 lattice_point(V
, sc
.rays
, t
->terms
[0], E_vertex
, options
);
1167 // the elements of factors are always lexpositive
1169 int s
= setup_factors(sc
.rays
, factors
, t
, sc
.sign
);
1171 t
->factors
[0] = new evalue
;
1172 value_init(t
->factors
[0]->d
);
1173 evalue_set_si(t
->factors
[0], s
, 1);
1174 reduce(factors
, v
, options
);
1176 for (int i
= 0; i
< enumerator_base::dim
; ++i
)
1178 evalue_free(E_vertex
[i
]);
1181 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1182 barvinok_options
*options
);
1185 static evalue
* barvinok_enumerate_cst(Polyhedron
*P
, Polyhedron
* C
,
1186 struct barvinok_options
*options
)
1192 return evalue_zero();
1195 ALLOC(evalue
, eres
);
1196 value_init(eres
->d
);
1197 value_set_si(eres
->d
, 0);
1198 eres
->x
.p
= new_enode(partition
, 2, C
->Dimension
);
1199 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[0],
1200 DomainConstraintSimplify(C
, options
->MaxRays
));
1201 value_set_si(eres
->x
.p
->arr
[1].d
, 1);
1202 value_init(eres
->x
.p
->arr
[1].x
.n
);
1204 value_set_si(eres
->x
.p
->arr
[1].x
.n
, 0);
1206 barvinok_count_with_options(P
, &eres
->x
.p
->arr
[1].x
.n
, options
);
1207 if (value_mone_p(eres
->x
.p
->arr
[1].x
.n
)) {
1208 value_clear(eres
->x
.p
->arr
[1].x
.n
);
1209 value_set_si(eres
->x
.p
->arr
[1].d
, -2); /* NaN */
1215 static evalue
* enumerate(Polyhedron
*P
, Polyhedron
* C
,
1216 struct barvinok_options
*options
)
1219 Polyhedron
*Porig
= P
;
1220 Polyhedron
*Corig
= C
;
1221 Polyhedron
*CEq
= NULL
;
1222 unsigned nparam
= C
->Dimension
;
1227 value_init(factor
.d
);
1228 evalue_set_si(&factor
, 1, 1);
1231 POL_ENSURE_FACETS(P
);
1232 POL_ENSURE_VERTICES(P
);
1233 POL_ENSURE_FACETS(C
);
1234 POL_ENSURE_VERTICES(C
);
1236 if (C
->Dimension
== 0 || emptyQ(P
) || emptyQ(C
)) {
1239 CEq
= Polyhedron_Copy(CEq
);
1240 eres
= barvinok_enumerate_cst(P
, CEq
? CEq
: Polyhedron_Copy(C
), options
);
1243 evalue_backsubstitute(eres
, CP
, options
->MaxRays
);
1247 emul(&factor
, eres
);
1248 if (options
->approx
->method
== BV_APPROX_DROP
) {
1249 if (options
->approx
->approximation
== BV_APPROX_SIGN_UPPER
)
1250 evalue_frac2polynomial(eres
, 1, options
->MaxRays
);
1251 if (options
->approx
->approximation
== BV_APPROX_SIGN_LOWER
)
1252 evalue_frac2polynomial(eres
, -1, options
->MaxRays
);
1253 if (options
->approx
->approximation
== BV_APPROX_SIGN_APPROX
)
1254 evalue_frac2polynomial(eres
, 0, options
->MaxRays
);
1256 reduce_evalue(eres
);
1257 free_evalue_refs(&factor
);
1265 if (Polyhedron_is_unbounded(P
, nparam
, options
->MaxRays
))
1268 if (P
->Dimension
== nparam
) {
1269 CEq
= DomainIntersection(P
, C
, options
->MaxRays
);
1270 P
= Universe_Polyhedron(0);
1273 if (P
->NbEq
!= 0 || C
->NbEq
!= 0) {
1276 remove_all_equalities(&P
, &C
, &CP
, NULL
, nparam
, options
->MaxRays
);
1277 if (C
!= D
&& D
!= Corig
)
1279 if (P
!= Q
&& Q
!= Porig
)
1281 eres
= enumerate(P
, C
, options
);
1285 Polyhedron
*T
= Polyhedron_Factor(P
, nparam
, NULL
, options
->MaxRays
);
1286 if (T
|| (P
->Dimension
== nparam
+1)) {
1288 Polyhedron
*FC
= Factor_Context(T
? T
: P
, nparam
, options
->MaxRays
);
1289 C
= DomainIntersection(C
, FC
, options
->MaxRays
);
1291 Polyhedron_Free(C2
);
1292 Polyhedron_Free(FC
);
1298 if (T
->Dimension
== C
->Dimension
) {
1307 eres
= barvinok_enumerate_ev_f(P
, C
, options
);
1314 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
1315 Polyhedron
*next
= Q
->next
;
1318 f
= barvinok_enumerate_ev_f(Q
, C
, options
);
1329 evalue
* barvinok_enumerate_with_options(Polyhedron
*P
, Polyhedron
* C
,
1330 struct barvinok_options
*options
)
1332 Polyhedron
*next
, *Cnext
, *C1
;
1333 Polyhedron
*Corig
= C
;
1338 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1342 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1346 C1
= Polyhedron_Project(P
, C
->Dimension
);
1347 C
= DomainIntersection(C
, C1
, options
->MaxRays
);
1348 Polyhedron_Free(C1
);
1352 if (options
->approx
->method
== BV_APPROX_BERNOULLI
||
1353 options
->summation
== BV_SUM_BERNOULLI
) {
1354 int summation
= options
->summation
;
1355 options
->summation
= BV_SUM_BERNOULLI
;
1356 eres
= barvinok_summate_unweighted(P
, C
, options
);
1357 options
->summation
= summation
;
1359 eres
= enumerate(P
, C
, options
);
1363 Corig
->next
= Cnext
;
1368 evalue
* barvinok_enumerate_ev(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1371 barvinok_options
*options
= barvinok_options_new_with_defaults();
1372 options
->MaxRays
= MaxRays
;
1373 E
= barvinok_enumerate_with_options(P
, C
, options
);
1374 barvinok_options_free(options
);
1378 evalue
*Param_Polyhedron_Enumerate(Param_Polyhedron
*PP
, Polyhedron
*P
,
1380 struct barvinok_options
*options
)
1384 unsigned nparam
= C
->Dimension
;
1385 unsigned dim
= P
->Dimension
- nparam
;
1388 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
1389 evalue_section
*s
= new evalue_section
[nd
];
1390 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
1392 enumerator_base
*et
= NULL
;
1397 et
= enumerator_base::create(P
, dim
, PP
, options
);
1399 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
1402 s
[i
].E
= evalue_zero();
1405 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
1408 et
->decompose_at(V
, _i
, options
);
1409 } catch (OrthogonalException
&e
) {
1410 FORALL_REDUCED_DOMAIN_RESET
;
1411 for (; i
>= 0; --i
) {
1412 evalue_free(s
[i
].E
);
1413 Domain_Free(s
[i
].D
);
1417 eadd(et
->vE
[_i
] , s
[i
].E
);
1418 END_FORALL_PVertex_in_ParamPolyhedron
;
1419 evalue_range_reduction_in_domain(s
[i
].E
, rVD
);
1420 END_FORALL_REDUCED_DOMAIN
1421 Polyhedron_Free(TC
);
1424 eres
= evalue_from_section_array(s
, nd
);
1430 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1431 barvinok_options
*options
)
1433 unsigned nparam
= C
->Dimension
;
1434 bool do_scale
= options
->approx
->method
== BV_APPROX_SCALE
;
1436 if (options
->summation
== BV_SUM_EULER
)
1437 return barvinok_summate_unweighted(P
, C
, options
);
1439 if (options
->approx
->method
== BV_APPROX_VOLUME
)
1440 return Param_Polyhedron_Volume(P
, C
, options
);
1442 if (P
->Dimension
- nparam
== 1 && !do_scale
)
1443 return ParamLine_Length(P
, C
, options
);
1445 Param_Polyhedron
*PP
= NULL
;
1449 eres
= scale_bound(P
, C
, options
);
1454 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
1457 eres
= scale(PP
, P
, C
, options
);
1459 eres
= Param_Polyhedron_Enumerate(PP
, P
, C
, options
);
1462 Param_Polyhedron_Free(PP
);
1467 Enumeration
* barvinok_enumerate(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1469 evalue
*EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
1471 return partition2enumeration(EP
);
1474 evalue
* barvinok_enumerate_union(Polyhedron
*D
, Polyhedron
* C
, unsigned MaxRays
)
1477 gen_fun
*gf
= barvinok_enumerate_union_series(D
, C
, MaxRays
);
1483 static __isl_give isl_pw_qpolynomial
*basic_set_card(
1484 __isl_take isl_basic_set
*bset
)
1488 isl_pw_qpolynomial
*pwqp
;
1489 unsigned nparam
= isl_basic_set_dim(bset
, isl_dim_param
);
1490 Polyhedron
*U
= Universe_Polyhedron(nparam
);
1493 barvinok_options
*options
;
1494 int options_allocated
= 0;
1496 ctx
= isl_basic_set_get_ctx(bset
);
1497 options
= isl_ctx_peek_barvinok_options(ctx
);
1499 options
= barvinok_options_new_with_defaults();
1500 options_allocated
= 1;
1503 dim
= isl_basic_set_get_space(bset
);
1504 dim
= isl_space_domain(dim
);
1506 P
= isl_basic_set_to_polylib(bset
);
1507 E
= enumerate(P
, U
, options
);
1509 pwqp
= isl_pw_qpolynomial_from_evalue(dim
, E
);
1510 isl_basic_set_free(bset
);
1515 if (options_allocated
)
1516 barvinok_options_free(options
);
1521 static isl_stat
basic_map_card(__isl_take isl_basic_map
*bmap
, void *user
)
1523 isl_pw_qpolynomial
**sum
= (isl_pw_qpolynomial
**)user
;
1524 isl_pw_qpolynomial
*pwqp
;
1525 unsigned nparam
= isl_basic_map_dim(bmap
, isl_dim_param
);
1526 unsigned n_in
= isl_basic_map_dim(bmap
, isl_dim_in
);
1527 isl_space
*target_dim
;
1528 isl_basic_set
*bset
;
1530 target_dim
= isl_basic_map_get_space(bmap
);
1531 target_dim
= isl_space_domain(target_dim
);
1533 bmap
= isl_basic_map_move_dims(bmap
, isl_dim_param
, nparam
,
1534 isl_dim_in
, 0, n_in
);
1536 bset
= isl_basic_map_range(bmap
);
1537 bset
= isl_basic_set_lift(bset
);
1538 pwqp
= isl_basic_set_multiplicative_call(bset
, &basic_set_card
);
1540 pwqp
= isl_pw_qpolynomial_move_dims(pwqp
, isl_dim_in
, 0,
1541 isl_dim_param
, nparam
, n_in
);
1542 pwqp
= isl_pw_qpolynomial_reset_domain_space(pwqp
, target_dim
);
1543 *sum
= isl_pw_qpolynomial_add(*sum
, pwqp
);
1548 static __isl_give isl_pw_qpolynomial
*card_as_sum(__isl_take isl_map
*map
,
1549 barvinok_options
*options
)
1555 isl_qpolynomial
*qp
;
1556 isl_pw_qpolynomial
*pwqp
;
1557 int summation
= options
->summation
;
1562 options
->summation
= BV_SUM_BERNOULLI
;
1564 set
= isl_map_wrap(map
);
1565 dim
= isl_set_get_space(set
);
1566 ctx
= isl_map_get_ctx(map
);
1567 one
= isl_val_one(ctx
);
1568 qp
= isl_qpolynomial_val_on_domain(dim
, one
);
1570 pwqp
= isl_pw_qpolynomial_alloc(set
, qp
);
1571 pwqp
= isl_pw_qpolynomial_sum(pwqp
);
1573 options
->summation
= summation
;
1578 __isl_give isl_pw_qpolynomial
*isl_map_card(__isl_take isl_map
*map
)
1582 isl_pw_qpolynomial
*sum
;
1583 barvinok_options
*options
;
1585 ctx
= isl_map_get_ctx(map
);
1586 options
= isl_ctx_peek_barvinok_options(ctx
);
1588 (options
->approx
->method
== BV_APPROX_BERNOULLI
||
1589 options
->summation
== BV_SUM_BERNOULLI
))
1590 return card_as_sum(map
, options
);
1592 dim
= isl_map_get_space(map
);
1593 dim
= isl_space_domain(dim
);
1594 dim
= isl_space_from_domain(dim
);
1595 dim
= isl_space_add_dims(dim
, isl_dim_out
, 1);
1596 sum
= isl_pw_qpolynomial_zero(dim
);
1598 map
= isl_map_make_disjoint(map
);
1599 map
= isl_map_compute_divs(map
);
1601 if (isl_map_foreach_basic_map(map
, &basic_map_card
, &sum
) < 0)
1609 isl_pw_qpolynomial_free(sum
);
1613 __isl_give isl_pw_qpolynomial
*isl_set_card(__isl_take isl_set
*set
)
1615 isl_pw_qpolynomial
*pwqp
;
1616 pwqp
= isl_map_card(isl_map_from_range(set
));
1617 pwqp
= isl_pw_qpolynomial_project_domain_on_params(pwqp
);
1621 __isl_give isl_pw_qpolynomial
*isl_basic_map_card(__isl_take isl_basic_map
*bmap
)
1623 return isl_map_card(isl_map_from_basic_map(bmap
));
1626 __isl_give isl_pw_qpolynomial
*isl_basic_set_card(__isl_take isl_basic_set
*bset
)
1628 isl_pw_qpolynomial
*pwqp
;
1629 pwqp
= isl_basic_map_card(isl_basic_map_from_range(bset
));
1630 pwqp
= isl_pw_qpolynomial_project_domain_on_params(pwqp
);
1634 static isl_stat
set_card(__isl_take isl_set
*set
, void *user
)
1636 isl_union_pw_qpolynomial
**res
= (isl_union_pw_qpolynomial
**)user
;
1637 isl_pw_qpolynomial
*pwqp
;
1638 isl_union_pw_qpolynomial
*upwqp
;
1640 pwqp
= isl_set_card(set
);
1641 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1642 *res
= isl_union_pw_qpolynomial_add(*res
, upwqp
);
1647 __isl_give isl_union_pw_qpolynomial
*isl_union_set_card(
1648 __isl_take isl_union_set
*uset
)
1651 isl_union_pw_qpolynomial
*res
;
1653 dim
= isl_union_set_get_space(uset
);
1654 res
= isl_union_pw_qpolynomial_zero(dim
);
1655 if (isl_union_set_foreach_set(uset
, &set_card
, &res
) < 0)
1657 isl_union_set_free(uset
);
1661 isl_union_set_free(uset
);
1662 isl_union_pw_qpolynomial_free(res
);
1666 static isl_stat
map_card(__isl_take isl_map
*map
, void *user
)
1668 isl_union_pw_qpolynomial
**res
= (isl_union_pw_qpolynomial
**)user
;
1669 isl_pw_qpolynomial
*pwqp
;
1670 isl_union_pw_qpolynomial
*upwqp
;
1672 pwqp
= isl_map_card(map
);
1673 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1674 *res
= isl_union_pw_qpolynomial_add(*res
, upwqp
);
1679 __isl_give isl_union_pw_qpolynomial
*isl_union_map_card(
1680 __isl_take isl_union_map
*umap
)
1683 isl_union_pw_qpolynomial
*res
;
1685 dim
= isl_union_map_get_space(umap
);
1686 res
= isl_union_pw_qpolynomial_zero(dim
);
1687 if (isl_union_map_foreach_map(umap
, &map_card
, &res
) < 0)
1689 isl_union_map_free(umap
);
1693 isl_union_map_free(umap
);
1694 isl_union_pw_qpolynomial_free(res
);