1 #include <barvinok/options.h>
2 #include <barvinok/util.h>
7 #include "section_array.h"
9 extern evalue
*evalue_outer_floor(evalue
*e
);
10 extern int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
);
11 extern void evalue_drop_floor(evalue
*e
, const evalue
*floor
);
13 #define ALLOC(type) (type*)malloc(sizeof(type))
14 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
16 /* Apply the variable transformation specified by T and CP on
17 * the polynomial e. T expresses the old variables in terms
18 * of the new variables (and optionally also the new parameters),
19 * while CP expresses the old parameters in terms of the new
22 static void transform_polynomial(evalue
*E
, Matrix
*T
, Matrix
*CP
,
23 unsigned nvar
, unsigned nparam
,
24 unsigned new_nvar
, unsigned new_nparam
)
29 subs
= ALLOCN(evalue
*, nvar
+nparam
);
31 for (j
= 0; j
< nvar
; ++j
) {
33 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[T
->NbRows
-1][T
->NbColumns
-1],
36 subs
[j
] = evalue_var(j
);
38 for (j
= 0; j
< nparam
; ++j
) {
40 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[nparam
][new_nparam
],
43 subs
[nvar
+j
] = evalue_var(j
);
44 evalue_shift_variables(subs
[nvar
+j
], 0, new_nvar
);
47 evalue_substitute(E
, subs
);
50 for (j
= 0; j
< nvar
+nparam
; ++j
)
55 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
57 struct evalue_section_array
*sections
,
58 struct barvinok_options
*options
)
60 unsigned dim
= P
->Dimension
;
61 unsigned new_dim
, new_nparam
;
62 Matrix
*T
= NULL
, *CP
= NULL
;
70 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
77 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
78 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
80 /* We can avoid these substitutions if E is a constant */
82 transform_polynomial(E
, T
, CP
, nvar
, dim
-nvar
,
83 new_dim
-new_nparam
, new_nparam
);
85 if (new_dim
-new_nparam
> 0) {
86 sum
= barvinok_sum_over_polytope(P
, E
, new_dim
-new_nparam
,
93 sum
->x
.p
= new_enode(partition
, 2, new_dim
);
94 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
95 value_clear(sum
->x
.p
->arr
[1].d
);
96 sum
->x
.p
->arr
[1] = *E
;
101 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
111 static evalue
*sum_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
112 struct barvinok_options
*options
)
114 if (options
->summation
== BV_SUM_EULER
)
115 return euler_summate(P
, E
, nvar
, options
);
116 else if (options
->summation
== BV_SUM_LAURENT
)
117 return laurent_summate(P
, E
, nvar
, options
);
121 /* Count the number of non-zero terms in e when viewed as a polynomial
122 * in only the first nvar variables. "count" is the number counted
125 static int evalue_count_terms(const evalue
*e
, unsigned nvar
, int count
)
129 if (EVALUE_IS_ZERO(*e
))
132 if (value_zero_p(e
->d
))
133 assert(e
->x
.p
->type
== polynomial
);
134 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1)
137 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
138 count
= evalue_count_terms(&e
->x
.p
->arr
[i
], nvar
, count
);
143 /* Create placeholder structure for unzipping.
144 * A "polynomial" is created with size terms in variable pos,
145 * with each term having itself as coefficient.
147 static evalue
*create_placeholder(int size
, int pos
)
150 evalue
*E
= ALLOC(evalue
);
152 E
->x
.p
= new_enode(polynomial
, size
, pos
+1);
153 for (i
= 0; i
< size
; ++i
) {
154 E
->x
.p
->arr
[i
].x
.p
= new_enode(polynomial
, i
+1, pos
+1);
155 for (j
= 0; j
< i
; ++j
)
156 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[j
], 0, 1);
157 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[i
], 1, 1);
162 /* Interchange each non-zero term in e (when viewed as a polynomial
163 * in only the first nvar variables) with a placeholder in ph (created
164 * by create_placeholder), resulting in two polynomials in the
165 * placeholder variable such that for each non-zero term in e
166 * there is a power of the placeholder variable such that the factors
167 * in the first nvar variables form the coefficient of that power in
168 * the first polynomial (e) and the factors in the remaining variables
169 * form the coefficient of that power in the second polynomial (ph).
171 static int evalue_unzip_terms(evalue
*e
, evalue
*ph
, unsigned nvar
, int count
)
175 if (EVALUE_IS_ZERO(*e
))
178 if (value_zero_p(e
->d
))
179 assert(e
->x
.p
->type
== polynomial
);
180 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1) {
182 *e
= ph
->x
.p
->arr
[count
];
183 ph
->x
.p
->arr
[count
] = t
;
187 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
188 count
= evalue_unzip_terms(&e
->x
.p
->arr
[i
], ph
, nvar
, count
);
193 /* Remove n variables at pos (0-based) from the polyhedron P.
194 * Each of these variables is assumed to be completely free,
195 * i.e., there is a line in the polyhedron corresponding to
196 * each of these variables.
198 static Polyhedron
*Polyhedron_Remove_Columns(Polyhedron
*P
, unsigned pos
,
202 unsigned NbConstraints
= 0;
211 assert(pos
<= P
->Dimension
);
213 if (POL_HAS(P
, POL_INEQUALITIES
))
214 NbConstraints
= P
->NbConstraints
;
215 if (POL_HAS(P
, POL_POINTS
))
216 NbRays
= P
->NbRays
- n
;
218 Q
= Polyhedron_Alloc(P
->Dimension
- n
, NbConstraints
, NbRays
);
219 if (POL_HAS(P
, POL_INEQUALITIES
)) {
221 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
222 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
223 Vector_Copy(P
->Constraint
[i
]+1+pos
+n
, Q
->Constraint
[i
]+1+pos
,
227 if (POL_HAS(P
, POL_POINTS
)) {
228 Q
->NbBid
= P
->NbBid
- n
;
229 for (i
= 0; i
< n
; ++i
)
230 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
231 for (i
= 0, j
= 0; i
< P
->NbRays
; ++i
) {
232 int line
= First_Non_Zero(P
->Ray
[i
], 1+P
->Dimension
+1);
234 if (line
-1 >= pos
&& line
-1 < pos
+n
) {
239 assert(i
-j
< Q
->NbRays
);
240 Vector_Copy(P
->Ray
[i
], Q
->Ray
[i
-j
], 1+pos
);
241 Vector_Copy(P
->Ray
[i
]+1+pos
+n
, Q
->Ray
[i
-j
]+1+pos
,
245 POL_SET(Q
, POL_VALID
);
246 if (POL_HAS(P
, POL_INEQUALITIES
))
247 POL_SET(Q
, POL_INEQUALITIES
);
248 if (POL_HAS(P
, POL_POINTS
))
249 POL_SET(Q
, POL_POINTS
);
250 if (POL_HAS(P
, POL_VERTICES
))
251 POL_SET(Q
, POL_VERTICES
);
255 /* Remove n variables at pos (0-based) from the union of polyhedra P.
256 * Each of these variables is assumed to be completely free,
257 * i.e., there is a line in the polyhedron corresponding to
258 * each of these variables.
260 static Polyhedron
*Domain_Remove_Columns(Polyhedron
*P
, unsigned pos
,
264 Polyhedron
**next
= &R
;
266 for (; P
; P
= P
->next
) {
267 *next
= Polyhedron_Remove_Columns(P
, pos
, n
);
268 next
= &(*next
)->next
;
273 /* Drop n parameters starting at first from partition evalue e */
274 static void drop_parameters(evalue
*e
, int first
, int n
)
278 if (EVALUE_IS_ZERO(*e
))
281 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
282 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
283 Polyhedron
*P
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
284 Polyhedron
*Q
= Domain_Remove_Columns(P
, first
, n
);
285 EVALUE_SET_DOMAIN(e
->x
.p
->arr
[2*i
], Q
);
287 evalue_shift_variables(&e
->x
.p
->arr
[2*i
+1], first
, -n
);
292 static void extract_term_into(const evalue
*src
, int var
, int exp
, evalue
*dst
)
296 if (value_notzero_p(src
->d
) ||
297 src
->x
.p
->type
!= polynomial
||
298 src
->x
.p
->pos
> var
+1) {
300 evalue_copy(dst
, src
);
302 evalue_set_si(dst
, 0, 1);
306 if (src
->x
.p
->pos
== var
+1) {
307 if (src
->x
.p
->size
> exp
)
308 evalue_copy(dst
, &src
->x
.p
->arr
[exp
]);
310 evalue_set_si(dst
, 0, 1);
314 dst
->x
.p
= new_enode(polynomial
, src
->x
.p
->size
, src
->x
.p
->pos
);
315 for (i
= 0; i
< src
->x
.p
->size
; ++i
)
316 extract_term_into(&src
->x
.p
->arr
[i
], var
, exp
,
320 /* Extract the coefficient of var^exp.
322 static evalue
*extract_term(const evalue
*e
, int var
, int exp
)
327 if (EVALUE_IS_ZERO(*e
))
328 return evalue_zero();
330 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
333 res
->x
.p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
334 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
335 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
336 Domain_Copy(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])));
337 extract_term_into(&e
->x
.p
->arr
[2*i
+1], var
, exp
,
338 &res
->x
.p
->arr
[2*i
+1]);
339 reduce_evalue(&res
->x
.p
->arr
[2*i
+1]);
344 /* Insert n free variables at pos (0-based) in the polyhedron P.
346 static Polyhedron
*Polyhedron_Insert_Columns(Polyhedron
*P
, unsigned pos
,
350 unsigned NbConstraints
= 0;
359 assert(pos
<= P
->Dimension
);
361 if (POL_HAS(P
, POL_INEQUALITIES
))
362 NbConstraints
= P
->NbConstraints
;
363 if (POL_HAS(P
, POL_POINTS
))
364 NbRays
= P
->NbRays
+ n
;
366 Q
= Polyhedron_Alloc(P
->Dimension
+n
, NbConstraints
, NbRays
);
367 if (POL_HAS(P
, POL_INEQUALITIES
)) {
369 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
370 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
371 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[i
]+1+pos
+n
,
375 if (POL_HAS(P
, POL_POINTS
)) {
376 Q
->NbBid
= P
->NbBid
+ n
;
377 for (i
= 0; i
< n
; ++i
)
378 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
379 for (i
= 0; i
< P
->NbRays
; ++i
) {
380 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[n
+i
], 1+pos
);
381 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[n
+i
]+1+pos
+n
,
385 POL_SET(Q
, POL_VALID
);
386 if (POL_HAS(P
, POL_INEQUALITIES
))
387 POL_SET(Q
, POL_INEQUALITIES
);
388 if (POL_HAS(P
, POL_POINTS
))
389 POL_SET(Q
, POL_POINTS
);
390 if (POL_HAS(P
, POL_VERTICES
))
391 POL_SET(Q
, POL_VERTICES
);
395 /* Perform summation of e over a list of 1 or more factors F, with context C.
396 * nvar is the total number of variables in the remaining factors.
397 * extra is the number of placeholder parameters introduced in e,
398 * but not (yet) in F or C.
400 * If there is only one factor left, F is intersected with the
401 * context C, the placeholder variables are added, and then
402 * e is summed over the resulting parametric polytope.
404 * If there is more than one factor left, we create to polynomials
405 * in a new placeholder variable (which is placed after the regular
406 * parameters, but before any previously introduced placeholder
407 * variables) that has the factors of the variables in the first
408 * factor of F and the factor of the remaining variables of
409 * each term as its coefficients.
410 * These two polynomials are then summed over their domains
411 * and afterwards the results are combined and the placeholder
412 * variable is removed again.
414 static evalue
*sum_factors(Polyhedron
*F
, Polyhedron
*C
, evalue
*e
,
415 unsigned nvar
, unsigned extra
,
416 struct barvinok_options
*options
)
419 unsigned nparam
= C
->Dimension
;
420 unsigned F_var
= F
->Dimension
- C
->Dimension
;
426 Polyhedron
*CA
= align_context(C
, nvar
+nparam
, options
->MaxRays
);
427 Polyhedron
*P
= DomainIntersection(F
, CA
, options
->MaxRays
);
428 Polyhedron
*Q
= Polyhedron_Insert_Columns(P
, nvar
+nparam
, extra
);
432 evalue
*sum
= sum_base(Q
, e
, nvar
, options
);
437 n
= evalue_count_terms(e
, F_var
, 0);
438 ph
= create_placeholder(n
, nvar
+nparam
);
439 evalue_shift_variables(e
, nvar
+nparam
, 1);
440 evalue_unzip_terms(e
, ph
, F_var
, 0);
441 evalue_shift_variables(e
, nvar
, -(nvar
-F_var
));
442 evalue_reorder_terms(ph
);
443 evalue_shift_variables(ph
, 0, -F_var
);
445 s2
= sum_factors(F
->next
, C
, ph
, nvar
-F_var
, extra
+1, options
);
448 s1
= sum_factors(F
, C
, e
, F_var
, extra
+1, options
);
451 /* remove placeholder "polynomial" */
455 drop_parameters(s2
, nparam
, 1);
460 for (i
= 0; i
< n
; ++i
) {
462 t1
= extract_term(s1
, nparam
, i
);
463 t2
= extract_term(s2
, nparam
, i
);
472 drop_parameters(s
, nparam
, 1);
476 /* Perform summation over a product of factors F, obtained using
477 * variable transformation T from the original problem specification.
479 * We first perform the corresponding transformation on the polynomial E,
480 * compute the common context over all factors and then perform
481 * the actual summation over the factors.
483 static evalue
*sum_product(Polyhedron
*F
, evalue
*E
, Matrix
*T
, unsigned nparam
,
484 struct barvinok_options
*options
)
488 unsigned nvar
= T
->NbRows
;
492 assert(nvar
== T
->NbColumns
);
493 T2
= Matrix_Alloc(nvar
+1, nvar
+1);
494 for (i
= 0; i
< nvar
; ++i
)
495 Vector_Copy(T
->p
[i
], T2
->p
[i
], nvar
);
496 value_set_si(T2
->p
[nvar
][nvar
], 1);
498 transform_polynomial(E
, T2
, NULL
, nvar
, nparam
, nvar
, nparam
);
500 C
= Factor_Context(F
, nparam
, options
->MaxRays
);
501 if (F
->Dimension
== nparam
) {
507 sum
= sum_factors(F
, C
, E
, nvar
, 0, options
);
515 /* Add two constraints corresponding to floor = floor(e/d),
518 * -e + d t + d-1 >= 0
520 * e is assumed to be an affine expression.
522 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
523 struct barvinok_options
*options
)
526 unsigned dim
= P
->Dimension
+1;
527 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
529 Value
*d
= &M
->p
[0][1+nvar
];
530 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
531 value_oppose(*d
, *d
);
532 value_set_si(M
->p
[0][0], 1);
533 value_set_si(M
->p
[1][0], 1);
534 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
535 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
536 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
538 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
539 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
540 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
543 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
548 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
559 /* Compute sum of a step-polynomial over a polytope by grouping
560 * terms containing the same floor-expressions and introducing
561 * new variables for each such expression.
562 * In particular, while there is any floor-expression left,
563 * the step-polynomial is split into a polynomial containing
564 * the expression, which is then converted to a new variable,
565 * and a polynomial not containing the expression.
567 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
568 struct barvinok_options
*options
)
575 while ((floor
= evalue_outer_floor(cur
))) {
578 evalue
*converted_floor
;
580 /* Ignore floors that do not depend on variables. */
581 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
584 converted
= evalue_dup(cur
);
585 converted_floor
= evalue_dup(floor
);
586 evalue_shift_variables(converted
, nvar
, 1);
587 evalue_shift_variables(converted_floor
, nvar
, 1);
588 evalue_replace_floor(converted
, converted_floor
, nvar
);
589 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
590 evalue_free(converted_floor
);
591 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
592 evalue_free(converted
);
594 sum
= evalue_add(t
, sum
);
597 cur
= evalue_dup(cur
);
598 evalue_drop_floor(cur
, floor
);
602 evalue_floor2frac(cur
);
606 if (EVALUE_IS_ZERO(*cur
))
610 unsigned nparam
= P
->Dimension
- nvar
;
611 Polyhedron
*F
= Polyhedron_Factor(P
, nparam
, &T
, options
->MaxRays
);
613 t
= sum_base(P
, cur
, nvar
, options
);
616 cur
= evalue_dup(cur
);
617 t
= sum_product(F
, cur
, T
, nparam
, options
);
624 return evalue_add(t
, sum
);
627 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
628 struct evalue_section_array
*sections
,
629 struct barvinok_options
*options
)
632 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
634 if (options
->summation
== BV_SUM_BERNOULLI
)
635 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
636 else if (options
->summation
== BV_SUM_BOX
)
637 return box_summate(P
, E
, nvar
, options
->MaxRays
);
639 evalue_frac2floor2(E
, 0);
641 return sum_step_polynomial(P
, E
, nvar
, options
);
644 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
647 struct evalue_section_array sections
;
651 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
652 return evalue_dup(e
);
654 assert(value_zero_p(e
->d
));
655 assert(e
->x
.p
->type
== partition
);
657 evalue_section_array_init(§ions
);
660 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
662 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
663 Polyhedron
*next
= D
->next
;
667 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
683 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
686 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
687 options
->MaxRays
= MaxRays
;
688 sum
= barvinok_summate(E
, nvar
, options
);
689 barvinok_options_free(options
);
693 evalue
*esum(evalue
*e
, int nvar
)
696 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
697 sum
= barvinok_summate(e
, nvar
, options
);
698 barvinok_options_free(options
);
702 /* Turn unweighted counting problem into "weighted" counting problem
703 * with weight equal to 1 and call barvinok_summate on this weighted problem.
705 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
706 struct barvinok_options
*options
)
712 if (emptyQ(P
) || emptyQ(C
))
713 return evalue_zero();
715 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
716 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
721 return evalue_zero();
725 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
726 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
727 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
728 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
729 free_evalue_refs(&e
);