4 #include <isl/union_set.h>
5 #include <isl/union_map.h>
6 #include <isl/polynomial.h>
7 #include <isl_set_polylib.h>
8 #include <barvinok/polylib.h>
9 #include <barvinok/options.h>
10 #include <barvinok/util.h>
11 #include "bernoulli.h"
14 #include "laurent_old.h"
16 #include "section_array.h"
17 #include "remove_equalities.h"
19 extern evalue
*evalue_outer_floor(evalue
*e
);
20 extern int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
);
21 extern void evalue_drop_floor(evalue
*e
, const evalue
*floor
);
23 #define ALLOC(type) (type*)malloc(sizeof(type))
24 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
26 /* Apply the variable transformation specified by T and CP on
27 * the polynomial e. T expresses the old variables in terms
28 * of the new variables (and optionally also the new parameters),
29 * while CP expresses the old parameters in terms of the new
32 static void transform_polynomial(evalue
*E
, Matrix
*T
, Matrix
*CP
,
33 unsigned nvar
, unsigned nparam
,
34 unsigned new_nvar
, unsigned new_nparam
)
39 subs
= ALLOCN(evalue
*, nvar
+nparam
);
41 for (j
= 0; j
< nvar
; ++j
) {
43 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[T
->NbRows
-1][T
->NbColumns
-1],
46 subs
[j
] = evalue_var(j
);
48 for (j
= 0; j
< nparam
; ++j
) {
50 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[nparam
][new_nparam
],
53 subs
[nvar
+j
] = evalue_var(j
);
54 evalue_shift_variables(subs
[nvar
+j
], 0, new_nvar
);
57 evalue_substitute(E
, subs
);
60 for (j
= 0; j
< nvar
+nparam
; ++j
)
65 /* Compute the sum of the quasi-polynomial E
66 * over a 0D (non-empty, but possibly parametric) polytope P.
70 * We simply return a partition evalue with P as domain and E as value.
72 static evalue
*sum_over_polytope_0D(Polyhedron
*P
, evalue
*E
)
78 sum
->x
.p
= new_enode(partition
, 2, P
->Dimension
);
79 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
80 value_clear(sum
->x
.p
->arr
[1].d
);
81 sum
->x
.p
->arr
[1] = *E
;
87 static evalue
*sum_with_equalities(Polyhedron
*P
, evalue
*E
,
88 unsigned nvar
, struct evalue_section_array
*sections
,
89 struct barvinok_options
*options
,
90 evalue
*(*base
)(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
91 struct evalue_section_array
*sections
,
92 struct barvinok_options
*options
))
94 unsigned dim
= P
->Dimension
;
95 unsigned new_dim
, new_nparam
;
96 Matrix
*T
= NULL
, *CP
= NULL
;
100 return evalue_zero();
104 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
108 return evalue_zero();
111 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
112 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
114 /* We can avoid these substitutions if E is a constant */
116 transform_polynomial(E
, T
, CP
, nvar
, dim
-nvar
,
117 new_dim
-new_nparam
, new_nparam
);
119 if (new_dim
-new_nparam
> 0) {
120 sum
= base(P
, E
, new_dim
-new_nparam
, sections
, options
);
124 sum
= sum_over_polytope_0D(P
, E
);
128 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
138 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
139 unsigned nvar
, struct evalue_section_array
*sections
,
140 struct barvinok_options
*options
)
142 return sum_with_equalities(P
, E
, nvar
, sections
, options
,
143 &barvinok_sum_over_polytope
);
146 static evalue
*sum_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
147 struct barvinok_options
*options
);
149 static evalue
*sum_base_wrap(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
150 struct evalue_section_array
*sections
, struct barvinok_options
*options
)
152 return sum_base(P
, E
, nvar
, options
);
155 static evalue
*sum_base_with_equalities(Polyhedron
*P
, evalue
*E
,
156 unsigned nvar
, struct barvinok_options
*options
)
158 return sum_with_equalities(P
, E
, nvar
, NULL
, options
, &sum_base_wrap
);
161 /* The substitutions in sum_step_polynomial may have reintroduced equalities
162 * (in particular, one of the floor expressions may be equal to one of
163 * the variables), so we need to check for them again.
165 static evalue
*sum_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
166 struct barvinok_options
*options
)
169 return sum_base_with_equalities(P
, E
, nvar
, options
);
170 if (options
->summation
== BV_SUM_EULER
)
171 return euler_summate(P
, E
, nvar
, options
);
172 else if (options
->summation
== BV_SUM_LAURENT
)
173 return laurent_summate(P
, E
, nvar
, options
);
174 else if (options
->summation
== BV_SUM_LAURENT_OLD
)
175 return laurent_summate_old(P
, E
, nvar
, options
);
179 /* Count the number of non-zero terms in e when viewed as a polynomial
180 * in only the first nvar variables. "count" is the number counted
183 static int evalue_count_terms(const evalue
*e
, unsigned nvar
, int count
)
187 if (EVALUE_IS_ZERO(*e
))
190 if (value_zero_p(e
->d
))
191 assert(e
->x
.p
->type
== polynomial
);
192 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1)
195 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
196 count
= evalue_count_terms(&e
->x
.p
->arr
[i
], nvar
, count
);
201 /* Create placeholder structure for unzipping.
202 * A "polynomial" is created with size terms in variable pos,
203 * with each term having itself as coefficient.
205 static evalue
*create_placeholder(int size
, int pos
)
208 evalue
*E
= ALLOC(evalue
);
210 E
->x
.p
= new_enode(polynomial
, size
, pos
+1);
211 for (i
= 0; i
< size
; ++i
) {
212 E
->x
.p
->arr
[i
].x
.p
= new_enode(polynomial
, i
+1, pos
+1);
213 for (j
= 0; j
< i
; ++j
)
214 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[j
], 0, 1);
215 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[i
], 1, 1);
220 /* Interchange each non-zero term in e (when viewed as a polynomial
221 * in only the first nvar variables) with a placeholder in ph (created
222 * by create_placeholder), resulting in two polynomials in the
223 * placeholder variable such that for each non-zero term in e
224 * there is a power of the placeholder variable such that the factors
225 * in the first nvar variables form the coefficient of that power in
226 * the first polynomial (e) and the factors in the remaining variables
227 * form the coefficient of that power in the second polynomial (ph).
229 static int evalue_unzip_terms(evalue
*e
, evalue
*ph
, unsigned nvar
, int count
)
233 if (EVALUE_IS_ZERO(*e
))
236 if (value_zero_p(e
->d
))
237 assert(e
->x
.p
->type
== polynomial
);
238 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1) {
240 *e
= ph
->x
.p
->arr
[count
];
241 ph
->x
.p
->arr
[count
] = t
;
245 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
246 count
= evalue_unzip_terms(&e
->x
.p
->arr
[i
], ph
, nvar
, count
);
251 /* Remove n variables at pos (0-based) from the polyhedron P.
252 * Each of these variables is assumed to be completely free,
253 * i.e., there is a line in the polyhedron corresponding to
254 * each of these variables.
256 static Polyhedron
*Polyhedron_Remove_Columns(Polyhedron
*P
, unsigned pos
,
260 unsigned NbConstraints
= 0;
267 assert(pos
<= P
->Dimension
);
269 if (POL_HAS(P
, POL_INEQUALITIES
))
270 NbConstraints
= P
->NbConstraints
;
271 if (POL_HAS(P
, POL_POINTS
))
272 NbRays
= P
->NbRays
- n
;
274 Q
= Polyhedron_Alloc(P
->Dimension
- n
, NbConstraints
, NbRays
);
275 if (POL_HAS(P
, POL_INEQUALITIES
)) {
277 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
278 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
279 Vector_Copy(P
->Constraint
[i
]+1+pos
+n
, Q
->Constraint
[i
]+1+pos
,
283 if (POL_HAS(P
, POL_POINTS
)) {
284 Q
->NbBid
= P
->NbBid
- n
;
285 for (i
= 0; i
< n
; ++i
)
286 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
287 for (i
= 0, j
= 0; i
< P
->NbRays
; ++i
) {
288 int line
= First_Non_Zero(P
->Ray
[i
], 1+P
->Dimension
+1);
290 if (line
-1 >= pos
&& line
-1 < pos
+n
) {
295 assert(i
-j
< Q
->NbRays
);
296 Vector_Copy(P
->Ray
[i
], Q
->Ray
[i
-j
], 1+pos
);
297 Vector_Copy(P
->Ray
[i
]+1+pos
+n
, Q
->Ray
[i
-j
]+1+pos
,
301 POL_SET(Q
, POL_VALID
);
302 if (POL_HAS(P
, POL_INEQUALITIES
))
303 POL_SET(Q
, POL_INEQUALITIES
);
304 if (POL_HAS(P
, POL_POINTS
))
305 POL_SET(Q
, POL_POINTS
);
306 if (POL_HAS(P
, POL_VERTICES
))
307 POL_SET(Q
, POL_VERTICES
);
311 /* Remove n variables at pos (0-based) from the union of polyhedra P.
312 * Each of these variables is assumed to be completely free,
313 * i.e., there is a line in the polyhedron corresponding to
314 * each of these variables.
316 static Polyhedron
*Domain_Remove_Columns(Polyhedron
*P
, unsigned pos
,
320 Polyhedron
**next
= &R
;
322 for (; P
; P
= P
->next
) {
323 *next
= Polyhedron_Remove_Columns(P
, pos
, n
);
324 next
= &(*next
)->next
;
329 /* Drop n parameters starting at first from partition evalue e */
330 static void drop_parameters(evalue
*e
, int first
, int n
)
334 if (EVALUE_IS_ZERO(*e
))
337 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
338 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
339 Polyhedron
*P
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
340 Polyhedron
*Q
= Domain_Remove_Columns(P
, first
, n
);
341 EVALUE_SET_DOMAIN(e
->x
.p
->arr
[2*i
], Q
);
343 evalue_shift_variables(&e
->x
.p
->arr
[2*i
+1], first
, -n
);
348 static void extract_term_into(const evalue
*src
, int var
, int exp
, evalue
*dst
)
352 if (value_notzero_p(src
->d
) ||
353 src
->x
.p
->type
!= polynomial
||
354 src
->x
.p
->pos
> var
+1) {
356 evalue_copy(dst
, src
);
358 evalue_set_si(dst
, 0, 1);
362 if (src
->x
.p
->pos
== var
+1) {
363 if (src
->x
.p
->size
> exp
)
364 evalue_copy(dst
, &src
->x
.p
->arr
[exp
]);
366 evalue_set_si(dst
, 0, 1);
370 dst
->x
.p
= new_enode(polynomial
, src
->x
.p
->size
, src
->x
.p
->pos
);
371 for (i
= 0; i
< src
->x
.p
->size
; ++i
)
372 extract_term_into(&src
->x
.p
->arr
[i
], var
, exp
,
376 /* Extract the coefficient of var^exp.
378 static evalue
*extract_term(const evalue
*e
, int var
, int exp
)
383 if (EVALUE_IS_ZERO(*e
))
384 return evalue_zero();
386 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
389 res
->x
.p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
390 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
391 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
392 Domain_Copy(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])));
393 extract_term_into(&e
->x
.p
->arr
[2*i
+1], var
, exp
,
394 &res
->x
.p
->arr
[2*i
+1]);
395 reduce_evalue(&res
->x
.p
->arr
[2*i
+1]);
400 /* Insert n free variables at pos (0-based) in the polyhedron P.
402 static Polyhedron
*Polyhedron_Insert_Columns(Polyhedron
*P
, unsigned pos
,
406 unsigned NbConstraints
= 0;
415 assert(pos
<= P
->Dimension
);
417 if (POL_HAS(P
, POL_INEQUALITIES
))
418 NbConstraints
= P
->NbConstraints
;
419 if (POL_HAS(P
, POL_POINTS
))
420 NbRays
= P
->NbRays
+ n
;
422 Q
= Polyhedron_Alloc(P
->Dimension
+n
, NbConstraints
, NbRays
);
423 if (POL_HAS(P
, POL_INEQUALITIES
)) {
425 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
426 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
427 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[i
]+1+pos
+n
,
431 if (POL_HAS(P
, POL_POINTS
)) {
432 Q
->NbBid
= P
->NbBid
+ n
;
433 for (i
= 0; i
< n
; ++i
)
434 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
435 for (i
= 0; i
< P
->NbRays
; ++i
) {
436 Vector_Copy(P
->Ray
[i
], Q
->Ray
[n
+i
], 1+pos
);
437 Vector_Copy(P
->Ray
[i
]+1+pos
, Q
->Ray
[n
+i
]+1+pos
+n
,
441 POL_SET(Q
, POL_VALID
);
442 if (POL_HAS(P
, POL_INEQUALITIES
))
443 POL_SET(Q
, POL_INEQUALITIES
);
444 if (POL_HAS(P
, POL_POINTS
))
445 POL_SET(Q
, POL_POINTS
);
446 if (POL_HAS(P
, POL_VERTICES
))
447 POL_SET(Q
, POL_VERTICES
);
451 /* Perform summation of e over a list of 1 or more factors F, with context C.
452 * nvar is the total number of variables in the remaining factors.
453 * extra is the number of placeholder parameters introduced in e,
454 * but not (yet) in F or C.
456 * If there is only one factor left, F is intersected with the
457 * context C, the placeholder variables are added, and then
458 * e is summed over the resulting parametric polytope.
460 * If there is more than one factor left, we create two polynomials
461 * in a new placeholder variable (which is placed after the regular
462 * parameters, but before any previously introduced placeholder
463 * variables) that has the factors of the variables in the first
464 * factor of F and the factor of the remaining variables of
465 * each term as its coefficients.
466 * These two polynomials are then summed over their domains
467 * and afterwards the results are combined and the placeholder
468 * variable is removed again.
470 static evalue
*sum_factors(Polyhedron
*F
, Polyhedron
*C
, evalue
*e
,
471 unsigned nvar
, unsigned extra
,
472 struct barvinok_options
*options
)
474 unsigned nparam
= C
->Dimension
;
475 unsigned F_var
= F
->Dimension
- C
->Dimension
;
481 Polyhedron
*CA
= align_context(C
, nvar
+nparam
, options
->MaxRays
);
482 Polyhedron
*P
= DomainIntersection(F
, CA
, options
->MaxRays
);
483 Polyhedron
*Q
= Polyhedron_Insert_Columns(P
, nvar
+nparam
, extra
);
487 evalue
*sum
= sum_base(Q
, e
, nvar
, options
);
492 n
= evalue_count_terms(e
, F_var
, 0);
493 ph
= create_placeholder(n
, nvar
+nparam
);
494 evalue_shift_variables(e
, nvar
+nparam
, 1);
495 evalue_unzip_terms(e
, ph
, F_var
, 0);
496 evalue_shift_variables(e
, nvar
, -(nvar
-F_var
));
497 evalue_reorder_terms(ph
);
498 evalue_shift_variables(ph
, 0, -F_var
);
500 s2
= sum_factors(F
->next
, C
, ph
, nvar
-F_var
, extra
+1, options
);
503 s1
= sum_factors(F
, C
, e
, F_var
, extra
+1, options
);
506 /* remove placeholder "polynomial" */
510 drop_parameters(s2
, nparam
, 1);
515 for (i
= 0; i
< n
; ++i
) {
517 t1
= extract_term(s1
, nparam
, i
);
518 t2
= extract_term(s2
, nparam
, i
);
527 drop_parameters(s
, nparam
, 1);
531 /* Perform summation over a product of factors F, obtained using
532 * variable transformation T from the original problem specification.
534 * We first perform the corresponding transformation on the polynomial E,
535 * compute the common context over all factors and then perform
536 * the actual summation over the factors.
538 static evalue
*sum_product(Polyhedron
*F
, evalue
*E
, Matrix
*T
, unsigned nparam
,
539 struct barvinok_options
*options
)
543 unsigned nvar
= T
->NbRows
;
547 assert(nvar
== T
->NbColumns
);
548 T2
= Matrix_Alloc(nvar
+1, nvar
+1);
549 for (i
= 0; i
< nvar
; ++i
)
550 Vector_Copy(T
->p
[i
], T2
->p
[i
], nvar
);
551 value_set_si(T2
->p
[nvar
][nvar
], 1);
553 transform_polynomial(E
, T2
, NULL
, nvar
, nparam
, nvar
, nparam
);
555 C
= Factor_Context(F
, nparam
, options
->MaxRays
);
556 if (F
->Dimension
== nparam
) {
562 sum
= sum_factors(F
, C
, E
, nvar
, 0, options
);
570 /* Add two constraints corresponding to floor = floor(e/d),
573 * -e + d t + d-1 >= 0
575 * e is assumed to be an affine expression.
577 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
578 struct barvinok_options
*options
)
581 unsigned dim
= P
->Dimension
+1;
582 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
584 Value
*d
= &M
->p
[0][1+nvar
];
585 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
586 value_oppose(*d
, *d
);
587 value_set_si(M
->p
[0][0], 1);
588 value_set_si(M
->p
[1][0], 1);
589 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
590 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
591 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
593 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
594 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
595 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
598 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
603 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
614 /* Compute sum of a step-polynomial over a polytope by grouping
615 * terms containing the same floor-expressions and introducing
616 * new variables for each such expression.
617 * In particular, while there is any floor-expression left,
618 * the step-polynomial is split into a polynomial containing
619 * the expression, which is then converted to a new variable,
620 * and a polynomial not containing the expression.
622 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
623 struct barvinok_options
*options
)
630 while ((floor
= evalue_outer_floor(cur
))) {
633 evalue
*converted_floor
;
635 /* Ignore floors that do not depend on variables. */
636 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
639 converted
= evalue_dup(cur
);
640 converted_floor
= evalue_dup(floor
);
641 evalue_shift_variables(converted
, nvar
, 1);
642 evalue_shift_variables(converted_floor
, nvar
, 1);
643 evalue_replace_floor(converted
, converted_floor
, nvar
);
644 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
645 evalue_free(converted_floor
);
646 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
647 evalue_free(converted
);
649 sum
= evalue_add(t
, sum
);
652 cur
= evalue_dup(cur
);
653 evalue_drop_floor(cur
, floor
);
657 evalue_floor2frac(cur
);
661 if (EVALUE_IS_ZERO(*cur
))
665 unsigned nparam
= P
->Dimension
- nvar
;
666 Polyhedron
*F
= Polyhedron_Factor(P
, nparam
, &T
, options
->MaxRays
);
668 t
= sum_base(P
, cur
, nvar
, options
);
671 cur
= evalue_dup(cur
);
672 t
= sum_product(F
, cur
, T
, nparam
, options
);
679 return evalue_add(t
, sum
);
682 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
683 struct evalue_section_array
*sections
,
684 struct barvinok_options
*options
)
687 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
690 return sum_over_polytope_0D(Polyhedron_Copy(P
), evalue_dup(E
));
692 if (options
->summation
== BV_SUM_BERNOULLI
)
693 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
694 else if (options
->summation
== BV_SUM_BOX
)
695 return box_summate(P
, E
, nvar
, options
->MaxRays
);
697 evalue_frac2floor2(E
, 0);
699 return sum_step_polynomial(P
, E
, nvar
, options
);
702 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
705 struct evalue_section_array sections
;
709 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
710 return evalue_dup(e
);
712 assert(value_zero_p(e
->d
));
713 assert(e
->x
.p
->type
== partition
);
715 evalue_section_array_init(§ions
);
718 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
720 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
721 Polyhedron
*next
= D
->next
;
725 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
741 static __isl_give isl_pw_qpolynomial
*add_unbounded_guarded_qp(
742 __isl_take isl_pw_qpolynomial
*sum
,
743 __isl_take isl_basic_set
*bset
, __isl_take isl_qpolynomial
*qp
)
747 if (!sum
|| !bset
|| !qp
)
750 zero
= isl_qpolynomial_is_zero(qp
);
757 isl_pw_qpolynomial
*pwqp
;
759 space
= isl_pw_qpolynomial_get_domain_space(sum
);
760 set
= isl_set_from_basic_set(isl_basic_set_copy(bset
));
761 set
= isl_map_domain(isl_map_from_range(set
));
762 set
= isl_set_reset_space(set
, isl_space_copy(space
));
763 pwqp
= isl_pw_qpolynomial_alloc(set
,
764 isl_qpolynomial_nan_on_domain(space
));
765 sum
= isl_pw_qpolynomial_add(sum
, pwqp
);
768 isl_basic_set_free(bset
);
769 isl_qpolynomial_free(qp
);
772 isl_basic_set_free(bset
);
773 isl_qpolynomial_free(qp
);
774 isl_pw_qpolynomial_free(sum
);
778 struct barvinok_summate_data
{
781 isl_pw_qpolynomial
*sum
;
785 struct evalue_section_array sections
;
786 struct barvinok_options
*options
;
789 static isl_stat
add_basic_guarded_qp(__isl_take isl_basic_set
*bset
, void *user
)
791 struct barvinok_summate_data
*data
= user
;
794 isl_pw_qpolynomial
*pwqp
;
796 unsigned nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
800 return isl_stat_error
;
802 bounded
= isl_basic_set_is_bounded(bset
);
807 data
->sum
= add_unbounded_guarded_qp(data
->sum
, bset
,
808 isl_qpolynomial_copy(data
->qp
));
812 space
= isl_space_params(isl_basic_set_get_space(bset
));
814 P
= isl_basic_set_to_polylib(bset
);
815 tmp
= barvinok_sum_over_polytope(P
, data
->e
, nvar
,
816 &data
->sections
, data
->options
);
819 pwqp
= isl_pw_qpolynomial_from_evalue(space
, tmp
);
821 pwqp
= isl_pw_qpolynomial_reset_domain_space(pwqp
,
822 isl_space_domain(isl_space_copy(data
->space
)));
823 data
->sum
= isl_pw_qpolynomial_add(data
->sum
, pwqp
);
825 isl_basic_set_free(bset
);
829 isl_basic_set_free(bset
);
830 return isl_stat_error
;
833 static isl_stat
add_guarded_qp(__isl_take isl_set
*set
,
834 __isl_take isl_qpolynomial
*qp
, void *user
)
837 struct barvinok_summate_data
*data
= user
;
844 if (data
->wrapping
) {
845 unsigned nparam
= isl_set_dim(set
, isl_dim_param
);
846 isl_qpolynomial
*qp2
= isl_qpolynomial_copy(qp
);
847 set
= isl_set_move_dims(set
, isl_dim_param
, nparam
,
848 isl_dim_set
, 0, data
->n_in
);
849 qp2
= isl_qpolynomial_move_dims(qp2
, isl_dim_param
, nparam
,
850 isl_dim_in
, 0, data
->n_in
);
851 data
->e
= isl_qpolynomial_to_evalue(qp2
);
852 isl_qpolynomial_free(qp2
);
854 data
->e
= isl_qpolynomial_to_evalue(qp
);
858 evalue_section_array_init(&data
->sections
);
860 set
= isl_set_make_disjoint(set
);
861 set
= isl_set_compute_divs(set
);
863 r
= isl_set_foreach_basic_set(set
, &add_basic_guarded_qp
, data
);
865 free(data
->sections
.s
);
867 evalue_free(data
->e
);
870 isl_qpolynomial_free(qp
);
875 isl_qpolynomial_free(qp
);
876 return isl_stat_error
;
879 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_sum(
880 __isl_take isl_pw_qpolynomial
*pwqp
)
883 struct barvinok_summate_data data
;
884 int options_allocated
= 0;
894 nvar
= isl_pw_qpolynomial_dim(pwqp
, isl_dim_set
);
896 data
.space
= isl_pw_qpolynomial_get_domain_space(pwqp
);
899 if (isl_space_is_params(data
.space
))
900 isl_die(isl_pw_qpolynomial_get_ctx(pwqp
), isl_error_invalid
,
901 "input polynomial has no domain", goto error
);
902 data
.wrapping
= isl_space_is_wrapping(data
.space
);
904 data
.space
= isl_space_unwrap(data
.space
);
905 data
.n_in
= isl_space_dim(data
.space
, isl_dim_in
);
906 nvar
= isl_space_dim(data
.space
, isl_dim_out
);
910 data
.space
= isl_space_domain(data
.space
);
912 return isl_pw_qpolynomial_reset_domain_space(pwqp
, data
.space
);
914 data
.space
= isl_space_from_domain(data
.space
);
915 data
.space
= isl_space_add_dims(data
.space
, isl_dim_out
, 1);
916 data
.sum
= isl_pw_qpolynomial_zero(isl_space_copy(data
.space
));
918 ctx
= isl_pw_qpolynomial_get_ctx(pwqp
);
919 data
.options
= isl_ctx_peek_barvinok_options(ctx
);
921 data
.options
= barvinok_options_new_with_defaults();
922 options_allocated
= 1;
925 if (isl_pw_qpolynomial_foreach_lifted_piece(pwqp
,
926 add_guarded_qp
, &data
) < 0)
929 if (options_allocated
)
930 barvinok_options_free(data
.options
);
932 isl_space_free(data
.space
);
934 isl_pw_qpolynomial_free(pwqp
);
938 if (options_allocated
)
939 barvinok_options_free(data
.options
);
940 isl_pw_qpolynomial_free(pwqp
);
941 isl_space_free(data
.space
);
942 isl_pw_qpolynomial_free(data
.sum
);
946 static isl_stat
pw_qpolynomial_sum(__isl_take isl_pw_qpolynomial
*pwqp
,
949 isl_union_pw_qpolynomial
**res
= (isl_union_pw_qpolynomial
**)user
;
950 isl_pw_qpolynomial
*sum
;
951 isl_union_pw_qpolynomial
*upwqp
;
953 sum
= isl_pw_qpolynomial_sum(pwqp
);
954 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(sum
);
955 *res
= isl_union_pw_qpolynomial_add(*res
, upwqp
);
960 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_sum(
961 __isl_take isl_union_pw_qpolynomial
*upwqp
)
964 isl_union_pw_qpolynomial
*res
;
966 space
= isl_union_pw_qpolynomial_get_space(upwqp
);
967 res
= isl_union_pw_qpolynomial_zero(space
);
968 if (isl_union_pw_qpolynomial_foreach_pw_qpolynomial(upwqp
,
969 &pw_qpolynomial_sum
, &res
) < 0)
971 isl_union_pw_qpolynomial_free(upwqp
);
975 isl_union_pw_qpolynomial_free(upwqp
);
976 isl_union_pw_qpolynomial_free(res
);
980 static int join_compatible(__isl_keep isl_space
*space1
,
981 __isl_keep isl_space
*space2
)
984 m
= isl_space_has_equal_params(space1
, space2
);
987 return isl_space_tuple_is_equal(space1
, isl_dim_out
,
991 /* Compute the intersection of the range of the map and the domain
992 * of the piecewise quasipolynomial and then sum the associated
993 * quasipolynomial over all elements in this intersection.
995 * We first introduce some unconstrained dimensions in the
996 * piecewise quasipolynomial, intersect the resulting domain
997 * with the wrapped map and then compute the sum.
999 __isl_give isl_pw_qpolynomial
*isl_map_apply_pw_qpolynomial(
1000 __isl_take isl_map
*map
, __isl_take isl_pw_qpolynomial
*pwqp
)
1004 isl_space
*map_space
;
1005 isl_space
*pwqp_space
;
1009 ctx
= isl_map_get_ctx(map
);
1013 map_space
= isl_map_get_space(map
);
1014 pwqp_space
= isl_pw_qpolynomial_get_space(pwqp
);
1015 ok
= join_compatible(map_space
, pwqp_space
);
1016 isl_space_free(map_space
);
1017 isl_space_free(pwqp_space
);
1019 isl_die(ctx
, isl_error_invalid
, "incompatible dimensions",
1022 n_in
= isl_map_dim(map
, isl_dim_in
);
1023 pwqp
= isl_pw_qpolynomial_insert_dims(pwqp
, isl_dim_in
, 0, n_in
);
1025 dom
= isl_map_wrap(map
);
1026 pwqp
= isl_pw_qpolynomial_reset_domain_space(pwqp
,
1027 isl_set_get_space(dom
));
1029 pwqp
= isl_pw_qpolynomial_intersect_domain(pwqp
, dom
);
1030 pwqp
= isl_pw_qpolynomial_sum(pwqp
);
1035 isl_pw_qpolynomial_free(pwqp
);
1039 __isl_give isl_pw_qpolynomial
*isl_set_apply_pw_qpolynomial(
1040 __isl_take isl_set
*set
, __isl_take isl_pw_qpolynomial
*pwqp
)
1044 map
= isl_map_from_range(set
);
1045 pwqp
= isl_map_apply_pw_qpolynomial(map
, pwqp
);
1046 pwqp
= isl_pw_qpolynomial_project_domain_on_params(pwqp
);
1050 struct barvinok_apply_data
{
1051 isl_union_pw_qpolynomial
*upwqp
;
1052 isl_union_pw_qpolynomial
*res
;
1056 static isl_stat
pw_qpolynomial_apply(__isl_take isl_pw_qpolynomial
*pwqp
,
1059 isl_space
*map_space
;
1060 isl_space
*pwqp_space
;
1061 struct barvinok_apply_data
*data
= user
;
1064 map_space
= isl_map_get_space(data
->map
);
1065 pwqp_space
= isl_pw_qpolynomial_get_space(pwqp
);
1066 ok
= join_compatible(map_space
, pwqp_space
);
1067 isl_space_free(map_space
);
1068 isl_space_free(pwqp_space
);
1071 isl_union_pw_qpolynomial
*upwqp
;
1073 pwqp
= isl_map_apply_pw_qpolynomial(isl_map_copy(data
->map
),
1075 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1076 data
->res
= isl_union_pw_qpolynomial_add(data
->res
, upwqp
);
1078 isl_pw_qpolynomial_free(pwqp
);
1083 static isl_stat
map_apply(__isl_take isl_map
*map
, void *user
)
1085 struct barvinok_apply_data
*data
= user
;
1089 r
= isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data
->upwqp
,
1090 &pw_qpolynomial_apply
, data
);
1096 __isl_give isl_union_pw_qpolynomial
*isl_union_map_apply_union_pw_qpolynomial(
1097 __isl_take isl_union_map
*umap
,
1098 __isl_take isl_union_pw_qpolynomial
*upwqp
)
1101 struct barvinok_apply_data data
;
1103 upwqp
= isl_union_pw_qpolynomial_align_params(upwqp
,
1104 isl_union_map_get_space(umap
));
1105 umap
= isl_union_map_align_params(umap
,
1106 isl_union_pw_qpolynomial_get_space(upwqp
));
1109 space
= isl_union_pw_qpolynomial_get_space(upwqp
);
1110 data
.res
= isl_union_pw_qpolynomial_zero(space
);
1111 if (isl_union_map_foreach_map(umap
, &map_apply
, &data
) < 0)
1114 isl_union_map_free(umap
);
1115 isl_union_pw_qpolynomial_free(upwqp
);
1119 isl_union_map_free(umap
);
1120 isl_union_pw_qpolynomial_free(upwqp
);
1121 isl_union_pw_qpolynomial_free(data
.res
);
1125 struct barvinok_apply_set_data
{
1126 isl_union_pw_qpolynomial
*upwqp
;
1127 isl_union_pw_qpolynomial
*res
;
1131 static isl_stat
pw_qpolynomial_apply_set(__isl_take isl_pw_qpolynomial
*pwqp
,
1134 isl_space
*set_space
;
1135 isl_space
*pwqp_space
;
1136 struct barvinok_apply_set_data
*data
= user
;
1139 set_space
= isl_set_get_space(data
->set
);
1140 pwqp_space
= isl_pw_qpolynomial_get_space(pwqp
);
1141 ok
= join_compatible(set_space
, pwqp_space
);
1142 isl_space_free(set_space
);
1143 isl_space_free(pwqp_space
);
1146 isl_union_pw_qpolynomial
*upwqp
;
1148 pwqp
= isl_set_apply_pw_qpolynomial(isl_set_copy(data
->set
),
1150 upwqp
= isl_union_pw_qpolynomial_from_pw_qpolynomial(pwqp
);
1151 data
->res
= isl_union_pw_qpolynomial_add(data
->res
, upwqp
);
1153 isl_pw_qpolynomial_free(pwqp
);
1158 static isl_stat
set_apply(__isl_take isl_set
*set
, void *user
)
1160 struct barvinok_apply_set_data
*data
= user
;
1164 r
= isl_union_pw_qpolynomial_foreach_pw_qpolynomial(data
->upwqp
,
1165 &pw_qpolynomial_apply_set
, data
);
1171 __isl_give isl_union_pw_qpolynomial
*isl_union_set_apply_union_pw_qpolynomial(
1172 __isl_take isl_union_set
*uset
,
1173 __isl_take isl_union_pw_qpolynomial
*upwqp
)
1176 struct barvinok_apply_set_data data
;
1178 upwqp
= isl_union_pw_qpolynomial_align_params(upwqp
,
1179 isl_union_set_get_space(uset
));
1180 uset
= isl_union_set_align_params(uset
,
1181 isl_union_pw_qpolynomial_get_space(upwqp
));
1184 space
= isl_union_pw_qpolynomial_get_space(upwqp
);
1185 data
.res
= isl_union_pw_qpolynomial_zero(space
);
1186 if (isl_union_set_foreach_set(uset
, &set_apply
, &data
) < 0)
1189 isl_union_set_free(uset
);
1190 isl_union_pw_qpolynomial_free(upwqp
);
1194 isl_union_set_free(uset
);
1195 isl_union_pw_qpolynomial_free(upwqp
);
1196 isl_union_pw_qpolynomial_free(data
.res
);
1200 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
1203 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
1204 options
->MaxRays
= MaxRays
;
1205 sum
= barvinok_summate(E
, nvar
, options
);
1206 barvinok_options_free(options
);
1210 evalue
*esum(evalue
*e
, int nvar
)
1213 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
1214 sum
= barvinok_summate(e
, nvar
, options
);
1215 barvinok_options_free(options
);
1219 /* Turn unweighted counting problem into "weighted" counting problem
1220 * with weight equal to 1 and call barvinok_summate on this weighted problem.
1222 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
1223 struct barvinok_options
*options
)
1229 if (emptyQ(P
) || emptyQ(C
))
1230 return evalue_zero();
1232 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
1233 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
1238 return evalue_zero();
1242 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
1243 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
1244 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
1245 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
1246 free_evalue_refs(&e
);