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
);
122 /* Count the number of non-zero terms in e when viewed as a polynomial
123 * in only the first nvar variables. "count" is the number counted
126 static int evalue_count_terms(const evalue
*e
, unsigned nvar
, int count
)
130 if (EVALUE_IS_ZERO(*e
))
133 if (value_zero_p(e
->d
))
134 assert(e
->x
.p
->type
== polynomial
);
135 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1)
138 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
139 count
= evalue_count_terms(&e
->x
.p
->arr
[i
], nvar
, count
);
144 /* Create placeholder structure for unzipping.
145 * A "polynomial" is created with size terms in variable pos,
146 * with each term having itself as coefficient.
148 static evalue
*create_placeholder(int size
, int pos
)
151 evalue
*E
= ALLOC(evalue
);
153 E
->x
.p
= new_enode(polynomial
, size
, pos
+1);
154 for (i
= 0; i
< size
; ++i
) {
155 E
->x
.p
->arr
[i
].x
.p
= new_enode(polynomial
, i
+1, pos
+1);
156 for (j
= 0; j
< i
; ++j
)
157 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[j
], 0, 1);
158 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[i
], 1, 1);
163 /* Interchange each non-zero term in e (when viewed as a polynomial
164 * in only the first nvar variables) with a placeholder in ph (created
165 * by create_placeholder), resulting in two polynomials in the
166 * placeholder variable such that for each non-zero term in e
167 * there is a power of the placeholder variable such that the factors
168 * in the first nvar variables form the coefficient of that power in
169 * the first polynomial (e) and the factors in the remaining variables
170 * form the coefficient of that power in the second polynomial (ph).
172 static int evalue_unzip_terms(evalue
*e
, evalue
*ph
, unsigned nvar
, int count
)
176 if (EVALUE_IS_ZERO(*e
))
179 if (value_zero_p(e
->d
))
180 assert(e
->x
.p
->type
== polynomial
);
181 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1) {
183 *e
= ph
->x
.p
->arr
[count
];
184 ph
->x
.p
->arr
[count
] = t
;
188 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
189 count
= evalue_unzip_terms(&e
->x
.p
->arr
[i
], ph
, nvar
, count
);
194 /* Remove n variables at pos (0-based) from the polyhedron P.
195 * Each of these variables is assumed to be completely free,
196 * i.e., there is a line in the polyhedron corresponding to
197 * each of these variables.
199 static Polyhedron
*Polyhedron_Remove_Columns(Polyhedron
*P
, unsigned pos
,
203 unsigned NbConstraints
= 0;
210 assert(pos
<= P
->Dimension
);
212 if (POL_HAS(P
, POL_INEQUALITIES
))
213 NbConstraints
= P
->NbConstraints
;
214 if (POL_HAS(P
, POL_POINTS
))
215 NbRays
= P
->NbRays
- n
;
217 Q
= Polyhedron_Alloc(P
->Dimension
- n
, NbConstraints
, NbRays
);
218 if (POL_HAS(P
, POL_INEQUALITIES
)) {
220 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
221 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
222 Vector_Copy(P
->Constraint
[i
]+1+pos
+n
, Q
->Constraint
[i
]+1+pos
,
226 if (POL_HAS(P
, POL_POINTS
)) {
227 Q
->NbBid
= P
->NbBid
- n
;
228 for (i
= 0; i
< n
; ++i
)
229 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
230 for (i
= 0, j
= 0; i
< P
->NbRays
; ++i
) {
231 int line
= First_Non_Zero(P
->Ray
[i
], 1+P
->Dimension
+1);
233 if (line
-1 >= pos
&& line
-1 < pos
+n
) {
238 assert(i
-j
< Q
->NbRays
);
239 Vector_Copy(P
->Ray
[i
], Q
->Ray
[i
-j
], 1+pos
);
240 Vector_Copy(P
->Ray
[i
]+1+pos
+n
, Q
->Ray
[i
-j
]+1+pos
,
244 POL_SET(Q
, POL_VALID
);
245 if (POL_HAS(P
, POL_INEQUALITIES
))
246 POL_SET(Q
, POL_INEQUALITIES
);
247 if (POL_HAS(P
, POL_POINTS
))
248 POL_SET(Q
, POL_POINTS
);
249 if (POL_HAS(P
, POL_VERTICES
))
250 POL_SET(Q
, POL_VERTICES
);
254 /* Remove n variables at pos (0-based) from the union of polyhedra P.
255 * Each of these variables is assumed to be completely free,
256 * i.e., there is a line in the polyhedron corresponding to
257 * each of these variables.
259 static Polyhedron
*Domain_Remove_Columns(Polyhedron
*P
, unsigned pos
,
263 Polyhedron
**next
= &R
;
265 for (; P
; P
= P
->next
) {
266 *next
= Polyhedron_Remove_Columns(P
, pos
, n
);
267 next
= &(*next
)->next
;
272 /* Drop n parameters starting at first from partition evalue e */
273 static void drop_parameters(evalue
*e
, int first
, int n
)
277 if (EVALUE_IS_ZERO(*e
))
280 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
281 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
282 Polyhedron
*P
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
283 Polyhedron
*Q
= Domain_Remove_Columns(P
, first
, n
);
284 EVALUE_SET_DOMAIN(e
->x
.p
->arr
[2*i
], Q
);
286 evalue_shift_variables(&e
->x
.p
->arr
[2*i
+1], first
, -n
);
291 static void extract_term_into(const evalue
*src
, int var
, int exp
, evalue
*dst
)
295 if (value_notzero_p(src
->d
) ||
296 src
->x
.p
->type
!= polynomial
||
297 src
->x
.p
->pos
> var
+1) {
299 evalue_copy(dst
, src
);
301 evalue_set_si(dst
, 0, 1);
305 if (src
->x
.p
->pos
== var
+1) {
306 if (src
->x
.p
->size
> exp
)
307 evalue_copy(dst
, &src
->x
.p
->arr
[exp
]);
309 evalue_set_si(dst
, 0, 1);
313 dst
->x
.p
= new_enode(polynomial
, src
->x
.p
->size
, src
->x
.p
->pos
);
314 for (i
= 0; i
< src
->x
.p
->size
; ++i
)
315 extract_term_into(&src
->x
.p
->arr
[i
], var
, exp
,
319 /* Extract the coefficient of var^exp.
321 static evalue
*extract_term(const evalue
*e
, int var
, int exp
)
326 if (EVALUE_IS_ZERO(*e
))
327 return evalue_zero();
329 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
332 res
->x
.p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
333 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
334 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
335 Domain_Copy(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])));
336 extract_term_into(&e
->x
.p
->arr
[2*i
+1], var
, exp
,
337 &res
->x
.p
->arr
[2*i
+1]);
338 reduce_evalue(&res
->x
.p
->arr
[2*i
+1]);
343 /* Insert n free variables at pos (0-based) in the polyhedron P.
345 static Polyhedron
*Polyhedron_Insert_Columns(Polyhedron
*P
, unsigned pos
,
349 unsigned NbConstraints
= 0;
358 assert(pos
<= P
->Dimension
);
360 if (POL_HAS(P
, POL_INEQUALITIES
))
361 NbConstraints
= P
->NbConstraints
;
362 if (POL_HAS(P
, POL_POINTS
))
363 NbRays
= P
->NbRays
+ n
;
365 Q
= Polyhedron_Alloc(P
->Dimension
+n
, NbConstraints
, NbRays
);
366 if (POL_HAS(P
, POL_INEQUALITIES
)) {
368 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
369 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
370 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[i
]+1+pos
+n
,
374 if (POL_HAS(P
, POL_POINTS
)) {
375 Q
->NbBid
= P
->NbBid
+ n
;
376 for (i
= 0; i
< n
; ++i
)
377 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
378 for (i
= 0; i
< P
->NbRays
; ++i
) {
379 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[n
+i
], 1+pos
);
380 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[n
+i
]+1+pos
+n
,
384 POL_SET(Q
, POL_VALID
);
385 if (POL_HAS(P
, POL_INEQUALITIES
))
386 POL_SET(Q
, POL_INEQUALITIES
);
387 if (POL_HAS(P
, POL_POINTS
))
388 POL_SET(Q
, POL_POINTS
);
389 if (POL_HAS(P
, POL_VERTICES
))
390 POL_SET(Q
, POL_VERTICES
);
394 /* Perform summation of e over a list of 1 or more factors F, with context C.
395 * nvar is the total number of variables in the remaining factors.
396 * extra is the number of placeholder parameters introduced in e,
397 * but not (yet) in F or C.
399 * If there is only one factor left, F is intersected with the
400 * context C, the placeholder variables are added, and then
401 * e is summed over the resulting parametric polytope.
403 * If there is more than one factor left, we create two polynomials
404 * in a new placeholder variable (which is placed after the regular
405 * parameters, but before any previously introduced placeholder
406 * variables) that has the factors of the variables in the first
407 * factor of F and the factor of the remaining variables of
408 * each term as its coefficients.
409 * These two polynomials are then summed over their domains
410 * and afterwards the results are combined and the placeholder
411 * variable is removed again.
413 static evalue
*sum_factors(Polyhedron
*F
, Polyhedron
*C
, evalue
*e
,
414 unsigned nvar
, unsigned extra
,
415 struct barvinok_options
*options
)
418 unsigned nparam
= C
->Dimension
;
419 unsigned F_var
= F
->Dimension
- C
->Dimension
;
425 Polyhedron
*CA
= align_context(C
, nvar
+nparam
, options
->MaxRays
);
426 Polyhedron
*P
= DomainIntersection(F
, CA
, options
->MaxRays
);
427 Polyhedron
*Q
= Polyhedron_Insert_Columns(P
, nvar
+nparam
, extra
);
431 evalue
*sum
= sum_base(Q
, e
, nvar
, options
);
436 n
= evalue_count_terms(e
, F_var
, 0);
437 ph
= create_placeholder(n
, nvar
+nparam
);
438 evalue_shift_variables(e
, nvar
+nparam
, 1);
439 evalue_unzip_terms(e
, ph
, F_var
, 0);
440 evalue_shift_variables(e
, nvar
, -(nvar
-F_var
));
441 evalue_reorder_terms(ph
);
442 evalue_shift_variables(ph
, 0, -F_var
);
444 s2
= sum_factors(F
->next
, C
, ph
, nvar
-F_var
, extra
+1, options
);
447 s1
= sum_factors(F
, C
, e
, F_var
, extra
+1, options
);
450 /* remove placeholder "polynomial" */
454 drop_parameters(s2
, nparam
, 1);
459 for (i
= 0; i
< n
; ++i
) {
461 t1
= extract_term(s1
, nparam
, i
);
462 t2
= extract_term(s2
, nparam
, i
);
471 drop_parameters(s
, nparam
, 1);
475 /* Perform summation over a product of factors F, obtained using
476 * variable transformation T from the original problem specification.
478 * We first perform the corresponding transformation on the polynomial E,
479 * compute the common context over all factors and then perform
480 * the actual summation over the factors.
482 static evalue
*sum_product(Polyhedron
*F
, evalue
*E
, Matrix
*T
, unsigned nparam
,
483 struct barvinok_options
*options
)
487 unsigned nvar
= T
->NbRows
;
491 assert(nvar
== T
->NbColumns
);
492 T2
= Matrix_Alloc(nvar
+1, nvar
+1);
493 for (i
= 0; i
< nvar
; ++i
)
494 Vector_Copy(T
->p
[i
], T2
->p
[i
], nvar
);
495 value_set_si(T2
->p
[nvar
][nvar
], 1);
497 transform_polynomial(E
, T2
, NULL
, nvar
, nparam
, nvar
, nparam
);
499 C
= Factor_Context(F
, nparam
, options
->MaxRays
);
500 if (F
->Dimension
== nparam
) {
506 sum
= sum_factors(F
, C
, E
, nvar
, 0, options
);
514 /* Add two constraints corresponding to floor = floor(e/d),
517 * -e + d t + d-1 >= 0
519 * e is assumed to be an affine expression.
521 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
522 struct barvinok_options
*options
)
525 unsigned dim
= P
->Dimension
+1;
526 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
528 Value
*d
= &M
->p
[0][1+nvar
];
529 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
530 value_oppose(*d
, *d
);
531 value_set_si(M
->p
[0][0], 1);
532 value_set_si(M
->p
[1][0], 1);
533 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
534 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
535 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
537 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
538 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
539 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
542 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
547 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
558 /* Compute sum of a step-polynomial over a polytope by grouping
559 * terms containing the same floor-expressions and introducing
560 * new variables for each such expression.
561 * In particular, while there is any floor-expression left,
562 * the step-polynomial is split into a polynomial containing
563 * the expression, which is then converted to a new variable,
564 * and a polynomial not containing the expression.
566 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
567 struct barvinok_options
*options
)
574 while ((floor
= evalue_outer_floor(cur
))) {
577 evalue
*converted_floor
;
579 /* Ignore floors that do not depend on variables. */
580 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
583 converted
= evalue_dup(cur
);
584 converted_floor
= evalue_dup(floor
);
585 evalue_shift_variables(converted
, nvar
, 1);
586 evalue_shift_variables(converted_floor
, nvar
, 1);
587 evalue_replace_floor(converted
, converted_floor
, nvar
);
588 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
589 evalue_free(converted_floor
);
590 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
591 evalue_free(converted
);
593 sum
= evalue_add(t
, sum
);
596 cur
= evalue_dup(cur
);
597 evalue_drop_floor(cur
, floor
);
601 evalue_floor2frac(cur
);
605 if (EVALUE_IS_ZERO(*cur
))
609 unsigned nparam
= P
->Dimension
- nvar
;
610 Polyhedron
*F
= Polyhedron_Factor(P
, nparam
, &T
, options
->MaxRays
);
612 t
= sum_base(P
, cur
, nvar
, options
);
615 cur
= evalue_dup(cur
);
616 t
= sum_product(F
, cur
, T
, nparam
, options
);
623 return evalue_add(t
, sum
);
626 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
627 struct evalue_section_array
*sections
,
628 struct barvinok_options
*options
)
631 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
633 if (options
->summation
== BV_SUM_BERNOULLI
)
634 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
635 else if (options
->summation
== BV_SUM_BOX
)
636 return box_summate(P
, E
, nvar
, options
->MaxRays
);
638 evalue_frac2floor2(E
, 0);
640 return sum_step_polynomial(P
, E
, nvar
, options
);
643 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
646 struct evalue_section_array sections
;
650 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
651 return evalue_dup(e
);
653 assert(value_zero_p(e
->d
));
654 assert(e
->x
.p
->type
== partition
);
656 evalue_section_array_init(§ions
);
659 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
661 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
662 Polyhedron
*next
= D
->next
;
666 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
682 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
685 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
686 options
->MaxRays
= MaxRays
;
687 sum
= barvinok_summate(E
, nvar
, options
);
688 barvinok_options_free(options
);
692 evalue
*esum(evalue
*e
, int nvar
)
695 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
696 sum
= barvinok_summate(e
, nvar
, options
);
697 barvinok_options_free(options
);
701 /* Turn unweighted counting problem into "weighted" counting problem
702 * with weight equal to 1 and call barvinok_summate on this weighted problem.
704 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
705 struct barvinok_options
*options
)
711 if (emptyQ(P
) || emptyQ(C
))
712 return evalue_zero();
714 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
715 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
720 return evalue_zero();
724 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
725 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
726 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
727 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
728 free_evalue_refs(&e
);