1 #include <barvinok/options.h>
2 #include <barvinok/util.h>
7 #include "section_array.h"
8 #include "remove_equalities.h"
10 extern evalue
*evalue_outer_floor(evalue
*e
);
11 extern int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
);
12 extern void evalue_drop_floor(evalue
*e
, const evalue
*floor
);
14 #define ALLOC(type) (type*)malloc(sizeof(type))
15 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
17 /* Apply the variable transformation specified by T and CP on
18 * the polynomial e. T expresses the old variables in terms
19 * of the new variables (and optionally also the new parameters),
20 * while CP expresses the old parameters in terms of the new
23 static void transform_polynomial(evalue
*E
, Matrix
*T
, Matrix
*CP
,
24 unsigned nvar
, unsigned nparam
,
25 unsigned new_nvar
, unsigned new_nparam
)
30 subs
= ALLOCN(evalue
*, nvar
+nparam
);
32 for (j
= 0; j
< nvar
; ++j
) {
34 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[T
->NbRows
-1][T
->NbColumns
-1],
37 subs
[j
] = evalue_var(j
);
39 for (j
= 0; j
< nparam
; ++j
) {
41 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[nparam
][new_nparam
],
44 subs
[nvar
+j
] = evalue_var(j
);
45 evalue_shift_variables(subs
[nvar
+j
], 0, new_nvar
);
48 evalue_substitute(E
, subs
);
51 for (j
= 0; j
< nvar
+nparam
; ++j
)
56 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
58 struct evalue_section_array
*sections
,
59 struct barvinok_options
*options
)
61 unsigned dim
= P
->Dimension
;
62 unsigned new_dim
, new_nparam
;
63 Matrix
*T
= NULL
, *CP
= NULL
;
71 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
78 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
79 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
81 /* We can avoid these substitutions if E is a constant */
83 transform_polynomial(E
, T
, CP
, nvar
, dim
-nvar
,
84 new_dim
-new_nparam
, new_nparam
);
86 if (new_dim
-new_nparam
> 0) {
87 sum
= barvinok_sum_over_polytope(P
, E
, new_dim
-new_nparam
,
94 sum
->x
.p
= new_enode(partition
, 2, new_dim
);
95 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
96 value_clear(sum
->x
.p
->arr
[1].d
);
97 sum
->x
.p
->arr
[1] = *E
;
102 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
112 static evalue
*sum_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
113 struct barvinok_options
*options
)
115 if (options
->summation
== BV_SUM_EULER
)
116 return euler_summate(P
, E
, nvar
, options
);
117 else if (options
->summation
== BV_SUM_LAURENT
)
118 return laurent_summate(P
, E
, nvar
, options
);
119 else if (options
->summation
== BV_SUM_LAURENT_OLD
)
120 return laurent_summate_old(P
, E
, nvar
, options
);
124 /* Count the number of non-zero terms in e when viewed as a polynomial
125 * in only the first nvar variables. "count" is the number counted
128 static int evalue_count_terms(const evalue
*e
, unsigned nvar
, int count
)
132 if (EVALUE_IS_ZERO(*e
))
135 if (value_zero_p(e
->d
))
136 assert(e
->x
.p
->type
== polynomial
);
137 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1)
140 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
141 count
= evalue_count_terms(&e
->x
.p
->arr
[i
], nvar
, count
);
146 /* Create placeholder structure for unzipping.
147 * A "polynomial" is created with size terms in variable pos,
148 * with each term having itself as coefficient.
150 static evalue
*create_placeholder(int size
, int pos
)
153 evalue
*E
= ALLOC(evalue
);
155 E
->x
.p
= new_enode(polynomial
, size
, pos
+1);
156 for (i
= 0; i
< size
; ++i
) {
157 E
->x
.p
->arr
[i
].x
.p
= new_enode(polynomial
, i
+1, pos
+1);
158 for (j
= 0; j
< i
; ++j
)
159 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[j
], 0, 1);
160 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[i
], 1, 1);
165 /* Interchange each non-zero term in e (when viewed as a polynomial
166 * in only the first nvar variables) with a placeholder in ph (created
167 * by create_placeholder), resulting in two polynomials in the
168 * placeholder variable such that for each non-zero term in e
169 * there is a power of the placeholder variable such that the factors
170 * in the first nvar variables form the coefficient of that power in
171 * the first polynomial (e) and the factors in the remaining variables
172 * form the coefficient of that power in the second polynomial (ph).
174 static int evalue_unzip_terms(evalue
*e
, evalue
*ph
, unsigned nvar
, int count
)
178 if (EVALUE_IS_ZERO(*e
))
181 if (value_zero_p(e
->d
))
182 assert(e
->x
.p
->type
== polynomial
);
183 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1) {
185 *e
= ph
->x
.p
->arr
[count
];
186 ph
->x
.p
->arr
[count
] = t
;
190 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
191 count
= evalue_unzip_terms(&e
->x
.p
->arr
[i
], ph
, nvar
, count
);
196 /* Remove n variables at pos (0-based) from the polyhedron P.
197 * Each of these variables is assumed to be completely free,
198 * i.e., there is a line in the polyhedron corresponding to
199 * each of these variables.
201 static Polyhedron
*Polyhedron_Remove_Columns(Polyhedron
*P
, unsigned pos
,
205 unsigned NbConstraints
= 0;
212 assert(pos
<= P
->Dimension
);
214 if (POL_HAS(P
, POL_INEQUALITIES
))
215 NbConstraints
= P
->NbConstraints
;
216 if (POL_HAS(P
, POL_POINTS
))
217 NbRays
= P
->NbRays
- n
;
219 Q
= Polyhedron_Alloc(P
->Dimension
- n
, NbConstraints
, NbRays
);
220 if (POL_HAS(P
, POL_INEQUALITIES
)) {
222 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
223 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
224 Vector_Copy(P
->Constraint
[i
]+1+pos
+n
, Q
->Constraint
[i
]+1+pos
,
228 if (POL_HAS(P
, POL_POINTS
)) {
229 Q
->NbBid
= P
->NbBid
- n
;
230 for (i
= 0; i
< n
; ++i
)
231 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
232 for (i
= 0, j
= 0; i
< P
->NbRays
; ++i
) {
233 int line
= First_Non_Zero(P
->Ray
[i
], 1+P
->Dimension
+1);
235 if (line
-1 >= pos
&& line
-1 < pos
+n
) {
240 assert(i
-j
< Q
->NbRays
);
241 Vector_Copy(P
->Ray
[i
], Q
->Ray
[i
-j
], 1+pos
);
242 Vector_Copy(P
->Ray
[i
]+1+pos
+n
, Q
->Ray
[i
-j
]+1+pos
,
246 POL_SET(Q
, POL_VALID
);
247 if (POL_HAS(P
, POL_INEQUALITIES
))
248 POL_SET(Q
, POL_INEQUALITIES
);
249 if (POL_HAS(P
, POL_POINTS
))
250 POL_SET(Q
, POL_POINTS
);
251 if (POL_HAS(P
, POL_VERTICES
))
252 POL_SET(Q
, POL_VERTICES
);
256 /* Remove n variables at pos (0-based) from the union of polyhedra P.
257 * Each of these variables is assumed to be completely free,
258 * i.e., there is a line in the polyhedron corresponding to
259 * each of these variables.
261 static Polyhedron
*Domain_Remove_Columns(Polyhedron
*P
, unsigned pos
,
265 Polyhedron
**next
= &R
;
267 for (; P
; P
= P
->next
) {
268 *next
= Polyhedron_Remove_Columns(P
, pos
, n
);
269 next
= &(*next
)->next
;
274 /* Drop n parameters starting at first from partition evalue e */
275 static void drop_parameters(evalue
*e
, int first
, int n
)
279 if (EVALUE_IS_ZERO(*e
))
282 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
283 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
284 Polyhedron
*P
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
285 Polyhedron
*Q
= Domain_Remove_Columns(P
, first
, n
);
286 EVALUE_SET_DOMAIN(e
->x
.p
->arr
[2*i
], Q
);
288 evalue_shift_variables(&e
->x
.p
->arr
[2*i
+1], first
, -n
);
293 static void extract_term_into(const evalue
*src
, int var
, int exp
, evalue
*dst
)
297 if (value_notzero_p(src
->d
) ||
298 src
->x
.p
->type
!= polynomial
||
299 src
->x
.p
->pos
> var
+1) {
301 evalue_copy(dst
, src
);
303 evalue_set_si(dst
, 0, 1);
307 if (src
->x
.p
->pos
== var
+1) {
308 if (src
->x
.p
->size
> exp
)
309 evalue_copy(dst
, &src
->x
.p
->arr
[exp
]);
311 evalue_set_si(dst
, 0, 1);
315 dst
->x
.p
= new_enode(polynomial
, src
->x
.p
->size
, src
->x
.p
->pos
);
316 for (i
= 0; i
< src
->x
.p
->size
; ++i
)
317 extract_term_into(&src
->x
.p
->arr
[i
], var
, exp
,
321 /* Extract the coefficient of var^exp.
323 static evalue
*extract_term(const evalue
*e
, int var
, int exp
)
328 if (EVALUE_IS_ZERO(*e
))
329 return evalue_zero();
331 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
334 res
->x
.p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
335 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
336 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
337 Domain_Copy(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])));
338 extract_term_into(&e
->x
.p
->arr
[2*i
+1], var
, exp
,
339 &res
->x
.p
->arr
[2*i
+1]);
340 reduce_evalue(&res
->x
.p
->arr
[2*i
+1]);
345 /* Insert n free variables at pos (0-based) in the polyhedron P.
347 static Polyhedron
*Polyhedron_Insert_Columns(Polyhedron
*P
, unsigned pos
,
351 unsigned NbConstraints
= 0;
360 assert(pos
<= P
->Dimension
);
362 if (POL_HAS(P
, POL_INEQUALITIES
))
363 NbConstraints
= P
->NbConstraints
;
364 if (POL_HAS(P
, POL_POINTS
))
365 NbRays
= P
->NbRays
+ n
;
367 Q
= Polyhedron_Alloc(P
->Dimension
+n
, NbConstraints
, NbRays
);
368 if (POL_HAS(P
, POL_INEQUALITIES
)) {
370 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
371 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
372 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[i
]+1+pos
+n
,
376 if (POL_HAS(P
, POL_POINTS
)) {
377 Q
->NbBid
= P
->NbBid
+ n
;
378 for (i
= 0; i
< n
; ++i
)
379 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
380 for (i
= 0; i
< P
->NbRays
; ++i
) {
381 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[n
+i
], 1+pos
);
382 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[n
+i
]+1+pos
+n
,
386 POL_SET(Q
, POL_VALID
);
387 if (POL_HAS(P
, POL_INEQUALITIES
))
388 POL_SET(Q
, POL_INEQUALITIES
);
389 if (POL_HAS(P
, POL_POINTS
))
390 POL_SET(Q
, POL_POINTS
);
391 if (POL_HAS(P
, POL_VERTICES
))
392 POL_SET(Q
, POL_VERTICES
);
396 /* Perform summation of e over a list of 1 or more factors F, with context C.
397 * nvar is the total number of variables in the remaining factors.
398 * extra is the number of placeholder parameters introduced in e,
399 * but not (yet) in F or C.
401 * If there is only one factor left, F is intersected with the
402 * context C, the placeholder variables are added, and then
403 * e is summed over the resulting parametric polytope.
405 * If there is more than one factor left, we create two polynomials
406 * in a new placeholder variable (which is placed after the regular
407 * parameters, but before any previously introduced placeholder
408 * variables) that has the factors of the variables in the first
409 * factor of F and the factor of the remaining variables of
410 * each term as its coefficients.
411 * These two polynomials are then summed over their domains
412 * and afterwards the results are combined and the placeholder
413 * variable is removed again.
415 static evalue
*sum_factors(Polyhedron
*F
, Polyhedron
*C
, evalue
*e
,
416 unsigned nvar
, unsigned extra
,
417 struct barvinok_options
*options
)
420 unsigned nparam
= C
->Dimension
;
421 unsigned F_var
= F
->Dimension
- C
->Dimension
;
427 Polyhedron
*CA
= align_context(C
, nvar
+nparam
, options
->MaxRays
);
428 Polyhedron
*P
= DomainIntersection(F
, CA
, options
->MaxRays
);
429 Polyhedron
*Q
= Polyhedron_Insert_Columns(P
, nvar
+nparam
, extra
);
433 evalue
*sum
= sum_base(Q
, e
, nvar
, options
);
438 n
= evalue_count_terms(e
, F_var
, 0);
439 ph
= create_placeholder(n
, nvar
+nparam
);
440 evalue_shift_variables(e
, nvar
+nparam
, 1);
441 evalue_unzip_terms(e
, ph
, F_var
, 0);
442 evalue_shift_variables(e
, nvar
, -(nvar
-F_var
));
443 evalue_reorder_terms(ph
);
444 evalue_shift_variables(ph
, 0, -F_var
);
446 s2
= sum_factors(F
->next
, C
, ph
, nvar
-F_var
, extra
+1, options
);
449 s1
= sum_factors(F
, C
, e
, F_var
, extra
+1, options
);
452 /* remove placeholder "polynomial" */
456 drop_parameters(s2
, nparam
, 1);
461 for (i
= 0; i
< n
; ++i
) {
463 t1
= extract_term(s1
, nparam
, i
);
464 t2
= extract_term(s2
, nparam
, i
);
473 drop_parameters(s
, nparam
, 1);
477 /* Perform summation over a product of factors F, obtained using
478 * variable transformation T from the original problem specification.
480 * We first perform the corresponding transformation on the polynomial E,
481 * compute the common context over all factors and then perform
482 * the actual summation over the factors.
484 static evalue
*sum_product(Polyhedron
*F
, evalue
*E
, Matrix
*T
, unsigned nparam
,
485 struct barvinok_options
*options
)
489 unsigned nvar
= T
->NbRows
;
493 assert(nvar
== T
->NbColumns
);
494 T2
= Matrix_Alloc(nvar
+1, nvar
+1);
495 for (i
= 0; i
< nvar
; ++i
)
496 Vector_Copy(T
->p
[i
], T2
->p
[i
], nvar
);
497 value_set_si(T2
->p
[nvar
][nvar
], 1);
499 transform_polynomial(E
, T2
, NULL
, nvar
, nparam
, nvar
, nparam
);
501 C
= Factor_Context(F
, nparam
, options
->MaxRays
);
502 if (F
->Dimension
== nparam
) {
508 sum
= sum_factors(F
, C
, E
, nvar
, 0, options
);
516 /* Add two constraints corresponding to floor = floor(e/d),
519 * -e + d t + d-1 >= 0
521 * e is assumed to be an affine expression.
523 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
524 struct barvinok_options
*options
)
527 unsigned dim
= P
->Dimension
+1;
528 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
530 Value
*d
= &M
->p
[0][1+nvar
];
531 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
532 value_oppose(*d
, *d
);
533 value_set_si(M
->p
[0][0], 1);
534 value_set_si(M
->p
[1][0], 1);
535 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
536 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
537 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
539 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
540 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
541 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
544 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
549 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
560 /* Compute sum of a step-polynomial over a polytope by grouping
561 * terms containing the same floor-expressions and introducing
562 * new variables for each such expression.
563 * In particular, while there is any floor-expression left,
564 * the step-polynomial is split into a polynomial containing
565 * the expression, which is then converted to a new variable,
566 * and a polynomial not containing the expression.
568 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
569 struct barvinok_options
*options
)
576 while ((floor
= evalue_outer_floor(cur
))) {
579 evalue
*converted_floor
;
581 /* Ignore floors that do not depend on variables. */
582 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
585 converted
= evalue_dup(cur
);
586 converted_floor
= evalue_dup(floor
);
587 evalue_shift_variables(converted
, nvar
, 1);
588 evalue_shift_variables(converted_floor
, nvar
, 1);
589 evalue_replace_floor(converted
, converted_floor
, nvar
);
590 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
591 evalue_free(converted_floor
);
592 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
593 evalue_free(converted
);
595 sum
= evalue_add(t
, sum
);
598 cur
= evalue_dup(cur
);
599 evalue_drop_floor(cur
, floor
);
603 evalue_floor2frac(cur
);
607 if (EVALUE_IS_ZERO(*cur
))
611 unsigned nparam
= P
->Dimension
- nvar
;
612 Polyhedron
*F
= Polyhedron_Factor(P
, nparam
, &T
, options
->MaxRays
);
614 t
= sum_base(P
, cur
, nvar
, options
);
617 cur
= evalue_dup(cur
);
618 t
= sum_product(F
, cur
, T
, nparam
, options
);
625 return evalue_add(t
, sum
);
628 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
629 struct evalue_section_array
*sections
,
630 struct barvinok_options
*options
)
633 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
635 if (options
->summation
== BV_SUM_BERNOULLI
)
636 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
637 else if (options
->summation
== BV_SUM_BOX
)
638 return box_summate(P
, E
, nvar
, options
->MaxRays
);
640 evalue_frac2floor2(E
, 0);
642 return sum_step_polynomial(P
, E
, nvar
, options
);
645 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
648 struct evalue_section_array sections
;
652 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
653 return evalue_dup(e
);
655 assert(value_zero_p(e
->d
));
656 assert(e
->x
.p
->type
== partition
);
658 evalue_section_array_init(§ions
);
661 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
663 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
664 Polyhedron
*next
= D
->next
;
668 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
684 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
687 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
688 options
->MaxRays
= MaxRays
;
689 sum
= barvinok_summate(E
, nvar
, options
);
690 barvinok_options_free(options
);
694 evalue
*esum(evalue
*e
, int nvar
)
697 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
698 sum
= barvinok_summate(e
, nvar
, options
);
699 barvinok_options_free(options
);
703 /* Turn unweighted counting problem into "weighted" counting problem
704 * with weight equal to 1 and call barvinok_summate on this weighted problem.
706 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
707 struct barvinok_options
*options
)
713 if (emptyQ(P
) || emptyQ(C
))
714 return evalue_zero();
716 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
717 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
722 return evalue_zero();
726 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
727 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
728 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
729 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
730 free_evalue_refs(&e
);