8 #include <NTL/mat_ZZ.h>
11 #include <isl_set_polylib.h>
12 #include <barvinok/util.h>
13 #include <barvinok/evalue.h>
15 #include <barvinok/barvinok.h>
16 #include <barvinok/genfun.h>
17 #include <barvinok/options.h>
18 #include <barvinok/sample.h>
19 #include "bfcounter.h"
20 #include "conversion.h"
22 #include "decomposer.h"
24 #include "lattice_point.h"
26 #include "reduce_domain.h"
27 #include "remove_equalities.h"
30 #include "bernoulli.h"
31 #include "param_util.h"
41 using std::ostringstream
;
43 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
56 coeff
= Matrix_Alloc(d
+1, d
+1+1);
57 value_set_si(coeff
->p
[0][0], 1);
58 value_set_si(coeff
->p
[0][d
+1], 1);
59 for (int i
= 1; i
<= d
; ++i
) {
60 value_multiply(coeff
->p
[i
][0], coeff
->p
[i
-1][0], d0
);
61 Vector_Combine(coeff
->p
[i
-1], coeff
->p
[i
-1]+1, coeff
->p
[i
]+1,
63 value_set_si(coeff
->p
[i
][d
+1], i
);
64 value_multiply(coeff
->p
[i
][d
+1], coeff
->p
[i
][d
+1], coeff
->p
[i
-1][d
+1]);
65 value_decrement(d0
, d0
);
70 void div(dpoly
& d
, Vector
*count
, int sign
) {
71 int len
= coeff
->NbRows
;
72 Matrix
* c
= Matrix_Alloc(coeff
->NbRows
, coeff
->NbColumns
);
75 for (int i
= 0; i
< len
; ++i
) {
76 Vector_Copy(coeff
->p
[i
], c
->p
[i
], len
+1);
77 for (int j
= 1; j
<= i
; ++j
) {
78 value_multiply(tmp
, d
.coeff
->p
[j
], c
->p
[i
][len
]);
79 value_oppose(tmp
, tmp
);
80 Vector_Combine(c
->p
[i
], c
->p
[i
-j
], c
->p
[i
],
81 c
->p
[i
-j
][len
], tmp
, len
);
82 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], c
->p
[i
-j
][len
]);
84 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], d
.coeff
->p
[0]);
87 value_set_si(tmp
, -1);
88 Vector_Scale(c
->p
[len
-1], count
->p
, tmp
, len
);
89 value_assign(count
->p
[len
], c
->p
[len
-1][len
]);
91 Vector_Copy(c
->p
[len
-1], count
->p
, len
+1);
92 Vector_Normalize(count
->p
, len
+1);
98 struct bfe_term
: public bfc_term_base
{
99 vector
<evalue
*> factors
;
101 bfe_term(int len
) : bfc_term_base(len
) {
105 for (int i
= 0; i
< factors
.size(); ++i
) {
108 free_evalue_refs(factors
[i
]);
114 static void print_int_vector(int *v
, int len
, const char *name
)
116 cerr
<< name
<< endl
;
117 for (int j
= 0; j
< len
; ++j
) {
123 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
126 cerr
<< "factors" << endl
;
127 cerr
<< factors
<< endl
;
128 for (int i
= 0; i
< v
.size(); ++i
) {
129 cerr
<< "term: " << i
<< endl
;
130 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
131 cerr
<< "terms" << endl
;
132 cerr
<< v
[i
]->terms
<< endl
;
133 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
134 cerr
<< bfct
->c
<< endl
;
138 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
141 cerr
<< "factors" << endl
;
142 cerr
<< factors
<< endl
;
143 for (int i
= 0; i
< v
.size(); ++i
) {
144 cerr
<< "term: " << i
<< endl
;
145 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
146 cerr
<< "terms" << endl
;
147 cerr
<< v
[i
]->terms
<< endl
;
148 bfe_term
* bfet
= static_cast<bfe_term
*>(v
[i
]);
149 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
150 const char * test
[] = {"a", "b"};
151 print_evalue(stderr
, bfet
->factors
[j
], test
);
152 fprintf(stderr
, "\n");
157 struct bfcounter
: public bfcounter_base
{
161 bfcounter(unsigned dim
) : bfcounter_base(dim
) {
170 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
171 virtual void get_count(Value
*result
) {
172 assert(value_one_p(&count
[0]._mp_den
));
173 value_assign(*result
, &count
[0]._mp_num
);
177 void bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
179 unsigned nf
= factors
.NumRows();
181 for (int i
= 0; i
< v
.size(); ++i
) {
182 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
184 // factor is always positive, so we always
186 for (int k
= 0; k
< nf
; ++k
)
187 total_power
+= v
[i
]->powers
[k
];
190 for (j
= 0; j
< nf
; ++j
)
191 if (v
[i
]->powers
[j
] > 0)
194 zz2value(factors
[j
][0], tz
);
195 dpoly
D(total_power
, tz
, 1);
196 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
197 zz2value(factors
[j
][0], tz
);
198 dpoly
fact(total_power
, tz
, 1);
202 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
203 zz2value(factors
[j
][0], tz
);
204 dpoly
fact(total_power
, tz
, 1);
208 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
209 zz2value(v
[i
]->terms
[k
][0], tz
);
210 dpoly
n(total_power
, tz
);
211 mpq_set_si(tcount
, 0, 1);
214 bfct
->c
[k
].n
= -bfct
->c
[k
].n
;
215 zz2value(bfct
->c
[k
].n
, tn
);
216 zz2value(bfct
->c
[k
].d
, td
);
218 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
219 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
220 mpq_canonicalize(tcount
);
221 mpq_add(count
, count
, tcount
);
228 /* Check whether the polyhedron is unbounded and if so,
229 * check whether it has any (and therefore an infinite number of)
231 * If one of the vertices is integer, then we are done.
232 * Otherwise, transform the polyhedron such that one of the rays
233 * is the first unit vector and cut it off at a height that ensures
234 * that if the whole polyhedron has any points, then the remaining part
235 * has integer points. In particular we add the largest coefficient
236 * of a ray to the highest vertex (rounded up).
238 static bool Polyhedron_is_infinite(Polyhedron
*P
, Value
* result
,
239 barvinok_options
*options
)
251 for (; r
< P
->NbRays
; ++r
)
252 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
254 if (P
->NbBid
== 0 && r
== P
->NbRays
)
257 if (options
->count_sample_infinite
) {
260 sample
= Polyhedron_Sample(P
, options
);
262 value_set_si(*result
, 0);
264 value_set_si(*result
, -1);
270 for (int i
= 0; i
< P
->NbRays
; ++i
)
271 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
272 value_set_si(*result
, -1);
277 M
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
278 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
279 Vector_AntiScale(P
->Ray
[r
]+1, M
->p
[0], g
, P
->Dimension
+1);
280 int ok
= unimodular_complete(M
, 1);
282 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
285 P
= Polyhedron_Preimage(P
, M2
, 0);
293 value_set_si(size
, 0);
295 for (int i
= 0; i
< P
->NbBid
; ++i
) {
296 value_absolute(tmp
, P
->Ray
[i
][1]);
297 if (value_gt(tmp
, size
))
298 value_assign(size
, tmp
);
300 for (int i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
301 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
302 if (value_gt(P
->Ray
[i
][1], size
))
303 value_assign(size
, P
->Ray
[i
][1]);
306 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
307 if (first
|| value_gt(tmp
, offset
)) {
308 value_assign(offset
, tmp
);
312 value_addto(offset
, offset
, size
);
316 v
= Vector_Alloc(P
->Dimension
+2);
317 value_set_si(v
->p
[0], 1);
318 value_set_si(v
->p
[1], -1);
319 value_assign(v
->p
[1+P
->Dimension
], offset
);
320 R
= AddConstraints(v
->p
, 1, P
, options
->MaxRays
);
328 barvinok_count_with_options(P
, &c
, options
);
331 value_set_si(*result
, 0);
333 value_set_si(*result
, -1);
339 static void evalue2value(evalue
*e
, Value
*v
)
341 if (EVALUE_IS_ZERO(*e
)) {
346 if (value_notzero_p(e
->d
)) {
347 assert(value_one_p(e
->d
));
348 value_assign(*v
, e
->x
.n
);
352 assert(e
->x
.p
->type
== partition
);
353 assert(e
->x
.p
->size
== 2);
354 assert(EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
== 0);
355 evalue2value(&e
->x
.p
->arr
[1], v
);
358 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
359 barvinok_options
*options
);
361 void barvinok_count_with_options(Polyhedron
*P
, Value
* result
,
362 struct barvinok_options
*options
)
367 bool infinite
= false;
371 "barvinok_count: input is a union; only first polyhedron is counted\n");
374 value_set_si(*result
, 0);
380 P
= remove_equalities(P
, options
->MaxRays
);
382 P
= DomainConstraintSimplify(P
, options
->MaxRays
);
386 } while (P
&& !emptyQ(P
) && P
->NbEq
!= 0);
387 if (!P
|| emptyQ(P
)) {
389 value_set_si(*result
, 0);
394 if (Polyhedron_is_infinite(P
, result
, options
)) {
399 if (P
->Dimension
== 0) {
400 /* Test whether the constraints are satisfied */
401 POL_ENSURE_VERTICES(P
);
402 value_set_si(*result
, !emptyQ(P
));
407 if (options
->summation
== BV_SUM_BERNOULLI
) {
408 Polyhedron
*C
= Universe_Polyhedron(0);
409 evalue
*sum
= barvinok_summate_unweighted(P
, C
, options
);
411 evalue2value(sum
, result
);
415 Q
= Polyhedron_Factor(P
, 0, NULL
, options
->MaxRays
);
423 barvinok_count_f(P
, result
, options
);
424 if (value_neg_p(*result
))
426 if (Q
&& P
->next
&& value_notzero_p(*result
)) {
430 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
431 barvinok_count_f(Q
, &factor
, options
);
432 if (value_neg_p(factor
)) {
435 } else if (Q
->next
&& value_zero_p(factor
)) {
436 value_set_si(*result
, 0);
439 value_multiply(*result
, *result
, factor
);
448 value_set_si(*result
, -1);
451 void barvinok_count(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
453 barvinok_options
*options
= barvinok_options_new_with_defaults();
454 options
->MaxRays
= NbMaxCons
;
455 barvinok_count_with_options(P
, result
, options
);
456 barvinok_options_free(options
);
459 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
460 barvinok_options
*options
)
463 value_set_si(*result
, 0);
467 if (P
->Dimension
== 1)
468 return Line_Length(P
, result
);
470 int c
= P
->NbConstraints
;
471 POL_ENSURE_FACETS(P
);
472 if (c
!= P
->NbConstraints
|| P
->NbEq
!= 0) {
473 Polyhedron
*next
= P
->next
;
475 barvinok_count_with_options(P
, result
, options
);
480 POL_ENSURE_VERTICES(P
);
482 if (Polyhedron_is_infinite(P
, result
, options
))
486 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
487 cnt
= new bfcounter(P
->Dimension
);
488 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
489 cnt
= new icounter(P
->Dimension
);
490 else if (options
->incremental_specialization
== BV_SPECIALIZATION_TODD
)
491 cnt
= new tcounter(P
->Dimension
, options
->max_index
);
493 cnt
= new counter(P
->Dimension
, options
->max_index
);
494 cnt
->start(P
, options
);
496 cnt
->get_count(result
);
500 typedef evalue
* evalue_p
;
502 struct enumerator_base
{
506 vertex_decomposer
*vpd
;
508 enumerator_base(unsigned dim
, vertex_decomposer
*vpd
)
513 vE
= new evalue_p
[vpd
->PP
->nbV
];
514 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
518 evalue_set_si(&mone
, -1, 1);
521 void decompose_at(Param_Vertices
*V
, int _i
, barvinok_options
*options
) {
525 value_init(vE
[_i
]->d
);
526 evalue_set_si(vE
[_i
], 0, 1);
528 vpd
->decompose_at_vertex(V
, _i
, options
);
531 virtual ~enumerator_base() {
532 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
534 free_evalue_refs(vE
[j
]);
539 free_evalue_refs(&mone
);
542 static enumerator_base
*create(Polyhedron
*P
, unsigned dim
,
543 Param_Polyhedron
*PP
,
544 barvinok_options
*options
);
547 struct enumerator
: public signed_cone_consumer
, public vertex_decomposer
,
548 public enumerator_base
{
556 enumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
557 vertex_decomposer(PP
, *this), enumerator_base(dim
, this) {
558 randomvector(P
, lambda
, dim
, 0);
560 c
= Vector_Alloc(dim
+2);
572 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
575 void enumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
579 assert(sc
.rays
.NumRows() == dim
);
580 for (int k
= 0; k
< dim
; ++k
) {
581 if (lambda
* sc
.rays
[k
] == 0)
585 lattice_point(V
, sc
.rays
, lambda
, &num
, sc
.det
, options
);
586 den
= sc
.rays
* lambda
;
591 zz2value(den
[0], tz
);
593 for (int k
= 1; k
< dim
; ++k
) {
594 zz2value(den
[k
], tz
);
595 dpoly
fact(dim
, tz
, 1);
601 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
602 evalue
*EV
= evalue_polynomial(c
, num
.E
[i
]);
605 evalue_free(num
.E
[i
]);
609 mpq_set_si(count
, 0, 1);
610 if (num
.constant
.length() == 1) {
611 zz2value(num
.constant
[0], tz
);
613 d
.div(n
, count
, sign
);
620 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
621 value_assign(acc
, c
->p
[dim
]);
622 zz2value(num
.constant
[i
], x
);
623 for (int j
= dim
-1; j
>= 0; --j
) {
624 value_multiply(acc
, acc
, x
);
625 value_addto(acc
, acc
, c
->p
[j
]);
627 value_addto(mpq_numref(count
), mpq_numref(count
), acc
);
629 mpz_set(mpq_denref(count
), c
->p
[dim
+1]);
635 evalue_set(&EV
, &count
[0]._mp_num
, &count
[0]._mp_den
);
637 free_evalue_refs(&EV
);
641 struct ienumerator_base
: enumerator_base
{
644 ienumerator_base(unsigned dim
, vertex_decomposer
*vpd
) :
645 enumerator_base(dim
,vpd
) {
646 E_vertex
= new evalue_p
[dim
];
649 virtual ~ienumerator_base() {
653 evalue
*E_num(int i
, int d
) {
654 return E_vertex
[i
+ (dim
-d
)];
663 cumulator(evalue
*factor
, evalue
*v
, dpoly_r
*r
) :
664 factor(factor
), v(v
), r(r
) {}
666 void cumulate(barvinok_options
*options
);
668 virtual void add_term(const vector
<int>& powers
, evalue
*f2
) = 0;
669 virtual ~cumulator() {}
672 void cumulator::cumulate(barvinok_options
*options
)
674 evalue cum
; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
676 evalue t
; // E_num[0] - (m-1)
680 if (options
->lookup_table
) {
682 evalue_set_si(&mone
, -1, 1);
686 evalue_copy(&cum
, factor
);
689 value_set_si(f
.d
, 1);
690 value_set_si(f
.x
.n
, 1);
694 if (!options
->lookup_table
) {
695 for (cst
= &t
; value_zero_p(cst
->d
); ) {
696 if (cst
->x
.p
->type
== fractional
)
697 cst
= &cst
->x
.p
->arr
[1];
699 cst
= &cst
->x
.p
->arr
[0];
703 for (int m
= 0; m
< r
->len
; ++m
) {
706 value_set_si(f
.d
, m
);
708 if (!options
->lookup_table
)
709 value_subtract(cst
->x
.n
, cst
->x
.n
, cst
->d
);
715 dpoly_r_term_list
& current
= r
->c
[r
->len
-1-m
];
716 dpoly_r_term_list::iterator j
;
717 for (j
= current
.begin(); j
!= current
.end(); ++j
) {
718 if ((*j
)->coeff
== 0)
720 evalue
*f2
= new evalue
;
723 zz2value((*j
)->coeff
, f2
->x
.n
);
724 zz2value(r
->denom
, f2
->d
);
727 add_term((*j
)->powers
, f2
);
730 free_evalue_refs(&f
);
731 free_evalue_refs(&t
);
732 free_evalue_refs(&cum
);
733 if (options
->lookup_table
)
734 free_evalue_refs(&mone
);
742 struct ie_cum
: public cumulator
{
743 vector
<E_poly_term
*> terms
;
745 ie_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
) : cumulator(factor
, v
, r
) {}
747 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
750 void ie_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
753 for (k
= 0; k
< terms
.size(); ++k
) {
754 if (terms
[k
]->powers
== powers
) {
755 eadd(f2
, terms
[k
]->E
);
756 free_evalue_refs(f2
);
761 if (k
>= terms
.size()) {
762 E_poly_term
*ET
= new E_poly_term
;
769 struct ienumerator
: public signed_cone_consumer
, public vertex_decomposer
,
770 public ienumerator_base
{
776 ienumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
777 vertex_decomposer(PP
, *this), ienumerator_base(dim
, this) {
778 vertex
.SetDims(1, dim
);
780 den
.SetDims(dim
, dim
);
790 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
791 void reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
792 barvinok_options
*options
);
795 void ienumerator::reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
796 barvinok_options
*options
)
798 unsigned len
= den_f
.NumRows(); // number of factors in den
799 unsigned dim
= num
.NumCols();
800 assert(num
.NumRows() == 1);
803 eadd(factor
, vE
[vert
]);
812 split_one(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
815 den_p
.SetLength(len
);
819 normalize(one
, num_s
, num_p
, den_s
, den_p
, den_r
);
825 for (int k
= 0; k
< len
; ++k
) {
828 else if (den_s
[k
] == 0)
832 reduce(factor
, num_p
, den_r
, options
);
836 pden
.SetDims(only_param
, dim
-1);
838 for (k
= 0, l
= 0; k
< len
; ++k
)
840 pden
[l
++] = den_r
[k
];
842 for (k
= 0; k
< len
; ++k
)
846 zz2value(num_s
[0], tz
);
847 dpoly
n(no_param
, tz
);
848 zz2value(den_s
[k
], tz
);
849 dpoly
D(no_param
, tz
, 1);
852 zz2value(den_s
[k
], tz
);
853 dpoly
fact(no_param
, tz
, 1);
858 // if no_param + only_param == len then all powers
859 // below will be all zero
860 if (no_param
+ only_param
== len
) {
861 if (E_num(0, dim
) != 0)
862 r
= new dpoly_r(n
, len
);
864 mpq_set_si(tcount
, 0, 1);
868 if (value_notzero_p(mpq_numref(tcount
))) {
872 value_assign(f
.x
.n
, mpq_numref(tcount
));
873 value_assign(f
.d
, mpq_denref(tcount
));
875 reduce(factor
, num_p
, pden
, options
);
876 free_evalue_refs(&f
);
881 for (k
= 0; k
< len
; ++k
) {
882 if (den_s
[k
] == 0 || den_p
[k
] == 0)
885 zz2value(den_s
[k
], tz
);
886 dpoly
pd(no_param
-1, tz
, 1);
889 for (l
= 0; l
< k
; ++l
)
890 if (den_r
[l
] == den_r
[k
])
894 r
= new dpoly_r(n
, pd
, l
, len
);
896 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
902 dpoly_r
*rc
= r
->div(D
);
905 if (E_num(0, dim
) == 0) {
906 int common
= pden
.NumRows();
907 dpoly_r_term_list
& final
= r
->c
[r
->len
-1];
913 zz2value(r
->denom
, f
.d
);
914 dpoly_r_term_list::iterator j
;
915 for (j
= final
.begin(); j
!= final
.end(); ++j
) {
916 if ((*j
)->coeff
== 0)
919 for (int k
= 0; k
< r
->dim
; ++k
) {
920 int n
= (*j
)->powers
[k
];
923 pden
.SetDims(rows
+n
, pden
.NumCols());
924 for (int l
= 0; l
< n
; ++l
)
925 pden
[rows
+l
] = den_r
[k
];
929 evalue_copy(&t
, factor
);
930 zz2value((*j
)->coeff
, f
.x
.n
);
932 reduce(&t
, num_p
, pden
, options
);
933 free_evalue_refs(&t
);
935 free_evalue_refs(&f
);
937 ie_cum
cum(factor
, E_num(0, dim
), r
);
938 cum
.cumulate(options
);
940 int common
= pden
.NumRows();
942 for (int j
= 0; j
< cum
.terms
.size(); ++j
) {
944 pden
.SetDims(rows
, pden
.NumCols());
945 for (int k
= 0; k
< r
->dim
; ++k
) {
946 int n
= cum
.terms
[j
]->powers
[k
];
949 pden
.SetDims(rows
+n
, pden
.NumCols());
950 for (int l
= 0; l
< n
; ++l
)
951 pden
[rows
+l
] = den_r
[k
];
954 reduce(cum
.terms
[j
]->E
, num_p
, pden
, options
);
955 free_evalue_refs(cum
.terms
[j
]->E
);
956 delete cum
.terms
[j
]->E
;
964 void ienumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
967 assert(sc
.rays
.NumRows() == dim
);
969 lattice_point(V
, sc
.rays
, vertex
[0], E_vertex
, options
);
975 evalue_set_si(&one
, sc
.sign
, 1);
976 reduce(&one
, vertex
, den
, options
);
977 free_evalue_refs(&one
);
979 for (int i
= 0; i
< dim
; ++i
)
981 evalue_free(E_vertex
[i
]);
984 struct bfenumerator
: public vertex_decomposer
, public bf_base
,
985 public ienumerator_base
{
988 bfenumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
989 vertex_decomposer(PP
, *this),
990 bf_base(dim
), ienumerator_base(dim
, this) {
998 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
999 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1001 bfc_term_base
* new_bf_term(int len
) {
1002 bfe_term
* t
= new bfe_term(len
);
1006 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
1007 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1008 factor
= bfet
->factors
[k
];
1009 assert(factor
!= NULL
);
1010 bfet
->factors
[k
] = NULL
;
1012 emul(&mone
, factor
);
1015 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&q
, int change
) {
1016 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1017 factor
= bfet
->factors
[k
];
1018 assert(factor
!= NULL
);
1019 bfet
->factors
[k
] = NULL
;
1025 value_oppose(f
.x
.n
, mpq_numref(q
));
1027 value_assign(f
.x
.n
, mpq_numref(q
));
1028 value_assign(f
.d
, mpq_denref(q
));
1030 free_evalue_refs(&f
);
1033 virtual void set_factor(bfc_term_base
*t
, int k
, const QQ
& c
, int change
) {
1034 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1036 factor
= new evalue
;
1041 zz2value(c
.n
, f
.x
.n
);
1043 value_oppose(f
.x
.n
, f
.x
.n
);
1046 value_init(factor
->d
);
1047 evalue_copy(factor
, bfet
->factors
[k
]);
1049 free_evalue_refs(&f
);
1052 void set_factor(evalue
*f
, int change
) {
1058 virtual void insert_term(bfc_term_base
*t
, int i
) {
1059 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1060 int len
= t
->terms
.NumRows()-1; // already increased by one
1062 bfet
->factors
.resize(len
+1);
1063 for (int j
= len
; j
> i
; --j
) {
1064 bfet
->factors
[j
] = bfet
->factors
[j
-1];
1065 t
->terms
[j
] = t
->terms
[j
-1];
1067 bfet
->factors
[i
] = factor
;
1071 virtual void update_term(bfc_term_base
*t
, int i
) {
1072 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1074 eadd(factor
, bfet
->factors
[i
]);
1075 free_evalue_refs(factor
);
1079 virtual bool constant_vertex(int dim
) { return E_num(0, dim
) == 0; }
1081 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
,
1082 barvinok_options
*options
);
1085 enumerator_base
*enumerator_base::create(Polyhedron
*P
, unsigned dim
,
1086 Param_Polyhedron
*PP
,
1087 barvinok_options
*options
)
1089 enumerator_base
*eb
;
1091 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
1092 eb
= new bfenumerator(P
, dim
, PP
);
1093 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
1094 eb
= new ienumerator(P
, dim
, PP
);
1096 eb
= new enumerator(P
, dim
, PP
);
1101 struct bfe_cum
: public cumulator
{
1103 bfc_term_base
*told
;
1107 bfe_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
, bf_reducer
*bfr
,
1108 bfc_term_base
*t
, int k
, bfenumerator
*e
) :
1109 cumulator(factor
, v
, r
), told(t
), k(k
),
1113 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
1116 void bfe_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
1118 bfr
->update_powers(powers
);
1120 bfc_term_base
* t
= bfe
->find_bfc_term(bfr
->vn
, bfr
->npowers
, bfr
->nnf
);
1121 bfe
->set_factor(f2
, bfr
->l_changes
% 2);
1122 bfe
->add_term(t
, told
->terms
[k
], bfr
->l_extra_num
);
1125 void bfenumerator::cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
1126 dpoly_r
*r
, barvinok_options
*options
)
1128 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1129 bfe_cum
cum(bfet
->factors
[k
], E_num(0, bfr
->d
), r
, bfr
, t
, k
, this);
1130 cum
.cumulate(options
);
1133 void bfenumerator::base(mat_ZZ
& factors
, bfc_vec
& v
)
1135 for (int i
= 0; i
< v
.size(); ++i
) {
1136 assert(v
[i
]->terms
.NumRows() == 1);
1137 evalue
*factor
= static_cast<bfe_term
*>(v
[i
])->factors
[0];
1138 eadd(factor
, vE
[vert
]);
1143 void bfenumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1145 assert(sc
.det
== 1);
1146 assert(sc
.rays
.NumRows() == enumerator_base::dim
);
1148 bfe_term
* t
= new bfe_term(enumerator_base::dim
);
1149 vector
< bfc_term_base
* > v
;
1152 t
->factors
.resize(1);
1154 t
->terms
.SetDims(1, enumerator_base::dim
);
1155 lattice_point(V
, sc
.rays
, t
->terms
[0], E_vertex
, options
);
1157 // the elements of factors are always lexpositive
1159 int s
= setup_factors(sc
.rays
, factors
, t
, sc
.sign
);
1161 t
->factors
[0] = new evalue
;
1162 value_init(t
->factors
[0]->d
);
1163 evalue_set_si(t
->factors
[0], s
, 1);
1164 reduce(factors
, v
, options
);
1166 for (int i
= 0; i
< enumerator_base::dim
; ++i
)
1168 evalue_free(E_vertex
[i
]);
1171 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1172 barvinok_options
*options
);
1175 static evalue
* barvinok_enumerate_cst(Polyhedron
*P
, Polyhedron
* C
,
1176 struct barvinok_options
*options
)
1182 return evalue_zero();
1185 ALLOC(evalue
, eres
);
1186 value_init(eres
->d
);
1187 value_set_si(eres
->d
, 0);
1188 eres
->x
.p
= new_enode(partition
, 2, C
->Dimension
);
1189 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[0],
1190 DomainConstraintSimplify(C
, options
->MaxRays
));
1191 value_set_si(eres
->x
.p
->arr
[1].d
, 1);
1192 value_init(eres
->x
.p
->arr
[1].x
.n
);
1194 value_set_si(eres
->x
.p
->arr
[1].x
.n
, 0);
1196 barvinok_count_with_options(P
, &eres
->x
.p
->arr
[1].x
.n
, options
);
1197 if (value_mone_p(eres
->x
.p
->arr
[1].x
.n
)) {
1198 value_clear(eres
->x
.p
->arr
[1].x
.n
);
1199 value_set_si(eres
->x
.p
->arr
[1].d
, -2); /* NaN */
1205 static evalue
* enumerate(Polyhedron
*P
, Polyhedron
* C
,
1206 struct barvinok_options
*options
)
1209 Polyhedron
*Porig
= P
;
1210 Polyhedron
*Corig
= C
;
1211 Polyhedron
*CEq
= NULL
, *rVD
;
1213 unsigned nparam
= C
->Dimension
;
1218 value_init(factor
.d
);
1219 evalue_set_si(&factor
, 1, 1);
1222 POL_ENSURE_FACETS(P
);
1223 POL_ENSURE_VERTICES(P
);
1224 POL_ENSURE_FACETS(C
);
1225 POL_ENSURE_VERTICES(C
);
1227 if (C
->Dimension
== 0 || emptyQ(P
) || emptyQ(C
)) {
1230 CEq
= Polyhedron_Copy(CEq
);
1231 eres
= barvinok_enumerate_cst(P
, CEq
? CEq
: Polyhedron_Copy(C
), options
);
1234 evalue_backsubstitute(eres
, CP
, options
->MaxRays
);
1238 emul(&factor
, eres
);
1239 if (options
->approx
->method
== BV_APPROX_DROP
) {
1240 if (options
->approx
->approximation
== BV_APPROX_SIGN_UPPER
)
1241 evalue_frac2polynomial(eres
, 1, options
->MaxRays
);
1242 if (options
->approx
->approximation
== BV_APPROX_SIGN_LOWER
)
1243 evalue_frac2polynomial(eres
, -1, options
->MaxRays
);
1244 if (options
->approx
->approximation
== BV_APPROX_SIGN_APPROX
)
1245 evalue_frac2polynomial(eres
, 0, options
->MaxRays
);
1247 reduce_evalue(eres
);
1248 free_evalue_refs(&factor
);
1256 if (Polyhedron_is_unbounded(P
, nparam
, options
->MaxRays
))
1259 if (P
->Dimension
== nparam
) {
1260 CEq
= DomainIntersection(P
, C
, options
->MaxRays
);
1261 P
= Universe_Polyhedron(0);
1264 if (P
->NbEq
!= 0 || C
->NbEq
!= 0) {
1267 remove_all_equalities(&P
, &C
, &CP
, NULL
, nparam
, options
->MaxRays
);
1268 if (C
!= D
&& D
!= Corig
)
1270 if (P
!= Q
&& Q
!= Porig
)
1272 eres
= enumerate(P
, C
, options
);
1276 Polyhedron
*T
= Polyhedron_Factor(P
, nparam
, NULL
, options
->MaxRays
);
1277 if (T
|| (P
->Dimension
== nparam
+1)) {
1279 Polyhedron
*FC
= Factor_Context(T
? T
: P
, nparam
, options
->MaxRays
);
1280 C
= DomainIntersection(C
, FC
, options
->MaxRays
);
1282 Polyhedron_Free(C2
);
1283 Polyhedron_Free(FC
);
1289 if (T
->Dimension
== C
->Dimension
) {
1298 eres
= barvinok_enumerate_ev_f(P
, C
, options
);
1305 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
1306 Polyhedron
*next
= Q
->next
;
1309 f
= barvinok_enumerate_ev_f(Q
, C
, options
);
1320 evalue
* barvinok_enumerate_with_options(Polyhedron
*P
, Polyhedron
* C
,
1321 struct barvinok_options
*options
)
1323 Polyhedron
*next
, *Cnext
, *C1
;
1324 Polyhedron
*Corig
= C
;
1329 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1333 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1337 C1
= Polyhedron_Project(P
, C
->Dimension
);
1338 C
= DomainIntersection(C
, C1
, options
->MaxRays
);
1339 Polyhedron_Free(C1
);
1343 if (options
->approx
->method
== BV_APPROX_BERNOULLI
||
1344 options
->summation
== BV_SUM_BERNOULLI
) {
1345 int summation
= options
->summation
;
1346 options
->summation
= BV_SUM_BERNOULLI
;
1347 eres
= barvinok_summate_unweighted(P
, C
, options
);
1348 options
->summation
= summation
;
1350 eres
= enumerate(P
, C
, options
);
1354 Corig
->next
= Cnext
;
1359 evalue
* barvinok_enumerate_ev(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1362 barvinok_options
*options
= barvinok_options_new_with_defaults();
1363 options
->MaxRays
= MaxRays
;
1364 E
= barvinok_enumerate_with_options(P
, C
, options
);
1365 barvinok_options_free(options
);
1369 evalue
*Param_Polyhedron_Enumerate(Param_Polyhedron
*PP
, Polyhedron
*P
,
1371 struct barvinok_options
*options
)
1375 unsigned nparam
= C
->Dimension
;
1376 unsigned dim
= P
->Dimension
- nparam
;
1379 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
1380 evalue_section
*s
= new evalue_section
[nd
];
1381 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
1383 enumerator_base
*et
= NULL
;
1388 et
= enumerator_base::create(P
, dim
, PP
, options
);
1390 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
1393 s
[i
].E
= evalue_zero();
1396 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
1399 et
->decompose_at(V
, _i
, options
);
1400 } catch (OrthogonalException
&e
) {
1401 FORALL_REDUCED_DOMAIN_RESET
;
1402 for (; i
>= 0; --i
) {
1403 evalue_free(s
[i
].E
);
1404 Domain_Free(s
[i
].D
);
1408 eadd(et
->vE
[_i
] , s
[i
].E
);
1409 END_FORALL_PVertex_in_ParamPolyhedron
;
1410 evalue_range_reduction_in_domain(s
[i
].E
, rVD
);
1411 END_FORALL_REDUCED_DOMAIN
1412 Polyhedron_Free(TC
);
1415 eres
= evalue_from_section_array(s
, nd
);
1421 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1422 barvinok_options
*options
)
1424 unsigned nparam
= C
->Dimension
;
1425 bool do_scale
= options
->approx
->method
== BV_APPROX_SCALE
;
1427 if (options
->summation
== BV_SUM_EULER
)
1428 return barvinok_summate_unweighted(P
, C
, options
);
1430 if (options
->approx
->method
== BV_APPROX_VOLUME
)
1431 return Param_Polyhedron_Volume(P
, C
, options
);
1433 if (P
->Dimension
- nparam
== 1 && !do_scale
)
1434 return ParamLine_Length(P
, C
, options
);
1436 Param_Polyhedron
*PP
= NULL
;
1440 eres
= scale_bound(P
, C
, options
);
1445 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
1448 eres
= scale(PP
, P
, C
, options
);
1450 eres
= Param_Polyhedron_Enumerate(PP
, P
, C
, options
);
1453 Param_Polyhedron_Free(PP
);
1458 Enumeration
* barvinok_enumerate(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1460 evalue
*EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
1462 return partition2enumeration(EP
);
1465 evalue
* barvinok_enumerate_union(Polyhedron
*D
, Polyhedron
* C
, unsigned MaxRays
)
1468 gen_fun
*gf
= barvinok_enumerate_union_series(D
, C
, MaxRays
);
1474 static __isl_give isl_pw_qpolynomial
*basic_set_card(
1475 __isl_take isl_basic_set
*bset
)
1479 isl_pw_qpolynomial
*pwqp
;
1480 unsigned nparam
= isl_basic_set_dim(bset
, isl_dim_param
);
1481 Polyhedron
*U
= Universe_Polyhedron(nparam
);
1484 barvinok_options
*options
;
1485 int options_allocated
= 0;
1487 ctx
= isl_basic_set_get_ctx(bset
);
1488 options
= isl_ctx_peek_barvinok_options(ctx
);
1490 options
= barvinok_options_new_with_defaults();
1491 options_allocated
= 1;
1494 dim
= isl_basic_set_get_space(bset
);
1495 dim
= isl_space_domain(dim
);
1497 P
= isl_basic_set_to_polylib(bset
);
1498 E
= enumerate(P
, U
, options
);
1500 pwqp
= isl_pw_qpolynomial_from_evalue(dim
, E
);
1501 isl_basic_set_free(bset
);
1506 if (options_allocated
)
1507 barvinok_options_free(options
);
1512 static isl_stat
basic_map_card(__isl_take isl_basic_map
*bmap
, void *user
)
1514 isl_pw_qpolynomial
**sum
= (isl_pw_qpolynomial
**)user
;
1515 isl_pw_qpolynomial
*pwqp
;
1516 unsigned nparam
= isl_basic_map_dim(bmap
, isl_dim_param
);
1517 unsigned n_in
= isl_basic_map_dim(bmap
, isl_dim_in
);
1518 isl_space
*target_dim
;
1519 isl_basic_set
*bset
;
1521 target_dim
= isl_basic_map_get_space(bmap
);
1522 target_dim
= isl_space_domain(target_dim
);
1524 bmap
= isl_basic_map_move_dims(bmap
, isl_dim_param
, nparam
,
1525 isl_dim_in
, 0, n_in
);
1527 bset
= isl_basic_map_range(bmap
);
1528 bset
= isl_basic_set_lift(bset
);
1529 pwqp
= isl_basic_set_multiplicative_call(bset
, &basic_set_card
);
1531 pwqp
= isl_pw_qpolynomial_move_dims(pwqp
, isl_dim_in
, 0,
1532 isl_dim_param
, nparam
, n_in
);
1533 pwqp
= isl_pw_qpolynomial_reset_domain_space(pwqp
, target_dim
);
1534 *sum
= isl_pw_qpolynomial_add(*sum
, pwqp
);
1539 static __isl_give isl_pw_qpolynomial
*card_as_sum(__isl_take isl_map
*map
,
1540 barvinok_options
*options
)
1546 isl_qpolynomial
*qp
;
1547 isl_pw_qpolynomial
*pwqp
;
1548 int summation
= options
->summation
;
1553 options
->summation
= BV_SUM_BERNOULLI
;
1555 set
= isl_map_wrap(map
);
1556 dim
= isl_set_get_space(set
);
1557 ctx
= isl_map_get_ctx(map
);
1558 one
= isl_val_one(ctx
);
1559 qp
= isl_qpolynomial_val_on_domain(dim
, one
);
1561 pwqp
= isl_pw_qpolynomial_alloc(set
, qp
);
1562 pwqp
= isl_pw_qpolynomial_sum(pwqp
);
1564 options
->summation
= summation
;
1569 __isl_give isl_pw_qpolynomial
*isl_map_card(__isl_take isl_map
*map
)
1573 isl_pw_qpolynomial
*sum
;
1574 barvinok_options
*options
;
1576 ctx
= isl_map_get_ctx(map
);
1577 options
= isl_ctx_peek_barvinok_options(ctx
);
1579 (options
->approx
->method
== BV_APPROX_BERNOULLI
||
1580 options
->summation
== BV_SUM_BERNOULLI
))
1581 return card_as_sum(map
, options
);
1583 dim
= isl_map_get_space(map
);
1584 dim
= isl_space_domain(dim
);
1585 dim
= isl_space_from_domain(dim
);
1586 dim
= isl_space_add_dims(dim
, isl_dim_out
, 1);
1587 sum
= isl_pw_qpolynomial_zero(dim
);
1589 map
= isl_map_make_disjoint(map
);
1590 map
= isl_map_compute_divs(map
);
1592 if (isl_map_foreach_basic_map(map
, &basic_map_card
, &sum
) < 0)
1600 isl_pw_qpolynomial_free(sum
);
1604 __isl_give isl_pw_qpolynomial
*isl_set_card(__isl_take isl_set
*set
)
1606 isl_pw_qpolynomial
*pwqp
;
1607 pwqp
= isl_map_card(isl_map_from_range(set
));
1608 pwqp
= isl_pw_qpolynomial_project_domain_on_params(pwqp
);
1612 __isl_give isl_pw_qpolynomial
*isl_basic_map_card(__isl_take isl_basic_map
*bmap
)
1614 return isl_map_card(isl_map_from_basic_map(bmap
));
1617 __isl_give isl_pw_qpolynomial
*isl_basic_set_card(__isl_take isl_basic_set
*bset
)
1619 isl_pw_qpolynomial
*pwqp
;
1620 pwqp
= isl_basic_map_card(isl_basic_map_from_range(bset
));
1621 pwqp
= isl_pw_qpolynomial_project_domain_on_params(pwqp
);
1625 static isl_stat
set_card(__isl_take isl_set
*set
, void *user
)
1627 isl_union_pw_qpolynomial
**res
= (isl_union_pw_qpolynomial
**)user
;
1628 isl_pw_qpolynomial
*pwqp
;
1629 isl_union_pw_qpolynomial
*upwqp
;
1631 pwqp
= isl_set_card(set
);
1632 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1633 *res
= isl_union_pw_qpolynomial_add(*res
, upwqp
);
1638 __isl_give isl_union_pw_qpolynomial
*isl_union_set_card(
1639 __isl_take isl_union_set
*uset
)
1642 isl_union_pw_qpolynomial
*res
;
1644 dim
= isl_union_set_get_space(uset
);
1645 res
= isl_union_pw_qpolynomial_zero(dim
);
1646 if (isl_union_set_foreach_set(uset
, &set_card
, &res
) < 0)
1648 isl_union_set_free(uset
);
1652 isl_union_set_free(uset
);
1653 isl_union_pw_qpolynomial_free(res
);
1657 static isl_stat
map_card(__isl_take isl_map
*map
, void *user
)
1659 isl_union_pw_qpolynomial
**res
= (isl_union_pw_qpolynomial
**)user
;
1660 isl_pw_qpolynomial
*pwqp
;
1661 isl_union_pw_qpolynomial
*upwqp
;
1663 pwqp
= isl_map_card(map
);
1664 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1665 *res
= isl_union_pw_qpolynomial_add(*res
, upwqp
);
1670 __isl_give isl_union_pw_qpolynomial
*isl_union_map_card(
1671 __isl_take isl_union_map
*umap
)
1674 isl_union_pw_qpolynomial
*res
;
1676 dim
= isl_union_map_get_space(umap
);
1677 res
= isl_union_pw_qpolynomial_zero(dim
);
1678 if (isl_union_map_foreach_map(umap
, &map_card
, &res
) < 0)
1680 isl_union_map_free(umap
);
1684 isl_union_map_free(umap
);
1685 isl_union_pw_qpolynomial_free(res
);