8 #include <NTL/mat_ZZ.h>
10 #include <barvinok/util.h>
11 #include <barvinok/evalue.h>
13 #include <barvinok/barvinok.h>
14 #include <barvinok/genfun.h>
15 #include <barvinok/options.h>
16 #include <barvinok/sample.h>
17 #include "bfcounter.h"
18 #include "conversion.h"
20 #include "decomposer.h"
22 #include "lattice_point.h"
24 #include "reduce_domain.h"
25 #include "remove_equalities.h"
28 #include "bernoulli.h"
29 #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);
100 * Searches for a vector that is not orthogonal to any
101 * of the rays in rays.
103 static void nonorthog(mat_ZZ
& rays
, vec_ZZ
& lambda
)
105 int dim
= rays
.NumCols();
107 lambda
.SetLength(dim
);
111 for (int i
= 2; !found
&& i
<= 50*dim
; i
+=4) {
112 for (int j
= 0; j
< MAX_TRY
; ++j
) {
113 for (int k
= 0; k
< dim
; ++k
) {
114 int r
= random_int(i
)+2;
115 int v
= (2*(r
%2)-1) * (r
>> 1);
119 for (; k
< rays
.NumRows(); ++k
)
120 if (lambda
* rays
[k
] == 0)
122 if (k
== rays
.NumRows()) {
131 static void add_rays(mat_ZZ
& rays
, Polyhedron
*i
, int *r
, int nvar
= -1,
134 unsigned dim
= i
->Dimension
;
137 for (int k
= 0; k
< i
->NbRays
; ++k
) {
138 if (!value_zero_p(i
->Ray
[k
][dim
+1]))
140 if (!all
&& nvar
!= dim
&& First_Non_Zero(i
->Ray
[k
]+1, nvar
) == -1)
142 values2zz(i
->Ray
[k
]+1, rays
[(*r
)++], nvar
);
146 struct bfe_term
: public bfc_term_base
{
147 vector
<evalue
*> factors
;
149 bfe_term(int len
) : bfc_term_base(len
) {
153 for (int i
= 0; i
< factors
.size(); ++i
) {
156 free_evalue_refs(factors
[i
]);
162 static void print_int_vector(int *v
, int len
, const char *name
)
164 cerr
<< name
<< endl
;
165 for (int j
= 0; j
< len
; ++j
) {
171 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
174 cerr
<< "factors" << endl
;
175 cerr
<< factors
<< endl
;
176 for (int i
= 0; i
< v
.size(); ++i
) {
177 cerr
<< "term: " << i
<< endl
;
178 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
179 cerr
<< "terms" << endl
;
180 cerr
<< v
[i
]->terms
<< endl
;
181 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
182 cerr
<< bfct
->c
<< endl
;
186 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
189 cerr
<< "factors" << endl
;
190 cerr
<< factors
<< endl
;
191 for (int i
= 0; i
< v
.size(); ++i
) {
192 cerr
<< "term: " << i
<< endl
;
193 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
194 cerr
<< "terms" << endl
;
195 cerr
<< v
[i
]->terms
<< endl
;
196 bfe_term
* bfet
= static_cast<bfe_term
*>(v
[i
]);
197 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
198 const char * test
[] = {"a", "b"};
199 print_evalue(stderr
, bfet
->factors
[j
], test
);
200 fprintf(stderr
, "\n");
205 struct bfcounter
: public bfcounter_base
{
209 bfcounter(unsigned dim
) : bfcounter_base(dim
) {
218 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
219 virtual void get_count(Value
*result
) {
220 assert(value_one_p(&count
[0]._mp_den
));
221 value_assign(*result
, &count
[0]._mp_num
);
225 void bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
227 unsigned nf
= factors
.NumRows();
229 for (int i
= 0; i
< v
.size(); ++i
) {
230 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
232 // factor is always positive, so we always
234 for (int k
= 0; k
< nf
; ++k
)
235 total_power
+= v
[i
]->powers
[k
];
238 for (j
= 0; j
< nf
; ++j
)
239 if (v
[i
]->powers
[j
] > 0)
242 zz2value(factors
[j
][0], tz
);
243 dpoly
D(total_power
, tz
, 1);
244 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
245 zz2value(factors
[j
][0], tz
);
246 dpoly
fact(total_power
, tz
, 1);
250 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
251 zz2value(factors
[j
][0], tz
);
252 dpoly
fact(total_power
, tz
, 1);
256 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
257 zz2value(v
[i
]->terms
[k
][0], tz
);
258 dpoly
n(total_power
, tz
);
259 mpq_set_si(tcount
, 0, 1);
262 bfct
->c
[k
].n
= -bfct
->c
[k
].n
;
263 zz2value(bfct
->c
[k
].n
, tn
);
264 zz2value(bfct
->c
[k
].d
, td
);
266 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
267 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
268 mpq_canonicalize(tcount
);
269 mpq_add(count
, count
, tcount
);
276 /* Check whether the polyhedron is unbounded and if so,
277 * check whether it has any (and therefore an infinite number of)
279 * If one of the vertices is integer, then we are done.
280 * Otherwise, transform the polyhedron such that one of the rays
281 * is the first unit vector and cut it off at a height that ensures
282 * that if the whole polyhedron has any points, then the remaining part
283 * has integer points. In particular we add the largest coefficient
284 * of a ray to the highest vertex (rounded up).
286 static bool Polyhedron_is_infinite(Polyhedron
*P
, Value
* result
,
287 barvinok_options
*options
)
299 for (; r
< P
->NbRays
; ++r
)
300 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
302 if (P
->NbBid
== 0 && r
== P
->NbRays
)
305 if (options
->count_sample_infinite
) {
308 sample
= Polyhedron_Sample(P
, options
);
310 value_set_si(*result
, 0);
312 value_set_si(*result
, -1);
318 for (int i
= 0; i
< P
->NbRays
; ++i
)
319 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
320 value_set_si(*result
, -1);
325 M
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
326 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
327 Vector_AntiScale(P
->Ray
[r
]+1, M
->p
[0], g
, P
->Dimension
+1);
328 int ok
= unimodular_complete(M
, 1);
330 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
333 P
= Polyhedron_Preimage(P
, M2
, 0);
341 value_set_si(size
, 0);
343 for (int i
= 0; i
< P
->NbBid
; ++i
) {
344 value_absolute(tmp
, P
->Ray
[i
][1]);
345 if (value_gt(tmp
, size
))
346 value_assign(size
, tmp
);
348 for (int i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
349 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
350 if (value_gt(P
->Ray
[i
][1], size
))
351 value_assign(size
, P
->Ray
[i
][1]);
354 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
355 if (first
|| value_gt(tmp
, offset
)) {
356 value_assign(offset
, tmp
);
360 value_addto(offset
, offset
, size
);
364 v
= Vector_Alloc(P
->Dimension
+2);
365 value_set_si(v
->p
[0], 1);
366 value_set_si(v
->p
[1], -1);
367 value_assign(v
->p
[1+P
->Dimension
], offset
);
368 R
= AddConstraints(v
->p
, 1, P
, options
->MaxRays
);
376 barvinok_count_with_options(P
, &c
, options
);
379 value_set_si(*result
, 0);
381 value_set_si(*result
, -1);
387 static void evalue2value(evalue
*e
, Value
*v
)
389 if (EVALUE_IS_ZERO(*e
)) {
394 if (value_notzero_p(e
->d
)) {
395 assert(value_one_p(e
->d
));
396 value_assign(*v
, e
->x
.n
);
400 assert(e
->x
.p
->type
== partition
);
401 assert(e
->x
.p
->size
== 2);
402 assert(EVALUE_DOMAIN(e
->x
.p
->arr
[0])->Dimension
== 0);
403 evalue2value(&e
->x
.p
->arr
[1], v
);
406 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
407 barvinok_options
*options
);
409 void barvinok_count_with_options(Polyhedron
*P
, Value
* result
,
410 struct barvinok_options
*options
)
415 bool infinite
= false;
419 "barvinok_count: input is a union; only first polyhedron is counted\n");
422 value_set_si(*result
, 0);
428 P
= remove_equalities(P
, options
->MaxRays
);
429 P
= DomainConstraintSimplify(P
, options
->MaxRays
);
433 } while (!emptyQ(P
) && P
->NbEq
!= 0);
436 value_set_si(*result
, 0);
441 if (Polyhedron_is_infinite(P
, result
, options
)) {
446 if (P
->Dimension
== 0) {
447 /* Test whether the constraints are satisfied */
448 POL_ENSURE_VERTICES(P
);
449 value_set_si(*result
, !emptyQ(P
));
454 if (options
->summation
== BV_SUM_BERNOULLI
) {
455 Polyhedron
*C
= Universe_Polyhedron(0);
456 evalue
*sum
= barvinok_summate_unweighted(P
, C
, options
);
458 evalue2value(sum
, result
);
462 Q
= Polyhedron_Factor(P
, 0, NULL
, options
->MaxRays
);
470 barvinok_count_f(P
, result
, options
);
471 if (value_neg_p(*result
))
473 if (Q
&& P
->next
&& value_notzero_p(*result
)) {
477 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
478 barvinok_count_f(Q
, &factor
, options
);
479 if (value_neg_p(factor
)) {
482 } else if (Q
->next
&& value_zero_p(factor
)) {
483 value_set_si(*result
, 0);
486 value_multiply(*result
, *result
, factor
);
495 value_set_si(*result
, -1);
498 void barvinok_count(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
500 barvinok_options
*options
= barvinok_options_new_with_defaults();
501 options
->MaxRays
= NbMaxCons
;
502 barvinok_count_with_options(P
, result
, options
);
503 barvinok_options_free(options
);
506 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
507 barvinok_options
*options
)
510 value_set_si(*result
, 0);
514 if (P
->Dimension
== 1)
515 return Line_Length(P
, result
);
517 int c
= P
->NbConstraints
;
518 POL_ENSURE_FACETS(P
);
519 if (c
!= P
->NbConstraints
|| P
->NbEq
!= 0) {
520 Polyhedron
*next
= P
->next
;
522 barvinok_count_with_options(P
, result
, options
);
527 POL_ENSURE_VERTICES(P
);
529 if (Polyhedron_is_infinite(P
, result
, options
))
533 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
534 cnt
= new bfcounter(P
->Dimension
);
535 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
536 cnt
= new icounter(P
->Dimension
);
537 else if (options
->incremental_specialization
== BV_SPECIALIZATION_TODD
)
538 cnt
= new tcounter(P
->Dimension
, options
->max_index
);
540 cnt
= new counter(P
->Dimension
, options
->max_index
);
541 cnt
->start(P
, options
);
543 cnt
->get_count(result
);
547 static void uni_polynom(int param
, Vector
*c
, evalue
*EP
)
549 unsigned dim
= c
->Size
-2;
551 value_set_si(EP
->d
,0);
552 EP
->x
.p
= new_enode(polynomial
, dim
+1, param
+1);
553 for (int j
= 0; j
<= dim
; ++j
)
554 evalue_set(&EP
->x
.p
->arr
[j
], c
->p
[j
], c
->p
[dim
+1]);
557 typedef evalue
* evalue_p
;
559 struct enumerator_base
{
563 vertex_decomposer
*vpd
;
565 enumerator_base(unsigned dim
, vertex_decomposer
*vpd
)
570 vE
= new evalue_p
[vpd
->PP
->nbV
];
571 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
575 evalue_set_si(&mone
, -1, 1);
578 void decompose_at(Param_Vertices
*V
, int _i
, barvinok_options
*options
) {
582 value_init(vE
[_i
]->d
);
583 evalue_set_si(vE
[_i
], 0, 1);
585 vpd
->decompose_at_vertex(V
, _i
, options
);
588 virtual ~enumerator_base() {
589 for (int j
= 0; j
< vpd
->PP
->nbV
; ++j
)
591 free_evalue_refs(vE
[j
]);
596 free_evalue_refs(&mone
);
599 static enumerator_base
*create(Polyhedron
*P
, unsigned dim
,
600 Param_Polyhedron
*PP
,
601 barvinok_options
*options
);
604 struct enumerator
: public signed_cone_consumer
, public vertex_decomposer
,
605 public enumerator_base
{
613 enumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
614 vertex_decomposer(PP
, *this), enumerator_base(dim
, this) {
615 randomvector(P
, lambda
, dim
);
617 c
= Vector_Alloc(dim
+2);
629 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
632 void enumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
636 assert(sc
.rays
.NumRows() == dim
);
637 for (int k
= 0; k
< dim
; ++k
) {
638 if (lambda
* sc
.rays
[k
] == 0)
642 lattice_point(V
, sc
.rays
, lambda
, &num
, sc
.det
, options
);
643 den
= sc
.rays
* lambda
;
648 zz2value(den
[0], tz
);
650 for (int k
= 1; k
< dim
; ++k
) {
651 zz2value(den
[k
], tz
);
652 dpoly
fact(dim
, tz
, 1);
658 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
659 evalue
*EV
= evalue_polynomial(c
, num
.E
[i
]);
662 free_evalue_refs(num
.E
[i
]);
667 mpq_set_si(count
, 0, 1);
668 if (num
.constant
.length() == 1) {
669 zz2value(num
.constant
[0], tz
);
671 d
.div(n
, count
, sign
);
678 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
679 value_assign(acc
, c
->p
[dim
]);
680 zz2value(num
.constant
[i
], x
);
681 for (int j
= dim
-1; j
>= 0; --j
) {
682 value_multiply(acc
, acc
, x
);
683 value_addto(acc
, acc
, c
->p
[j
]);
685 value_addto(mpq_numref(count
), mpq_numref(count
), acc
);
687 mpz_set(mpq_denref(count
), c
->p
[dim
+1]);
693 evalue_set(&EV
, &count
[0]._mp_num
, &count
[0]._mp_den
);
695 free_evalue_refs(&EV
);
699 struct ienumerator_base
: enumerator_base
{
702 ienumerator_base(unsigned dim
, vertex_decomposer
*vpd
) :
703 enumerator_base(dim
,vpd
) {
704 E_vertex
= new evalue_p
[dim
];
707 virtual ~ienumerator_base() {
711 evalue
*E_num(int i
, int d
) {
712 return E_vertex
[i
+ (dim
-d
)];
721 cumulator(evalue
*factor
, evalue
*v
, dpoly_r
*r
) :
722 factor(factor
), v(v
), r(r
) {}
724 void cumulate(barvinok_options
*options
);
726 virtual void add_term(const vector
<int>& powers
, evalue
*f2
) = 0;
727 virtual ~cumulator() {}
730 void cumulator::cumulate(barvinok_options
*options
)
732 evalue cum
; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
734 evalue t
; // E_num[0] - (m-1)
738 if (options
->lookup_table
) {
740 evalue_set_si(&mone
, -1, 1);
744 evalue_copy(&cum
, factor
);
747 value_set_si(f
.d
, 1);
748 value_set_si(f
.x
.n
, 1);
752 if (!options
->lookup_table
) {
753 for (cst
= &t
; value_zero_p(cst
->d
); ) {
754 if (cst
->x
.p
->type
== fractional
)
755 cst
= &cst
->x
.p
->arr
[1];
757 cst
= &cst
->x
.p
->arr
[0];
761 for (int m
= 0; m
< r
->len
; ++m
) {
764 value_set_si(f
.d
, m
);
766 if (!options
->lookup_table
)
767 value_subtract(cst
->x
.n
, cst
->x
.n
, cst
->d
);
773 dpoly_r_term_list
& current
= r
->c
[r
->len
-1-m
];
774 dpoly_r_term_list::iterator j
;
775 for (j
= current
.begin(); j
!= current
.end(); ++j
) {
776 if ((*j
)->coeff
== 0)
778 evalue
*f2
= new evalue
;
781 zz2value((*j
)->coeff
, f2
->x
.n
);
782 zz2value(r
->denom
, f2
->d
);
785 add_term((*j
)->powers
, f2
);
788 free_evalue_refs(&f
);
789 free_evalue_refs(&t
);
790 free_evalue_refs(&cum
);
791 if (options
->lookup_table
)
792 free_evalue_refs(&mone
);
800 struct ie_cum
: public cumulator
{
801 vector
<E_poly_term
*> terms
;
803 ie_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
) : cumulator(factor
, v
, r
) {}
805 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
808 void ie_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
811 for (k
= 0; k
< terms
.size(); ++k
) {
812 if (terms
[k
]->powers
== powers
) {
813 eadd(f2
, terms
[k
]->E
);
814 free_evalue_refs(f2
);
819 if (k
>= terms
.size()) {
820 E_poly_term
*ET
= new E_poly_term
;
827 struct ienumerator
: public signed_cone_consumer
, public vertex_decomposer
,
828 public ienumerator_base
{
835 ienumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
836 vertex_decomposer(PP
, *this), ienumerator_base(dim
, this) {
837 vertex
.SetDims(1, dim
);
839 den
.SetDims(dim
, dim
);
849 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
850 void reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
851 barvinok_options
*options
);
854 void ienumerator::reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
855 barvinok_options
*options
)
857 unsigned len
= den_f
.NumRows(); // number of factors in den
858 unsigned dim
= num
.NumCols();
859 assert(num
.NumRows() == 1);
862 eadd(factor
, vE
[vert
]);
871 split_one(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
874 den_p
.SetLength(len
);
878 normalize(one
, num_s
, num_p
, den_s
, den_p
, den_r
);
884 for (int k
= 0; k
< len
; ++k
) {
887 else if (den_s
[k
] == 0)
891 reduce(factor
, num_p
, den_r
, options
);
895 pden
.SetDims(only_param
, dim
-1);
897 for (k
= 0, l
= 0; k
< len
; ++k
)
899 pden
[l
++] = den_r
[k
];
901 for (k
= 0; k
< len
; ++k
)
905 zz2value(num_s
[0], tz
);
906 dpoly
n(no_param
, tz
);
907 zz2value(den_s
[k
], tz
);
908 dpoly
D(no_param
, tz
, 1);
911 zz2value(den_s
[k
], tz
);
912 dpoly
fact(no_param
, tz
, 1);
917 // if no_param + only_param == len then all powers
918 // below will be all zero
919 if (no_param
+ only_param
== len
) {
920 if (E_num(0, dim
) != 0)
921 r
= new dpoly_r(n
, len
);
923 mpq_set_si(tcount
, 0, 1);
927 if (value_notzero_p(mpq_numref(tcount
))) {
931 value_assign(f
.x
.n
, mpq_numref(tcount
));
932 value_assign(f
.d
, mpq_denref(tcount
));
934 reduce(factor
, num_p
, pden
, options
);
935 free_evalue_refs(&f
);
940 for (k
= 0; k
< len
; ++k
) {
941 if (den_s
[k
] == 0 || den_p
[k
] == 0)
944 zz2value(den_s
[k
], tz
);
945 dpoly
pd(no_param
-1, tz
, 1);
948 for (l
= 0; l
< k
; ++l
)
949 if (den_r
[l
] == den_r
[k
])
953 r
= new dpoly_r(n
, pd
, l
, len
);
955 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
961 dpoly_r
*rc
= r
->div(D
);
964 if (E_num(0, dim
) == 0) {
965 int common
= pden
.NumRows();
966 dpoly_r_term_list
& final
= r
->c
[r
->len
-1];
972 zz2value(r
->denom
, f
.d
);
973 dpoly_r_term_list::iterator j
;
974 for (j
= final
.begin(); j
!= final
.end(); ++j
) {
975 if ((*j
)->coeff
== 0)
978 for (int k
= 0; k
< r
->dim
; ++k
) {
979 int n
= (*j
)->powers
[k
];
982 pden
.SetDims(rows
+n
, pden
.NumCols());
983 for (int l
= 0; l
< n
; ++l
)
984 pden
[rows
+l
] = den_r
[k
];
988 evalue_copy(&t
, factor
);
989 zz2value((*j
)->coeff
, f
.x
.n
);
991 reduce(&t
, num_p
, pden
, options
);
992 free_evalue_refs(&t
);
994 free_evalue_refs(&f
);
996 ie_cum
cum(factor
, E_num(0, dim
), r
);
997 cum
.cumulate(options
);
999 int common
= pden
.NumRows();
1001 for (int j
= 0; j
< cum
.terms
.size(); ++j
) {
1003 pden
.SetDims(rows
, pden
.NumCols());
1004 for (int k
= 0; k
< r
->dim
; ++k
) {
1005 int n
= cum
.terms
[j
]->powers
[k
];
1008 pden
.SetDims(rows
+n
, pden
.NumCols());
1009 for (int l
= 0; l
< n
; ++l
)
1010 pden
[rows
+l
] = den_r
[k
];
1013 reduce(cum
.terms
[j
]->E
, num_p
, pden
, options
);
1014 free_evalue_refs(cum
.terms
[j
]->E
);
1015 delete cum
.terms
[j
]->E
;
1016 delete cum
.terms
[j
];
1023 static int type_offset(enode
*p
)
1025 return p
->type
== fractional
? 1 :
1026 p
->type
== flooring
? 1 : 0;
1029 static int edegree(evalue
*e
)
1034 if (value_notzero_p(e
->d
))
1038 int i
= type_offset(p
);
1039 if (p
->size
-i
-1 > d
)
1040 d
= p
->size
- i
- 1;
1041 for (; i
< p
->size
; i
++) {
1042 int d2
= edegree(&p
->arr
[i
]);
1049 void ienumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1051 assert(sc
.det
== 1);
1052 assert(sc
.rays
.NumRows() == dim
);
1054 lattice_point(V
, sc
.rays
, vertex
[0], E_vertex
, options
);
1060 evalue_set_si(&one
, sc
.sign
, 1);
1061 reduce(&one
, vertex
, den
, options
);
1062 free_evalue_refs(&one
);
1064 for (int i
= 0; i
< dim
; ++i
)
1066 evalue_free(E_vertex
[i
]);
1069 struct bfenumerator
: public vertex_decomposer
, public bf_base
,
1070 public ienumerator_base
{
1073 bfenumerator(Polyhedron
*P
, unsigned dim
, Param_Polyhedron
*PP
) :
1074 vertex_decomposer(PP
, *this),
1075 bf_base(dim
), ienumerator_base(dim
, this) {
1083 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
1084 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1086 bfc_term_base
* new_bf_term(int len
) {
1087 bfe_term
* t
= new bfe_term(len
);
1091 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
1092 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1093 factor
= bfet
->factors
[k
];
1094 assert(factor
!= NULL
);
1095 bfet
->factors
[k
] = NULL
;
1097 emul(&mone
, factor
);
1100 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&q
, int change
) {
1101 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1102 factor
= bfet
->factors
[k
];
1103 assert(factor
!= NULL
);
1104 bfet
->factors
[k
] = NULL
;
1110 value_oppose(f
.x
.n
, mpq_numref(q
));
1112 value_assign(f
.x
.n
, mpq_numref(q
));
1113 value_assign(f
.d
, mpq_denref(q
));
1115 free_evalue_refs(&f
);
1118 virtual void set_factor(bfc_term_base
*t
, int k
, const QQ
& c
, int change
) {
1119 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1121 factor
= new evalue
;
1126 zz2value(c
.n
, f
.x
.n
);
1128 value_oppose(f
.x
.n
, f
.x
.n
);
1131 value_init(factor
->d
);
1132 evalue_copy(factor
, bfet
->factors
[k
]);
1134 free_evalue_refs(&f
);
1137 void set_factor(evalue
*f
, int change
) {
1143 virtual void insert_term(bfc_term_base
*t
, int i
) {
1144 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1145 int len
= t
->terms
.NumRows()-1; // already increased by one
1147 bfet
->factors
.resize(len
+1);
1148 for (int j
= len
; j
> i
; --j
) {
1149 bfet
->factors
[j
] = bfet
->factors
[j
-1];
1150 t
->terms
[j
] = t
->terms
[j
-1];
1152 bfet
->factors
[i
] = factor
;
1156 virtual void update_term(bfc_term_base
*t
, int i
) {
1157 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1159 eadd(factor
, bfet
->factors
[i
]);
1160 free_evalue_refs(factor
);
1164 virtual bool constant_vertex(int dim
) { return E_num(0, dim
) == 0; }
1166 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
,
1167 barvinok_options
*options
);
1170 enumerator_base
*enumerator_base::create(Polyhedron
*P
, unsigned dim
,
1171 Param_Polyhedron
*PP
,
1172 barvinok_options
*options
)
1174 enumerator_base
*eb
;
1176 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
1177 eb
= new bfenumerator(P
, dim
, PP
);
1178 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
1179 eb
= new ienumerator(P
, dim
, PP
);
1181 eb
= new enumerator(P
, dim
, PP
);
1186 struct bfe_cum
: public cumulator
{
1188 bfc_term_base
*told
;
1192 bfe_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
, bf_reducer
*bfr
,
1193 bfc_term_base
*t
, int k
, bfenumerator
*e
) :
1194 cumulator(factor
, v
, r
), told(t
), k(k
),
1198 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
1201 void bfe_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
1203 bfr
->update_powers(powers
);
1205 bfc_term_base
* t
= bfe
->find_bfc_term(bfr
->vn
, bfr
->npowers
, bfr
->nnf
);
1206 bfe
->set_factor(f2
, bfr
->l_changes
% 2);
1207 bfe
->add_term(t
, told
->terms
[k
], bfr
->l_extra_num
);
1210 void bfenumerator::cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
1211 dpoly_r
*r
, barvinok_options
*options
)
1213 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1214 bfe_cum
cum(bfet
->factors
[k
], E_num(0, bfr
->d
), r
, bfr
, t
, k
, this);
1215 cum
.cumulate(options
);
1218 void bfenumerator::base(mat_ZZ
& factors
, bfc_vec
& v
)
1220 for (int i
= 0; i
< v
.size(); ++i
) {
1221 assert(v
[i
]->terms
.NumRows() == 1);
1222 evalue
*factor
= static_cast<bfe_term
*>(v
[i
])->factors
[0];
1223 eadd(factor
, vE
[vert
]);
1228 void bfenumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1230 assert(sc
.det
== 1);
1231 assert(sc
.rays
.NumRows() == enumerator_base::dim
);
1233 bfe_term
* t
= new bfe_term(enumerator_base::dim
);
1234 vector
< bfc_term_base
* > v
;
1237 t
->factors
.resize(1);
1239 t
->terms
.SetDims(1, enumerator_base::dim
);
1240 lattice_point(V
, sc
.rays
, t
->terms
[0], E_vertex
, options
);
1242 // the elements of factors are always lexpositive
1244 int s
= setup_factors(sc
.rays
, factors
, t
, sc
.sign
);
1246 t
->factors
[0] = new evalue
;
1247 value_init(t
->factors
[0]->d
);
1248 evalue_set_si(t
->factors
[0], s
, 1);
1249 reduce(factors
, v
, options
);
1251 for (int i
= 0; i
< enumerator_base::dim
; ++i
)
1253 evalue_free(E_vertex
[i
]);
1256 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1257 barvinok_options
*options
);
1260 static evalue
* barvinok_enumerate_cst(Polyhedron
*P
, Polyhedron
* C
,
1261 struct barvinok_options
*options
)
1267 return evalue_zero();
1270 ALLOC(evalue
, eres
);
1271 value_init(eres
->d
);
1272 value_set_si(eres
->d
, 0);
1273 eres
->x
.p
= new_enode(partition
, 2, C
->Dimension
);
1274 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[0],
1275 DomainConstraintSimplify(C
, options
->MaxRays
));
1276 value_set_si(eres
->x
.p
->arr
[1].d
, 1);
1277 value_init(eres
->x
.p
->arr
[1].x
.n
);
1279 value_set_si(eres
->x
.p
->arr
[1].x
.n
, 0);
1281 barvinok_count_with_options(P
, &eres
->x
.p
->arr
[1].x
.n
, options
);
1286 static evalue
* enumerate(Polyhedron
*P
, Polyhedron
* C
,
1287 struct barvinok_options
*options
)
1290 Polyhedron
*Porig
= P
;
1291 Polyhedron
*Corig
= C
;
1292 Polyhedron
*CEq
= NULL
, *rVD
;
1294 unsigned nparam
= C
->Dimension
;
1299 value_init(factor
.d
);
1300 evalue_set_si(&factor
, 1, 1);
1303 POL_ENSURE_FACETS(P
);
1304 POL_ENSURE_VERTICES(P
);
1305 POL_ENSURE_FACETS(C
);
1306 POL_ENSURE_VERTICES(C
);
1308 if (C
->Dimension
== 0 || emptyQ(P
) || emptyQ(C
)) {
1311 CEq
= Polyhedron_Copy(CEq
);
1312 eres
= barvinok_enumerate_cst(P
, CEq
? CEq
: Polyhedron_Copy(C
), options
);
1315 evalue_backsubstitute(eres
, CP
, options
->MaxRays
);
1319 emul(&factor
, eres
);
1320 if (options
->approximation_method
== BV_APPROX_DROP
) {
1321 if (options
->polynomial_approximation
== BV_APPROX_SIGN_UPPER
)
1322 evalue_frac2polynomial(eres
, 1, options
->MaxRays
);
1323 if (options
->polynomial_approximation
== BV_APPROX_SIGN_LOWER
)
1324 evalue_frac2polynomial(eres
, -1, options
->MaxRays
);
1325 if (options
->polynomial_approximation
== BV_APPROX_SIGN_APPROX
)
1326 evalue_frac2polynomial(eres
, 0, options
->MaxRays
);
1328 reduce_evalue(eres
);
1329 free_evalue_refs(&factor
);
1337 if (Polyhedron_is_unbounded(P
, nparam
, options
->MaxRays
))
1340 if (P
->Dimension
== nparam
) {
1342 P
= Universe_Polyhedron(0);
1345 if (P
->NbEq
!= 0 || C
->NbEq
!= 0) {
1348 remove_all_equalities(&P
, &C
, &CP
, NULL
, nparam
, options
->MaxRays
);
1349 if (C
!= D
&& D
!= Corig
)
1351 if (P
!= Q
&& Q
!= Porig
)
1353 eres
= enumerate(P
, C
, options
);
1357 Polyhedron
*T
= Polyhedron_Factor(P
, nparam
, NULL
, options
->MaxRays
);
1358 if (T
|| (P
->Dimension
== nparam
+1)) {
1361 for (Q
= T
? T
: P
; Q
; Q
= Q
->next
) {
1362 Polyhedron
*next
= Q
->next
;
1366 if (Q
->Dimension
!= C
->Dimension
)
1367 QC
= Polyhedron_Project(Q
, nparam
);
1370 C
= DomainIntersection(C
, QC
, options
->MaxRays
);
1372 Polyhedron_Free(C2
);
1374 Polyhedron_Free(QC
);
1383 if (T
->Dimension
== C
->Dimension
) {
1392 eres
= barvinok_enumerate_ev_f(P
, C
, options
);
1399 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
1400 Polyhedron
*next
= Q
->next
;
1403 f
= barvinok_enumerate_ev_f(Q
, C
, options
);
1414 evalue
* barvinok_enumerate_with_options(Polyhedron
*P
, Polyhedron
* C
,
1415 struct barvinok_options
*options
)
1417 Polyhedron
*next
, *Cnext
, *C1
;
1418 Polyhedron
*Corig
= C
;
1423 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1427 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1431 C1
= Polyhedron_Project(P
, C
->Dimension
);
1432 C
= DomainIntersection(C
, C1
, options
->MaxRays
);
1433 Polyhedron_Free(C1
);
1437 if (options
->approximation_method
== BV_APPROX_BERNOULLI
||
1438 options
->summation
== BV_SUM_BERNOULLI
) {
1439 int summation
= options
->summation
;
1440 options
->summation
= BV_SUM_BERNOULLI
;
1441 eres
= barvinok_summate_unweighted(P
, C
, options
);
1442 options
->summation
= summation
;
1444 eres
= enumerate(P
, C
, options
);
1448 Corig
->next
= Cnext
;
1453 evalue
* barvinok_enumerate_ev(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1456 barvinok_options
*options
= barvinok_options_new_with_defaults();
1457 options
->MaxRays
= MaxRays
;
1458 E
= barvinok_enumerate_with_options(P
, C
, options
);
1459 barvinok_options_free(options
);
1463 evalue
*Param_Polyhedron_Enumerate(Param_Polyhedron
*PP
, Polyhedron
*P
,
1465 struct barvinok_options
*options
)
1469 unsigned nparam
= C
->Dimension
;
1470 unsigned dim
= P
->Dimension
- nparam
;
1473 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
1474 evalue_section
*s
= new evalue_section
[nd
];
1476 enumerator_base
*et
= NULL
;
1481 et
= enumerator_base::create(P
, dim
, PP
, options
);
1483 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
1484 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
1487 s
[i
].E
= evalue_zero();
1490 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
1493 et
->decompose_at(V
, _i
, options
);
1494 } catch (OrthogonalException
&e
) {
1495 FORALL_REDUCED_DOMAIN_RESET
;
1496 for (; i
>= 0; --i
) {
1497 evalue_free(s
[i
].E
);
1498 Domain_Free(s
[i
].D
);
1502 eadd(et
->vE
[_i
] , s
[i
].E
);
1503 END_FORALL_PVertex_in_ParamPolyhedron
;
1504 evalue_range_reduction_in_domain(s
[i
].E
, rVD
);
1505 END_FORALL_REDUCED_DOMAIN
1506 Polyhedron_Free(TC
);
1509 eres
= evalue_from_section_array(s
, nd
);
1515 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1516 barvinok_options
*options
)
1518 unsigned nparam
= C
->Dimension
;
1519 bool do_scale
= options
->approximation_method
== BV_APPROX_SCALE
;
1521 if (options
->summation
== BV_SUM_EULER
)
1522 return barvinok_summate_unweighted(P
, C
, options
);
1524 if (options
->approximation_method
== BV_APPROX_VOLUME
)
1525 return Param_Polyhedron_Volume(P
, C
, options
);
1527 if (P
->Dimension
- nparam
== 1 && !do_scale
)
1528 return ParamLine_Length(P
, C
, options
);
1530 Param_Polyhedron
*PP
= NULL
;
1534 eres
= scale_bound(P
, C
, options
);
1539 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
1542 eres
= scale(PP
, P
, C
, options
);
1544 eres
= Param_Polyhedron_Enumerate(PP
, P
, C
, options
);
1547 Param_Polyhedron_Free(PP
);
1552 Enumeration
* barvinok_enumerate(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1554 evalue
*EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
1556 return partition2enumeration(EP
);
1559 evalue
* barvinok_enumerate_union(Polyhedron
*D
, Polyhedron
* C
, unsigned MaxRays
)
1562 gen_fun
*gf
= barvinok_enumerate_union_series(D
, C
, MaxRays
);