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"
21 #include "decomposer.h"
23 #include "lattice_point.h"
24 #include "reduce_domain.h"
25 #include "remove_equalities.h"
28 #include "bernoulli.h"
29 #include "param_util.h"
40 using std::ostringstream
;
42 #define ALLOC(t,p) p = (t*)malloc(sizeof(*p))
55 coeff
= Matrix_Alloc(d
+1, d
+1+1);
56 value_set_si(coeff
->p
[0][0], 1);
57 value_set_si(coeff
->p
[0][d
+1], 1);
58 for (int i
= 1; i
<= d
; ++i
) {
59 value_multiply(coeff
->p
[i
][0], coeff
->p
[i
-1][0], d0
);
60 Vector_Combine(coeff
->p
[i
-1], coeff
->p
[i
-1]+1, coeff
->p
[i
]+1,
62 value_set_si(coeff
->p
[i
][d
+1], i
);
63 value_multiply(coeff
->p
[i
][d
+1], coeff
->p
[i
][d
+1], coeff
->p
[i
-1][d
+1]);
64 value_decrement(d0
, d0
);
69 void div(dpoly
& d
, Vector
*count
, int sign
) {
70 int len
= coeff
->NbRows
;
71 Matrix
* c
= Matrix_Alloc(coeff
->NbRows
, coeff
->NbColumns
);
74 for (int i
= 0; i
< len
; ++i
) {
75 Vector_Copy(coeff
->p
[i
], c
->p
[i
], len
+1);
76 for (int j
= 1; j
<= i
; ++j
) {
77 value_multiply(tmp
, d
.coeff
->p
[j
], c
->p
[i
][len
]);
78 value_oppose(tmp
, tmp
);
79 Vector_Combine(c
->p
[i
], c
->p
[i
-j
], c
->p
[i
],
80 c
->p
[i
-j
][len
], tmp
, len
);
81 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], c
->p
[i
-j
][len
]);
83 value_multiply(c
->p
[i
][len
], c
->p
[i
][len
], d
.coeff
->p
[0]);
86 value_set_si(tmp
, -1);
87 Vector_Scale(c
->p
[len
-1], count
->p
, tmp
, len
);
88 value_assign(count
->p
[len
], c
->p
[len
-1][len
]);
90 Vector_Copy(c
->p
[len
-1], count
->p
, len
+1);
91 Vector_Normalize(count
->p
, len
+1);
99 * Searches for a vector that is not orthogonal to any
100 * of the rays in rays.
102 static void nonorthog(mat_ZZ
& rays
, vec_ZZ
& lambda
)
104 int dim
= rays
.NumCols();
106 lambda
.SetLength(dim
);
110 for (int i
= 2; !found
&& i
<= 50*dim
; i
+=4) {
111 for (int j
= 0; j
< MAX_TRY
; ++j
) {
112 for (int k
= 0; k
< dim
; ++k
) {
113 int r
= random_int(i
)+2;
114 int v
= (2*(r
%2)-1) * (r
>> 1);
118 for (; k
< rays
.NumRows(); ++k
)
119 if (lambda
* rays
[k
] == 0)
121 if (k
== rays
.NumRows()) {
130 static void add_rays(mat_ZZ
& rays
, Polyhedron
*i
, int *r
, int nvar
= -1,
133 unsigned dim
= i
->Dimension
;
136 for (int k
= 0; k
< i
->NbRays
; ++k
) {
137 if (!value_zero_p(i
->Ray
[k
][dim
+1]))
139 if (!all
&& nvar
!= dim
&& First_Non_Zero(i
->Ray
[k
]+1, nvar
) == -1)
141 values2zz(i
->Ray
[k
]+1, rays
[(*r
)++], nvar
);
145 static void mask_r(Matrix
*f
, int nr
, Vector
*lcm
, int p
, Vector
*val
, evalue
*ev
)
147 unsigned nparam
= lcm
->Size
;
150 Vector
* prod
= Vector_Alloc(f
->NbRows
);
151 Matrix_Vector_Product(f
, val
->p
, prod
->p
);
153 for (int i
= 0; i
< nr
; ++i
) {
154 value_modulus(prod
->p
[i
], prod
->p
[i
], f
->p
[i
][nparam
+1]);
155 isint
&= value_zero_p(prod
->p
[i
]);
157 value_set_si(ev
->d
, 1);
159 value_set_si(ev
->x
.n
, isint
);
166 if (value_one_p(lcm
->p
[p
]))
167 mask_r(f
, nr
, lcm
, p
+1, val
, ev
);
169 value_assign(tmp
, lcm
->p
[p
]);
170 value_set_si(ev
->d
, 0);
171 ev
->x
.p
= new_enode(periodic
, VALUE_TO_INT(tmp
), p
+1);
173 value_decrement(tmp
, tmp
);
174 value_assign(val
->p
[p
], tmp
);
175 mask_r(f
, nr
, lcm
, p
+1, val
, &ev
->x
.p
->arr
[VALUE_TO_INT(tmp
)]);
176 } while (value_pos_p(tmp
));
181 static void mask_fractional(Matrix
*f
, evalue
*factor
)
183 int nr
= f
->NbRows
, nc
= f
->NbColumns
;
186 for (n
= 0; n
< nr
&& value_notzero_p(f
->p
[n
][nc
-1]); ++n
)
187 if (value_notone_p(f
->p
[n
][nc
-1]) &&
188 value_notmone_p(f
->p
[n
][nc
-1]))
202 value_set_si(EV
.x
.n
, 1);
204 for (n
= 0; n
< nr
; ++n
) {
205 value_assign(m
, f
->p
[n
][nc
-1]);
206 if (value_one_p(m
) || value_mone_p(m
))
209 int j
= normal_mod(f
->p
[n
], nc
-1, &m
);
211 free_evalue_refs(factor
);
212 value_init(factor
->d
);
213 evalue_set_si(factor
, 0, 1);
217 values2zz(f
->p
[n
], row
, nc
-1);
220 if (j
< (nc
-1)-1 && row
[j
] > g
/2) {
221 for (int k
= j
; k
< (nc
-1); ++k
)
227 value_set_si(EP
.d
, 0);
228 EP
.x
.p
= new_enode(relation
, 2, 0);
229 value_clear(EP
.x
.p
->arr
[1].d
);
230 EP
.x
.p
->arr
[1] = *factor
;
231 evalue
*ev
= &EP
.x
.p
->arr
[0];
232 value_set_si(ev
->d
, 0);
233 ev
->x
.p
= new_enode(fractional
, 3, -1);
234 evalue_set_si(&ev
->x
.p
->arr
[1], 0, 1);
235 evalue_set_si(&ev
->x
.p
->arr
[2], 1, 1);
236 evalue
*E
= multi_monom(row
);
237 value_assign(EV
.d
, m
);
239 value_clear(ev
->x
.p
->arr
[0].d
);
240 ev
->x
.p
->arr
[0] = *E
;
246 free_evalue_refs(&EV
);
252 static void mask_table(Matrix
*f
, evalue
*factor
)
254 int nr
= f
->NbRows
, nc
= f
->NbColumns
;
257 for (n
= 0; n
< nr
&& value_notzero_p(f
->p
[n
][nc
-1]); ++n
)
258 if (value_notone_p(f
->p
[n
][nc
-1]) &&
259 value_notmone_p(f
->p
[n
][nc
-1]))
267 unsigned np
= nc
- 2;
268 Vector
*lcm
= Vector_Alloc(np
);
269 Vector
*val
= Vector_Alloc(nc
);
270 Vector_Set(val
->p
, 0, nc
);
271 value_set_si(val
->p
[np
], 1);
272 Vector_Set(lcm
->p
, 1, np
);
273 for (n
= 0; n
< nr
; ++n
) {
274 if (value_one_p(f
->p
[n
][nc
-1]) ||
275 value_mone_p(f
->p
[n
][nc
-1]))
277 for (int j
= 0; j
< np
; ++j
)
278 if (value_notzero_p(f
->p
[n
][j
])) {
279 value_gcd(tmp
, f
->p
[n
][j
], f
->p
[n
][nc
-1]);
280 value_division(tmp
, f
->p
[n
][nc
-1], tmp
);
281 value_lcm(lcm
->p
[j
], tmp
, lcm
->p
[j
]);
286 mask_r(f
, nr
, lcm
, 0, val
, &EP
);
291 free_evalue_refs(&EP
);
294 static void mask(Matrix
*f
, evalue
*factor
, barvinok_options
*options
)
296 if (options
->lookup_table
)
297 mask_table(f
, factor
);
299 mask_fractional(f
, factor
);
302 struct bfe_term
: public bfc_term_base
{
303 vector
<evalue
*> factors
;
305 bfe_term(int len
) : bfc_term_base(len
) {
309 for (int i
= 0; i
< factors
.size(); ++i
) {
312 free_evalue_refs(factors
[i
]);
318 static void print_int_vector(int *v
, int len
, const char *name
)
320 cerr
<< name
<< endl
;
321 for (int j
= 0; j
< len
; ++j
) {
327 static void print_bfc_terms(mat_ZZ
& factors
, bfc_vec
& v
)
330 cerr
<< "factors" << endl
;
331 cerr
<< factors
<< endl
;
332 for (int i
= 0; i
< v
.size(); ++i
) {
333 cerr
<< "term: " << i
<< endl
;
334 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
335 cerr
<< "terms" << endl
;
336 cerr
<< v
[i
]->terms
<< endl
;
337 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
338 cerr
<< bfct
->c
<< endl
;
342 static void print_bfe_terms(mat_ZZ
& factors
, bfc_vec
& v
)
345 cerr
<< "factors" << endl
;
346 cerr
<< factors
<< endl
;
347 for (int i
= 0; i
< v
.size(); ++i
) {
348 cerr
<< "term: " << i
<< endl
;
349 print_int_vector(v
[i
]->powers
, factors
.NumRows(), "powers");
350 cerr
<< "terms" << endl
;
351 cerr
<< v
[i
]->terms
<< endl
;
352 bfe_term
* bfet
= static_cast<bfe_term
*>(v
[i
]);
353 for (int j
= 0; j
< v
[i
]->terms
.NumRows(); ++j
) {
354 const char * test
[] = {"a", "b"};
355 print_evalue(stderr
, bfet
->factors
[j
], test
);
356 fprintf(stderr
, "\n");
361 struct bfcounter
: public bfcounter_base
{
365 bfcounter(unsigned dim
) : bfcounter_base(dim
) {
374 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
375 virtual void get_count(Value
*result
) {
376 assert(value_one_p(&count
[0]._mp_den
));
377 value_assign(*result
, &count
[0]._mp_num
);
381 void bfcounter::base(mat_ZZ
& factors
, bfc_vec
& v
)
383 unsigned nf
= factors
.NumRows();
385 for (int i
= 0; i
< v
.size(); ++i
) {
386 bfc_term
* bfct
= static_cast<bfc_term
*>(v
[i
]);
388 // factor is always positive, so we always
390 for (int k
= 0; k
< nf
; ++k
)
391 total_power
+= v
[i
]->powers
[k
];
394 for (j
= 0; j
< nf
; ++j
)
395 if (v
[i
]->powers
[j
] > 0)
398 zz2value(factors
[j
][0], tz
);
399 dpoly
D(total_power
, tz
, 1);
400 for (int k
= 1; k
< v
[i
]->powers
[j
]; ++k
) {
401 zz2value(factors
[j
][0], tz
);
402 dpoly
fact(total_power
, tz
, 1);
406 for (int k
= 0; k
< v
[i
]->powers
[j
]; ++k
) {
407 zz2value(factors
[j
][0], tz
);
408 dpoly
fact(total_power
, tz
, 1);
412 for (int k
= 0; k
< v
[i
]->terms
.NumRows(); ++k
) {
413 zz2value(v
[i
]->terms
[k
][0], tz
);
414 dpoly
n(total_power
, tz
);
415 mpq_set_si(tcount
, 0, 1);
418 bfct
->c
[k
].n
= -bfct
->c
[k
].n
;
419 zz2value(bfct
->c
[k
].n
, tn
);
420 zz2value(bfct
->c
[k
].d
, td
);
422 mpz_mul(mpq_numref(tcount
), mpq_numref(tcount
), tn
);
423 mpz_mul(mpq_denref(tcount
), mpq_denref(tcount
), td
);
424 mpq_canonicalize(tcount
);
425 mpq_add(count
, count
, tcount
);
432 /* Check whether the polyhedron is unbounded and if so,
433 * check whether it has any (and therefore an infinite number of)
435 * If one of the vertices is integer, then we are done.
436 * Otherwise, transform the polyhedron such that one of the rays
437 * is the first unit vector and cut it off at a height that ensures
438 * that if the whole polyhedron has any points, then the remaining part
439 * has integer points. In particular we add the largest coefficient
440 * of a ray to the highest vertex (rounded up).
442 static bool Polyhedron_is_infinite(Polyhedron
*P
, Value
* result
,
443 barvinok_options
*options
)
455 for (; r
< P
->NbRays
; ++r
)
456 if (value_zero_p(P
->Ray
[r
][P
->Dimension
+1]))
458 if (P
->NbBid
== 0 && r
== P
->NbRays
)
461 if (options
->count_sample_infinite
) {
464 sample
= Polyhedron_Sample(P
, options
);
466 value_set_si(*result
, 0);
468 value_set_si(*result
, -1);
474 for (int i
= 0; i
< P
->NbRays
; ++i
)
475 if (value_one_p(P
->Ray
[i
][1+P
->Dimension
])) {
476 value_set_si(*result
, -1);
481 M
= Matrix_Alloc(P
->Dimension
+1, P
->Dimension
+1);
482 Vector_Gcd(P
->Ray
[r
]+1, P
->Dimension
, &g
);
483 Vector_AntiScale(P
->Ray
[r
]+1, M
->p
[0], g
, P
->Dimension
+1);
484 int ok
= unimodular_complete(M
, 1);
486 value_set_si(M
->p
[P
->Dimension
][P
->Dimension
], 1);
489 P
= Polyhedron_Preimage(P
, M2
, 0);
497 value_set_si(size
, 0);
499 for (int i
= 0; i
< P
->NbBid
; ++i
) {
500 value_absolute(tmp
, P
->Ray
[i
][1]);
501 if (value_gt(tmp
, size
))
502 value_assign(size
, tmp
);
504 for (int i
= P
->NbBid
; i
< P
->NbRays
; ++i
) {
505 if (value_zero_p(P
->Ray
[i
][P
->Dimension
+1])) {
506 if (value_gt(P
->Ray
[i
][1], size
))
507 value_assign(size
, P
->Ray
[i
][1]);
510 mpz_cdiv_q(tmp
, P
->Ray
[i
][1], P
->Ray
[i
][P
->Dimension
+1]);
511 if (first
|| value_gt(tmp
, offset
)) {
512 value_assign(offset
, tmp
);
516 value_addto(offset
, offset
, size
);
520 v
= Vector_Alloc(P
->Dimension
+2);
521 value_set_si(v
->p
[0], 1);
522 value_set_si(v
->p
[1], -1);
523 value_assign(v
->p
[1+P
->Dimension
], offset
);
524 R
= AddConstraints(v
->p
, 1, P
, options
->MaxRays
);
532 barvinok_count_with_options(P
, &c
, options
);
535 value_set_si(*result
, 0);
537 value_set_si(*result
, -1);
543 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
544 barvinok_options
*options
);
546 void barvinok_count_with_options(Polyhedron
*P
, Value
* result
,
547 struct barvinok_options
*options
)
552 bool infinite
= false;
556 "barvinok_count: input is a union; only first polyhedron is counted\n");
559 value_set_si(*result
, 0);
565 P
= remove_equalities(P
, options
->MaxRays
);
566 P
= DomainConstraintSimplify(P
, options
->MaxRays
);
570 } while (!emptyQ(P
) && P
->NbEq
!= 0);
573 value_set_si(*result
, 0);
578 if (Polyhedron_is_infinite(P
, result
, options
)) {
583 if (P
->Dimension
== 0) {
584 /* Test whether the constraints are satisfied */
585 POL_ENSURE_VERTICES(P
);
586 value_set_si(*result
, !emptyQ(P
));
591 Q
= Polyhedron_Factor(P
, 0, NULL
, options
->MaxRays
);
599 barvinok_count_f(P
, result
, options
);
600 if (value_neg_p(*result
))
602 if (Q
&& P
->next
&& value_notzero_p(*result
)) {
606 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
607 barvinok_count_f(Q
, &factor
, options
);
608 if (value_neg_p(factor
)) {
611 } else if (Q
->next
&& value_zero_p(factor
)) {
612 value_set_si(*result
, 0);
615 value_multiply(*result
, *result
, factor
);
624 value_set_si(*result
, -1);
627 void barvinok_count(Polyhedron
*P
, Value
* result
, unsigned NbMaxCons
)
629 barvinok_options
*options
= barvinok_options_new_with_defaults();
630 options
->MaxRays
= NbMaxCons
;
631 barvinok_count_with_options(P
, result
, options
);
632 barvinok_options_free(options
);
635 static void barvinok_count_f(Polyhedron
*P
, Value
* result
,
636 barvinok_options
*options
)
639 value_set_si(*result
, 0);
643 if (P
->Dimension
== 1)
644 return Line_Length(P
, result
);
646 int c
= P
->NbConstraints
;
647 POL_ENSURE_FACETS(P
);
648 if (c
!= P
->NbConstraints
|| P
->NbEq
!= 0) {
649 Polyhedron
*next
= P
->next
;
651 barvinok_count_with_options(P
, result
, options
);
656 POL_ENSURE_VERTICES(P
);
658 if (Polyhedron_is_infinite(P
, result
, options
))
662 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
663 cnt
= new bfcounter(P
->Dimension
);
664 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
665 cnt
= new icounter(P
->Dimension
);
666 else if (options
->incremental_specialization
== BV_SPECIALIZATION_TODD
)
667 cnt
= new tcounter(P
->Dimension
, options
->max_index
);
669 cnt
= new counter(P
->Dimension
, options
->max_index
);
670 cnt
->start(P
, options
);
672 cnt
->get_count(result
);
676 static void uni_polynom(int param
, Vector
*c
, evalue
*EP
)
678 unsigned dim
= c
->Size
-2;
680 value_set_si(EP
->d
,0);
681 EP
->x
.p
= new_enode(polynomial
, dim
+1, param
+1);
682 for (int j
= 0; j
<= dim
; ++j
)
683 evalue_set(&EP
->x
.p
->arr
[j
], c
->p
[j
], c
->p
[dim
+1]);
686 typedef evalue
* evalue_p
;
688 struct enumerator_base
{
692 vertex_decomposer
*vpd
;
694 enumerator_base(unsigned dim
, vertex_decomposer
*vpd
)
699 vE
= new evalue_p
[vpd
->nbV
];
700 for (int j
= 0; j
< vpd
->nbV
; ++j
)
704 evalue_set_si(&mone
, -1, 1);
707 void decompose_at(Param_Vertices
*V
, int _i
, barvinok_options
*options
) {
711 value_init(vE
[_i
]->d
);
712 evalue_set_si(vE
[_i
], 0, 1);
714 vpd
->decompose_at_vertex(V
, _i
, options
);
717 virtual ~enumerator_base() {
718 for (int j
= 0; j
< vpd
->nbV
; ++j
)
720 free_evalue_refs(vE
[j
]);
725 free_evalue_refs(&mone
);
728 static enumerator_base
*create(Polyhedron
*P
, unsigned dim
, unsigned nbV
,
729 barvinok_options
*options
);
732 struct enumerator
: public signed_cone_consumer
, public vertex_decomposer
,
733 public enumerator_base
{
741 enumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) :
742 vertex_decomposer(P
, nbV
, *this), enumerator_base(dim
, this) {
745 randomvector(P
, lambda
, dim
);
747 c
= Vector_Alloc(dim
+2);
759 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
762 void enumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
766 assert(sc
.rays
.NumRows() == dim
);
767 for (int k
= 0; k
< dim
; ++k
) {
768 if (lambda
* sc
.rays
[k
] == 0)
772 lattice_point(V
, sc
.rays
, lambda
, &num
, sc
.det
, options
);
773 den
= sc
.rays
* lambda
;
778 zz2value(den
[0], tz
);
780 for (int k
= 1; k
< dim
; ++k
) {
781 zz2value(den
[k
], tz
);
782 dpoly
fact(dim
, tz
, 1);
788 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
789 evalue
*EV
= evalue_polynomial(c
, num
.E
[i
]);
792 free_evalue_refs(num
.E
[i
]);
797 mpq_set_si(count
, 0, 1);
798 if (num
.constant
.length() == 1) {
799 zz2value(num
.constant
[0], tz
);
801 d
.div(n
, count
, sign
);
808 for (unsigned long i
= 0; i
< sc
.det
; ++i
) {
809 value_assign(acc
, c
->p
[dim
]);
810 zz2value(num
.constant
[i
], x
);
811 for (int j
= dim
-1; j
>= 0; --j
) {
812 value_multiply(acc
, acc
, x
);
813 value_addto(acc
, acc
, c
->p
[j
]);
815 value_addto(mpq_numref(count
), mpq_numref(count
), acc
);
817 mpz_set(mpq_denref(count
), c
->p
[dim
+1]);
823 evalue_set(&EV
, &count
[0]._mp_num
, &count
[0]._mp_den
);
825 free_evalue_refs(&EV
);
829 struct ienumerator_base
: enumerator_base
{
832 ienumerator_base(unsigned dim
, vertex_decomposer
*vpd
) :
833 enumerator_base(dim
,vpd
) {
834 E_vertex
= new evalue_p
[dim
];
837 virtual ~ienumerator_base() {
841 evalue
*E_num(int i
, int d
) {
842 return E_vertex
[i
+ (dim
-d
)];
851 cumulator(evalue
*factor
, evalue
*v
, dpoly_r
*r
) :
852 factor(factor
), v(v
), r(r
) {}
854 void cumulate(barvinok_options
*options
);
856 virtual void add_term(const vector
<int>& powers
, evalue
*f2
) = 0;
857 virtual ~cumulator() {}
860 void cumulator::cumulate(barvinok_options
*options
)
862 evalue cum
; // factor * 1 * E_num[0]/1 * (E_num[0]-1)/2 *...
864 evalue t
; // E_num[0] - (m-1)
868 if (options
->lookup_table
) {
870 evalue_set_si(&mone
, -1, 1);
874 evalue_copy(&cum
, factor
);
877 value_set_si(f
.d
, 1);
878 value_set_si(f
.x
.n
, 1);
882 if (!options
->lookup_table
) {
883 for (cst
= &t
; value_zero_p(cst
->d
); ) {
884 if (cst
->x
.p
->type
== fractional
)
885 cst
= &cst
->x
.p
->arr
[1];
887 cst
= &cst
->x
.p
->arr
[0];
891 for (int m
= 0; m
< r
->len
; ++m
) {
894 value_set_si(f
.d
, m
);
896 if (!options
->lookup_table
)
897 value_subtract(cst
->x
.n
, cst
->x
.n
, cst
->d
);
903 dpoly_r_term_list
& current
= r
->c
[r
->len
-1-m
];
904 dpoly_r_term_list::iterator j
;
905 for (j
= current
.begin(); j
!= current
.end(); ++j
) {
906 if ((*j
)->coeff
== 0)
908 evalue
*f2
= new evalue
;
911 zz2value((*j
)->coeff
, f2
->x
.n
);
912 zz2value(r
->denom
, f2
->d
);
915 add_term((*j
)->powers
, f2
);
918 free_evalue_refs(&f
);
919 free_evalue_refs(&t
);
920 free_evalue_refs(&cum
);
921 if (options
->lookup_table
)
922 free_evalue_refs(&mone
);
930 struct ie_cum
: public cumulator
{
931 vector
<E_poly_term
*> terms
;
933 ie_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
) : cumulator(factor
, v
, r
) {}
935 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
938 void ie_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
941 for (k
= 0; k
< terms
.size(); ++k
) {
942 if (terms
[k
]->powers
== powers
) {
943 eadd(f2
, terms
[k
]->E
);
944 free_evalue_refs(f2
);
949 if (k
>= terms
.size()) {
950 E_poly_term
*ET
= new E_poly_term
;
957 struct ienumerator
: public signed_cone_consumer
, public vertex_decomposer
,
958 public ienumerator_base
{
965 ienumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) :
966 vertex_decomposer(P
, nbV
, *this), ienumerator_base(dim
, this) {
967 vertex
.SetDims(1, dim
);
969 den
.SetDims(dim
, dim
);
979 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
980 void reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
981 barvinok_options
*options
);
984 void ienumerator::reduce(evalue
*factor
, const mat_ZZ
& num
, const mat_ZZ
& den_f
,
985 barvinok_options
*options
)
987 unsigned len
= den_f
.NumRows(); // number of factors in den
988 unsigned dim
= num
.NumCols();
989 assert(num
.NumRows() == 1);
992 eadd(factor
, vE
[vert
]);
1001 split_one(num
, num_s
, num_p
, den_f
, den_s
, den_r
);
1004 den_p
.SetLength(len
);
1008 normalize(one
, num_s
, num_p
, den_s
, den_p
, den_r
);
1010 emul(&mone
, factor
);
1014 for (int k
= 0; k
< len
; ++k
) {
1017 else if (den_s
[k
] == 0)
1020 if (no_param
== 0) {
1021 reduce(factor
, num_p
, den_r
, options
);
1025 pden
.SetDims(only_param
, dim
-1);
1027 for (k
= 0, l
= 0; k
< len
; ++k
)
1029 pden
[l
++] = den_r
[k
];
1031 for (k
= 0; k
< len
; ++k
)
1035 zz2value(num_s
[0], tz
);
1036 dpoly
n(no_param
, tz
);
1037 zz2value(den_s
[k
], tz
);
1038 dpoly
D(no_param
, tz
, 1);
1039 for ( ; ++k
< len
; )
1040 if (den_p
[k
] == 0) {
1041 zz2value(den_s
[k
], tz
);
1042 dpoly
fact(no_param
, tz
, 1);
1047 // if no_param + only_param == len then all powers
1048 // below will be all zero
1049 if (no_param
+ only_param
== len
) {
1050 if (E_num(0, dim
) != 0)
1051 r
= new dpoly_r(n
, len
);
1053 mpq_set_si(tcount
, 0, 1);
1055 n
.div(D
, tcount
, 1);
1057 if (value_notzero_p(mpq_numref(tcount
))) {
1061 value_assign(f
.x
.n
, mpq_numref(tcount
));
1062 value_assign(f
.d
, mpq_denref(tcount
));
1064 reduce(factor
, num_p
, pden
, options
);
1065 free_evalue_refs(&f
);
1070 for (k
= 0; k
< len
; ++k
) {
1071 if (den_s
[k
] == 0 || den_p
[k
] == 0)
1074 zz2value(den_s
[k
], tz
);
1075 dpoly
pd(no_param
-1, tz
, 1);
1078 for (l
= 0; l
< k
; ++l
)
1079 if (den_r
[l
] == den_r
[k
])
1083 r
= new dpoly_r(n
, pd
, l
, len
);
1085 dpoly_r
*nr
= new dpoly_r(r
, pd
, l
, len
);
1091 dpoly_r
*rc
= r
->div(D
);
1094 if (E_num(0, dim
) == 0) {
1095 int common
= pden
.NumRows();
1096 dpoly_r_term_list
& final
= r
->c
[r
->len
-1];
1102 zz2value(r
->denom
, f
.d
);
1103 dpoly_r_term_list::iterator j
;
1104 for (j
= final
.begin(); j
!= final
.end(); ++j
) {
1105 if ((*j
)->coeff
== 0)
1108 for (int k
= 0; k
< r
->dim
; ++k
) {
1109 int n
= (*j
)->powers
[k
];
1112 pden
.SetDims(rows
+n
, pden
.NumCols());
1113 for (int l
= 0; l
< n
; ++l
)
1114 pden
[rows
+l
] = den_r
[k
];
1118 evalue_copy(&t
, factor
);
1119 zz2value((*j
)->coeff
, f
.x
.n
);
1121 reduce(&t
, num_p
, pden
, options
);
1122 free_evalue_refs(&t
);
1124 free_evalue_refs(&f
);
1126 ie_cum
cum(factor
, E_num(0, dim
), r
);
1127 cum
.cumulate(options
);
1129 int common
= pden
.NumRows();
1131 for (int j
= 0; j
< cum
.terms
.size(); ++j
) {
1133 pden
.SetDims(rows
, pden
.NumCols());
1134 for (int k
= 0; k
< r
->dim
; ++k
) {
1135 int n
= cum
.terms
[j
]->powers
[k
];
1138 pden
.SetDims(rows
+n
, pden
.NumCols());
1139 for (int l
= 0; l
< n
; ++l
)
1140 pden
[rows
+l
] = den_r
[k
];
1143 reduce(cum
.terms
[j
]->E
, num_p
, pden
, options
);
1144 free_evalue_refs(cum
.terms
[j
]->E
);
1145 delete cum
.terms
[j
]->E
;
1146 delete cum
.terms
[j
];
1153 static int type_offset(enode
*p
)
1155 return p
->type
== fractional
? 1 :
1156 p
->type
== flooring
? 1 : 0;
1159 static int edegree(evalue
*e
)
1164 if (value_notzero_p(e
->d
))
1168 int i
= type_offset(p
);
1169 if (p
->size
-i
-1 > d
)
1170 d
= p
->size
- i
- 1;
1171 for (; i
< p
->size
; i
++) {
1172 int d2
= edegree(&p
->arr
[i
]);
1179 void ienumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1181 assert(sc
.det
== 1);
1182 assert(sc
.rays
.NumRows() == dim
);
1184 lattice_point(V
, sc
.rays
, vertex
[0], E_vertex
, options
);
1190 evalue_set_si(&one
, sc
.sign
, 1);
1191 reduce(&one
, vertex
, den
, options
);
1192 free_evalue_refs(&one
);
1194 for (int i
= 0; i
< dim
; ++i
)
1196 free_evalue_refs(E_vertex
[i
]);
1201 struct bfenumerator
: public vertex_decomposer
, public bf_base
,
1202 public ienumerator_base
{
1205 bfenumerator(Polyhedron
*P
, unsigned dim
, unsigned nbV
) :
1206 vertex_decomposer(P
, nbV
, *this),
1207 bf_base(dim
), ienumerator_base(dim
, this) {
1215 virtual void handle(const signed_cone
& sc
, barvinok_options
*options
);
1216 virtual void base(mat_ZZ
& factors
, bfc_vec
& v
);
1218 bfc_term_base
* new_bf_term(int len
) {
1219 bfe_term
* t
= new bfe_term(len
);
1223 virtual void set_factor(bfc_term_base
*t
, int k
, int change
) {
1224 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1225 factor
= bfet
->factors
[k
];
1226 assert(factor
!= NULL
);
1227 bfet
->factors
[k
] = NULL
;
1229 emul(&mone
, factor
);
1232 virtual void set_factor(bfc_term_base
*t
, int k
, mpq_t
&q
, int change
) {
1233 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1234 factor
= bfet
->factors
[k
];
1235 assert(factor
!= NULL
);
1236 bfet
->factors
[k
] = NULL
;
1242 value_oppose(f
.x
.n
, mpq_numref(q
));
1244 value_assign(f
.x
.n
, mpq_numref(q
));
1245 value_assign(f
.d
, mpq_denref(q
));
1247 free_evalue_refs(&f
);
1250 virtual void set_factor(bfc_term_base
*t
, int k
, const QQ
& c
, int change
) {
1251 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1253 factor
= new evalue
;
1258 zz2value(c
.n
, f
.x
.n
);
1260 value_oppose(f
.x
.n
, f
.x
.n
);
1263 value_init(factor
->d
);
1264 evalue_copy(factor
, bfet
->factors
[k
]);
1266 free_evalue_refs(&f
);
1269 void set_factor(evalue
*f
, int change
) {
1275 virtual void insert_term(bfc_term_base
*t
, int i
) {
1276 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1277 int len
= t
->terms
.NumRows()-1; // already increased by one
1279 bfet
->factors
.resize(len
+1);
1280 for (int j
= len
; j
> i
; --j
) {
1281 bfet
->factors
[j
] = bfet
->factors
[j
-1];
1282 t
->terms
[j
] = t
->terms
[j
-1];
1284 bfet
->factors
[i
] = factor
;
1288 virtual void update_term(bfc_term_base
*t
, int i
) {
1289 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1291 eadd(factor
, bfet
->factors
[i
]);
1292 free_evalue_refs(factor
);
1296 virtual bool constant_vertex(int dim
) { return E_num(0, dim
) == 0; }
1298 virtual void cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
, dpoly_r
*r
,
1299 barvinok_options
*options
);
1302 enumerator_base
*enumerator_base::create(Polyhedron
*P
, unsigned dim
, unsigned nbV
,
1303 barvinok_options
*options
)
1305 enumerator_base
*eb
;
1307 if (options
->incremental_specialization
== BV_SPECIALIZATION_BF
)
1308 eb
= new bfenumerator(P
, dim
, nbV
);
1309 else if (options
->incremental_specialization
== BV_SPECIALIZATION_DF
)
1310 eb
= new ienumerator(P
, dim
, nbV
);
1312 eb
= new enumerator(P
, dim
, nbV
);
1317 struct bfe_cum
: public cumulator
{
1319 bfc_term_base
*told
;
1323 bfe_cum(evalue
*factor
, evalue
*v
, dpoly_r
*r
, bf_reducer
*bfr
,
1324 bfc_term_base
*t
, int k
, bfenumerator
*e
) :
1325 cumulator(factor
, v
, r
), told(t
), k(k
),
1329 virtual void add_term(const vector
<int>& powers
, evalue
*f2
);
1332 void bfe_cum::add_term(const vector
<int>& powers
, evalue
*f2
)
1334 bfr
->update_powers(powers
);
1336 bfc_term_base
* t
= bfe
->find_bfc_term(bfr
->vn
, bfr
->npowers
, bfr
->nnf
);
1337 bfe
->set_factor(f2
, bfr
->l_changes
% 2);
1338 bfe
->add_term(t
, told
->terms
[k
], bfr
->l_extra_num
);
1341 void bfenumerator::cum(bf_reducer
*bfr
, bfc_term_base
*t
, int k
,
1342 dpoly_r
*r
, barvinok_options
*options
)
1344 bfe_term
* bfet
= static_cast<bfe_term
*>(t
);
1345 bfe_cum
cum(bfet
->factors
[k
], E_num(0, bfr
->d
), r
, bfr
, t
, k
, this);
1346 cum
.cumulate(options
);
1349 void bfenumerator::base(mat_ZZ
& factors
, bfc_vec
& v
)
1351 for (int i
= 0; i
< v
.size(); ++i
) {
1352 assert(v
[i
]->terms
.NumRows() == 1);
1353 evalue
*factor
= static_cast<bfe_term
*>(v
[i
])->factors
[0];
1354 eadd(factor
, vE
[vert
]);
1359 void bfenumerator::handle(const signed_cone
& sc
, barvinok_options
*options
)
1361 assert(sc
.det
== 1);
1362 assert(sc
.rays
.NumRows() == enumerator_base::dim
);
1364 bfe_term
* t
= new bfe_term(enumerator_base::dim
);
1365 vector
< bfc_term_base
* > v
;
1368 t
->factors
.resize(1);
1370 t
->terms
.SetDims(1, enumerator_base::dim
);
1371 lattice_point(V
, sc
.rays
, t
->terms
[0], E_vertex
, options
);
1373 // the elements of factors are always lexpositive
1375 int s
= setup_factors(sc
.rays
, factors
, t
, sc
.sign
);
1377 t
->factors
[0] = new evalue
;
1378 value_init(t
->factors
[0]->d
);
1379 evalue_set_si(t
->factors
[0], s
, 1);
1380 reduce(factors
, v
, options
);
1382 for (int i
= 0; i
< enumerator_base::dim
; ++i
)
1384 free_evalue_refs(E_vertex
[i
]);
1389 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1390 barvinok_options
*options
);
1393 static evalue
* barvinok_enumerate_cst(Polyhedron
*P
, Polyhedron
* C
,
1394 struct barvinok_options
*options
)
1400 return evalue_zero();
1403 ALLOC(evalue
, eres
);
1404 value_init(eres
->d
);
1405 value_set_si(eres
->d
, 0);
1406 eres
->x
.p
= new_enode(partition
, 2, C
->Dimension
);
1407 EVALUE_SET_DOMAIN(eres
->x
.p
->arr
[0],
1408 DomainConstraintSimplify(C
, options
->MaxRays
));
1409 value_set_si(eres
->x
.p
->arr
[1].d
, 1);
1410 value_init(eres
->x
.p
->arr
[1].x
.n
);
1412 value_set_si(eres
->x
.p
->arr
[1].x
.n
, 0);
1414 barvinok_count_with_options(P
, &eres
->x
.p
->arr
[1].x
.n
, options
);
1419 static evalue
* enumerate(Polyhedron
*P
, Polyhedron
* C
,
1420 struct barvinok_options
*options
)
1423 Polyhedron
*Porig
= P
;
1424 Polyhedron
*Corig
= C
;
1425 Polyhedron
*CEq
= NULL
, *rVD
;
1427 unsigned nparam
= C
->Dimension
;
1432 value_init(factor
.d
);
1433 evalue_set_si(&factor
, 1, 1);
1436 POL_ENSURE_FACETS(P
);
1437 POL_ENSURE_VERTICES(P
);
1438 POL_ENSURE_FACETS(C
);
1439 POL_ENSURE_VERTICES(C
);
1441 if (C
->Dimension
== 0 || emptyQ(P
) || emptyQ(C
)) {
1444 CEq
= Polyhedron_Copy(CEq
);
1445 eres
= barvinok_enumerate_cst(P
, CEq
? CEq
: Polyhedron_Copy(C
), options
);
1448 evalue_backsubstitute(eres
, CP
, options
->MaxRays
);
1452 emul(&factor
, eres
);
1453 if (options
->approximation_method
== BV_APPROX_DROP
) {
1454 if (options
->polynomial_approximation
== BV_APPROX_SIGN_UPPER
)
1455 evalue_frac2polynomial(eres
, 1, options
->MaxRays
);
1456 if (options
->polynomial_approximation
== BV_APPROX_SIGN_LOWER
)
1457 evalue_frac2polynomial(eres
, -1, options
->MaxRays
);
1458 if (options
->polynomial_approximation
== BV_APPROX_SIGN_APPROX
)
1459 evalue_frac2polynomial(eres
, 0, options
->MaxRays
);
1461 reduce_evalue(eres
);
1462 free_evalue_refs(&factor
);
1470 if (Polyhedron_is_unbounded(P
, nparam
, options
->MaxRays
))
1475 P
= remove_equalities_p(Polyhedron_Copy(P
), P
->Dimension
-nparam
, &f
,
1477 mask(f
, &factor
, options
);
1480 if (P
->Dimension
== nparam
) {
1482 P
= Universe_Polyhedron(0);
1485 if (P
->NbEq
!= 0 || C
->NbEq
!= 0) {
1488 remove_all_equalities(&P
, &C
, &CP
, NULL
, nparam
, options
->MaxRays
);
1489 if (C
!= D
&& D
!= Corig
)
1491 if (P
!= Q
&& Q
!= Porig
)
1493 eres
= enumerate(P
, C
, options
);
1497 Polyhedron
*T
= Polyhedron_Factor(P
, nparam
, NULL
, options
->MaxRays
);
1498 if (T
|| (P
->Dimension
== nparam
+1)) {
1501 for (Q
= T
? T
: P
; Q
; Q
= Q
->next
) {
1502 Polyhedron
*next
= Q
->next
;
1506 if (Q
->Dimension
!= C
->Dimension
)
1507 QC
= Polyhedron_Project(Q
, nparam
);
1510 C
= DomainIntersection(C
, QC
, options
->MaxRays
);
1512 Polyhedron_Free(C2
);
1514 Polyhedron_Free(QC
);
1523 if (T
->Dimension
== C
->Dimension
) {
1532 eres
= barvinok_enumerate_ev_f(P
, C
, options
);
1539 for (Q
= P
->next
; Q
; Q
= Q
->next
) {
1540 Polyhedron
*next
= Q
->next
;
1543 f
= barvinok_enumerate_ev_f(Q
, C
, options
);
1554 evalue
* barvinok_enumerate_with_options(Polyhedron
*P
, Polyhedron
* C
,
1555 struct barvinok_options
*options
)
1557 Polyhedron
*next
, *Cnext
, *C1
;
1558 Polyhedron
*Corig
= C
;
1563 "barvinok_enumerate: input is a union; only first polyhedron is enumerated\n");
1567 "barvinok_enumerate: context is a union; only first polyhedron is considered\n");
1571 C1
= Polyhedron_Project(P
, C
->Dimension
);
1572 C
= DomainIntersection(C
, C1
, options
->MaxRays
);
1573 Polyhedron_Free(C1
);
1577 if (options
->approximation_method
== BV_APPROX_BERNOULLI
)
1578 eres
= Bernoulli_sum(P
, C
, options
);
1580 eres
= enumerate(P
, C
, options
);
1584 Corig
->next
= Cnext
;
1589 evalue
* barvinok_enumerate_ev(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1592 barvinok_options
*options
= barvinok_options_new_with_defaults();
1593 options
->MaxRays
= MaxRays
;
1594 E
= barvinok_enumerate_with_options(P
, C
, options
);
1595 barvinok_options_free(options
);
1599 evalue
*Param_Polyhedron_Enumerate(Param_Polyhedron
*PP
, Polyhedron
*P
,
1601 struct barvinok_options
*options
)
1605 unsigned nparam
= C
->Dimension
;
1606 unsigned dim
= P
->Dimension
- nparam
;
1609 for (nd
= 0, D
=PP
->D
; D
; ++nd
, D
=D
->next
);
1610 evalue_section
*s
= new evalue_section
[nd
];
1612 enumerator_base
*et
= NULL
;
1617 et
= enumerator_base::create(P
, dim
, PP
->nbV
, options
);
1619 Polyhedron
*TC
= true_context(P
, C
, options
->MaxRays
);
1620 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
1623 s
[i
].E
= evalue_zero();
1626 FORALL_PVertex_in_ParamPolyhedron(V
,D
,PP
) // _i is internal counter
1629 et
->decompose_at(V
, _i
, options
);
1630 } catch (OrthogonalException
&e
) {
1631 FORALL_REDUCED_DOMAIN_RESET
;
1632 for (; i
>= 0; --i
) {
1633 evalue_free(s
[i
].E
);
1634 Domain_Free(s
[i
].D
);
1638 eadd(et
->vE
[_i
] , s
[i
].E
);
1639 END_FORALL_PVertex_in_ParamPolyhedron
;
1640 evalue_range_reduction_in_domain(s
[i
].E
, rVD
);
1641 END_FORALL_REDUCED_DOMAIN
1642 Polyhedron_Free(TC
);
1645 eres
= evalue_from_section_array(s
, nd
);
1651 static evalue
* barvinok_enumerate_ev_f(Polyhedron
*P
, Polyhedron
* C
,
1652 barvinok_options
*options
)
1654 unsigned nparam
= C
->Dimension
;
1655 bool do_scale
= options
->approximation_method
== BV_APPROX_SCALE
;
1657 if (options
->approximation_method
== BV_APPROX_VOLUME
)
1658 return Param_Polyhedron_Volume(P
, C
, options
);
1660 if (P
->Dimension
- nparam
== 1 && !do_scale
)
1661 return ParamLine_Length(P
, C
, options
);
1663 Param_Polyhedron
*PP
= NULL
;
1667 eres
= scale_bound(P
, C
, options
);
1672 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
1675 eres
= scale(PP
, P
, C
, options
);
1677 eres
= Param_Polyhedron_Enumerate(PP
, P
, C
, options
);
1680 Param_Polyhedron_Free(PP
);
1685 Enumeration
* barvinok_enumerate(Polyhedron
*P
, Polyhedron
* C
, unsigned MaxRays
)
1687 evalue
*EP
= barvinok_enumerate_ev(P
, C
, MaxRays
);
1689 return partition2enumeration(EP
);
1692 evalue
* barvinok_enumerate_union(Polyhedron
*D
, Polyhedron
* C
, unsigned MaxRays
)
1695 gen_fun
*gf
= barvinok_enumerate_union_series(D
, C
, MaxRays
);
1701 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
1703 if (options
->summation
== BV_SUM_EULER
)
1704 return euler_summate(e
, nvar
, options
);
1705 else if (options
->summation
== BV_SUM_BERNOULLI
)
1706 return Bernoulli_sum_evalue(e
, nvar
, options
);
1708 return evalue_sum(e
, nvar
, options
->MaxRays
);