1 #include <barvinok/options.h>
2 #include <barvinok/util.h>
6 #include "laurent_old.h"
8 #include "section_array.h"
9 #include "remove_equalities.h"
11 extern evalue
*evalue_outer_floor(evalue
*e
);
12 extern int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
);
13 extern void evalue_drop_floor(evalue
*e
, const evalue
*floor
);
15 #define ALLOC(type) (type*)malloc(sizeof(type))
16 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
18 /* Apply the variable transformation specified by T and CP on
19 * the polynomial e. T expresses the old variables in terms
20 * of the new variables (and optionally also the new parameters),
21 * while CP expresses the old parameters in terms of the new
24 static void transform_polynomial(evalue
*E
, Matrix
*T
, Matrix
*CP
,
25 unsigned nvar
, unsigned nparam
,
26 unsigned new_nvar
, unsigned new_nparam
)
31 subs
= ALLOCN(evalue
*, nvar
+nparam
);
33 for (j
= 0; j
< nvar
; ++j
) {
35 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[T
->NbRows
-1][T
->NbColumns
-1],
38 subs
[j
] = evalue_var(j
);
40 for (j
= 0; j
< nparam
; ++j
) {
42 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[nparam
][new_nparam
],
45 subs
[nvar
+j
] = evalue_var(j
);
46 evalue_shift_variables(subs
[nvar
+j
], 0, new_nvar
);
49 evalue_substitute(E
, subs
);
52 for (j
= 0; j
< nvar
+nparam
; ++j
)
57 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
59 struct evalue_section_array
*sections
,
60 struct barvinok_options
*options
)
62 unsigned dim
= P
->Dimension
;
63 unsigned new_dim
, new_nparam
;
64 Matrix
*T
= NULL
, *CP
= NULL
;
72 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
79 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
80 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
82 /* We can avoid these substitutions if E is a constant */
84 transform_polynomial(E
, T
, CP
, nvar
, dim
-nvar
,
85 new_dim
-new_nparam
, new_nparam
);
87 if (new_dim
-new_nparam
> 0) {
88 sum
= barvinok_sum_over_polytope(P
, E
, new_dim
-new_nparam
,
95 sum
->x
.p
= new_enode(partition
, 2, new_dim
);
96 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
97 value_clear(sum
->x
.p
->arr
[1].d
);
98 sum
->x
.p
->arr
[1] = *E
;
103 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
113 static evalue
*sum_base(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
114 struct barvinok_options
*options
)
116 if (options
->summation
== BV_SUM_EULER
)
117 return euler_summate(P
, E
, nvar
, options
);
118 else if (options
->summation
== BV_SUM_LAURENT
)
119 return laurent_summate(P
, E
, nvar
, options
);
120 else if (options
->summation
== BV_SUM_LAURENT_OLD
)
121 return laurent_summate_old(P
, E
, nvar
, options
);
125 /* Count the number of non-zero terms in e when viewed as a polynomial
126 * in only the first nvar variables. "count" is the number counted
129 static int evalue_count_terms(const evalue
*e
, unsigned nvar
, int count
)
133 if (EVALUE_IS_ZERO(*e
))
136 if (value_zero_p(e
->d
))
137 assert(e
->x
.p
->type
== polynomial
);
138 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1)
141 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
142 count
= evalue_count_terms(&e
->x
.p
->arr
[i
], nvar
, count
);
147 /* Create placeholder structure for unzipping.
148 * A "polynomial" is created with size terms in variable pos,
149 * with each term having itself as coefficient.
151 static evalue
*create_placeholder(int size
, int pos
)
154 evalue
*E
= ALLOC(evalue
);
156 E
->x
.p
= new_enode(polynomial
, size
, pos
+1);
157 for (i
= 0; i
< size
; ++i
) {
158 E
->x
.p
->arr
[i
].x
.p
= new_enode(polynomial
, i
+1, pos
+1);
159 for (j
= 0; j
< i
; ++j
)
160 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[j
], 0, 1);
161 evalue_set_si(&E
->x
.p
->arr
[i
].x
.p
->arr
[i
], 1, 1);
166 /* Interchange each non-zero term in e (when viewed as a polynomial
167 * in only the first nvar variables) with a placeholder in ph (created
168 * by create_placeholder), resulting in two polynomials in the
169 * placeholder variable such that for each non-zero term in e
170 * there is a power of the placeholder variable such that the factors
171 * in the first nvar variables form the coefficient of that power in
172 * the first polynomial (e) and the factors in the remaining variables
173 * form the coefficient of that power in the second polynomial (ph).
175 static int evalue_unzip_terms(evalue
*e
, evalue
*ph
, unsigned nvar
, int count
)
179 if (EVALUE_IS_ZERO(*e
))
182 if (value_zero_p(e
->d
))
183 assert(e
->x
.p
->type
== polynomial
);
184 if (value_notzero_p(e
->d
) || e
->x
.p
->pos
>= nvar
+1) {
186 *e
= ph
->x
.p
->arr
[count
];
187 ph
->x
.p
->arr
[count
] = t
;
191 for (i
= 0; i
< e
->x
.p
->size
; ++i
)
192 count
= evalue_unzip_terms(&e
->x
.p
->arr
[i
], ph
, nvar
, count
);
197 /* Remove n variables at pos (0-based) from the polyhedron P.
198 * Each of these variables is assumed to be completely free,
199 * i.e., there is a line in the polyhedron corresponding to
200 * each of these variables.
202 static Polyhedron
*Polyhedron_Remove_Columns(Polyhedron
*P
, unsigned pos
,
206 unsigned NbConstraints
= 0;
213 assert(pos
<= P
->Dimension
);
215 if (POL_HAS(P
, POL_INEQUALITIES
))
216 NbConstraints
= P
->NbConstraints
;
217 if (POL_HAS(P
, POL_POINTS
))
218 NbRays
= P
->NbRays
- n
;
220 Q
= Polyhedron_Alloc(P
->Dimension
- n
, NbConstraints
, NbRays
);
221 if (POL_HAS(P
, POL_INEQUALITIES
)) {
223 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
224 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
225 Vector_Copy(P
->Constraint
[i
]+1+pos
+n
, Q
->Constraint
[i
]+1+pos
,
229 if (POL_HAS(P
, POL_POINTS
)) {
230 Q
->NbBid
= P
->NbBid
- n
;
231 for (i
= 0; i
< n
; ++i
)
232 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
233 for (i
= 0, j
= 0; i
< P
->NbRays
; ++i
) {
234 int line
= First_Non_Zero(P
->Ray
[i
], 1+P
->Dimension
+1);
236 if (line
-1 >= pos
&& line
-1 < pos
+n
) {
241 assert(i
-j
< Q
->NbRays
);
242 Vector_Copy(P
->Ray
[i
], Q
->Ray
[i
-j
], 1+pos
);
243 Vector_Copy(P
->Ray
[i
]+1+pos
+n
, Q
->Ray
[i
-j
]+1+pos
,
247 POL_SET(Q
, POL_VALID
);
248 if (POL_HAS(P
, POL_INEQUALITIES
))
249 POL_SET(Q
, POL_INEQUALITIES
);
250 if (POL_HAS(P
, POL_POINTS
))
251 POL_SET(Q
, POL_POINTS
);
252 if (POL_HAS(P
, POL_VERTICES
))
253 POL_SET(Q
, POL_VERTICES
);
257 /* Remove n variables at pos (0-based) from the union of polyhedra P.
258 * Each of these variables is assumed to be completely free,
259 * i.e., there is a line in the polyhedron corresponding to
260 * each of these variables.
262 static Polyhedron
*Domain_Remove_Columns(Polyhedron
*P
, unsigned pos
,
266 Polyhedron
**next
= &R
;
268 for (; P
; P
= P
->next
) {
269 *next
= Polyhedron_Remove_Columns(P
, pos
, n
);
270 next
= &(*next
)->next
;
275 /* Drop n parameters starting at first from partition evalue e */
276 static void drop_parameters(evalue
*e
, int first
, int n
)
280 if (EVALUE_IS_ZERO(*e
))
283 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
284 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
285 Polyhedron
*P
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]);
286 Polyhedron
*Q
= Domain_Remove_Columns(P
, first
, n
);
287 EVALUE_SET_DOMAIN(e
->x
.p
->arr
[2*i
], Q
);
289 evalue_shift_variables(&e
->x
.p
->arr
[2*i
+1], first
, -n
);
294 static void extract_term_into(const evalue
*src
, int var
, int exp
, evalue
*dst
)
298 if (value_notzero_p(src
->d
) ||
299 src
->x
.p
->type
!= polynomial
||
300 src
->x
.p
->pos
> var
+1) {
302 evalue_copy(dst
, src
);
304 evalue_set_si(dst
, 0, 1);
308 if (src
->x
.p
->pos
== var
+1) {
309 if (src
->x
.p
->size
> exp
)
310 evalue_copy(dst
, &src
->x
.p
->arr
[exp
]);
312 evalue_set_si(dst
, 0, 1);
316 dst
->x
.p
= new_enode(polynomial
, src
->x
.p
->size
, src
->x
.p
->pos
);
317 for (i
= 0; i
< src
->x
.p
->size
; ++i
)
318 extract_term_into(&src
->x
.p
->arr
[i
], var
, exp
,
322 /* Extract the coefficient of var^exp.
324 static evalue
*extract_term(const evalue
*e
, int var
, int exp
)
329 if (EVALUE_IS_ZERO(*e
))
330 return evalue_zero();
332 assert(value_zero_p(e
->d
) && e
->x
.p
->type
== partition
);
335 res
->x
.p
= new_enode(partition
, e
->x
.p
->size
, e
->x
.p
->pos
);
336 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
337 EVALUE_SET_DOMAIN(res
->x
.p
->arr
[2*i
],
338 Domain_Copy(EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
])));
339 extract_term_into(&e
->x
.p
->arr
[2*i
+1], var
, exp
,
340 &res
->x
.p
->arr
[2*i
+1]);
341 reduce_evalue(&res
->x
.p
->arr
[2*i
+1]);
346 /* Insert n free variables at pos (0-based) in the polyhedron P.
348 static Polyhedron
*Polyhedron_Insert_Columns(Polyhedron
*P
, unsigned pos
,
352 unsigned NbConstraints
= 0;
361 assert(pos
<= P
->Dimension
);
363 if (POL_HAS(P
, POL_INEQUALITIES
))
364 NbConstraints
= P
->NbConstraints
;
365 if (POL_HAS(P
, POL_POINTS
))
366 NbRays
= P
->NbRays
+ n
;
368 Q
= Polyhedron_Alloc(P
->Dimension
+n
, NbConstraints
, NbRays
);
369 if (POL_HAS(P
, POL_INEQUALITIES
)) {
371 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
372 Vector_Copy(P
->Constraint
[i
], Q
->Constraint
[i
], 1+pos
);
373 Vector_Copy(P
->Constraint
[i
]+1+pos
, Q
->Constraint
[i
]+1+pos
+n
,
377 if (POL_HAS(P
, POL_POINTS
)) {
378 Q
->NbBid
= P
->NbBid
+ n
;
379 for (i
= 0; i
< n
; ++i
)
380 value_set_si(Q
->Ray
[i
][1+pos
+i
], 1);
381 for (i
= 0; i
< P
->NbRays
; ++i
) {
382 Vector_Copy(P
->Ray
[i
], Q
->Ray
[n
+i
], 1+pos
);
383 Vector_Copy(P
->Ray
[i
]+1+pos
, Q
->Ray
[n
+i
]+1+pos
+n
,
387 POL_SET(Q
, POL_VALID
);
388 if (POL_HAS(P
, POL_INEQUALITIES
))
389 POL_SET(Q
, POL_INEQUALITIES
);
390 if (POL_HAS(P
, POL_POINTS
))
391 POL_SET(Q
, POL_POINTS
);
392 if (POL_HAS(P
, POL_VERTICES
))
393 POL_SET(Q
, POL_VERTICES
);
397 /* Perform summation of e over a list of 1 or more factors F, with context C.
398 * nvar is the total number of variables in the remaining factors.
399 * extra is the number of placeholder parameters introduced in e,
400 * but not (yet) in F or C.
402 * If there is only one factor left, F is intersected with the
403 * context C, the placeholder variables are added, and then
404 * e is summed over the resulting parametric polytope.
406 * If there is more than one factor left, we create two polynomials
407 * in a new placeholder variable (which is placed after the regular
408 * parameters, but before any previously introduced placeholder
409 * variables) that has the factors of the variables in the first
410 * factor of F and the factor of the remaining variables of
411 * each term as its coefficients.
412 * These two polynomials are then summed over their domains
413 * and afterwards the results are combined and the placeholder
414 * variable is removed again.
416 static evalue
*sum_factors(Polyhedron
*F
, Polyhedron
*C
, evalue
*e
,
417 unsigned nvar
, unsigned extra
,
418 struct barvinok_options
*options
)
421 unsigned nparam
= C
->Dimension
;
422 unsigned F_var
= F
->Dimension
- C
->Dimension
;
428 Polyhedron
*CA
= align_context(C
, nvar
+nparam
, options
->MaxRays
);
429 Polyhedron
*P
= DomainIntersection(F
, CA
, options
->MaxRays
);
430 Polyhedron
*Q
= Polyhedron_Insert_Columns(P
, nvar
+nparam
, extra
);
434 evalue
*sum
= sum_base(Q
, e
, nvar
, options
);
439 n
= evalue_count_terms(e
, F_var
, 0);
440 ph
= create_placeholder(n
, nvar
+nparam
);
441 evalue_shift_variables(e
, nvar
+nparam
, 1);
442 evalue_unzip_terms(e
, ph
, F_var
, 0);
443 evalue_shift_variables(e
, nvar
, -(nvar
-F_var
));
444 evalue_reorder_terms(ph
);
445 evalue_shift_variables(ph
, 0, -F_var
);
447 s2
= sum_factors(F
->next
, C
, ph
, nvar
-F_var
, extra
+1, options
);
450 s1
= sum_factors(F
, C
, e
, F_var
, extra
+1, options
);
453 /* remove placeholder "polynomial" */
457 drop_parameters(s2
, nparam
, 1);
462 for (i
= 0; i
< n
; ++i
) {
464 t1
= extract_term(s1
, nparam
, i
);
465 t2
= extract_term(s2
, nparam
, i
);
474 drop_parameters(s
, nparam
, 1);
478 /* Perform summation over a product of factors F, obtained using
479 * variable transformation T from the original problem specification.
481 * We first perform the corresponding transformation on the polynomial E,
482 * compute the common context over all factors and then perform
483 * the actual summation over the factors.
485 static evalue
*sum_product(Polyhedron
*F
, evalue
*E
, Matrix
*T
, unsigned nparam
,
486 struct barvinok_options
*options
)
490 unsigned nvar
= T
->NbRows
;
494 assert(nvar
== T
->NbColumns
);
495 T2
= Matrix_Alloc(nvar
+1, nvar
+1);
496 for (i
= 0; i
< nvar
; ++i
)
497 Vector_Copy(T
->p
[i
], T2
->p
[i
], nvar
);
498 value_set_si(T2
->p
[nvar
][nvar
], 1);
500 transform_polynomial(E
, T2
, NULL
, nvar
, nparam
, nvar
, nparam
);
502 C
= Factor_Context(F
, nparam
, options
->MaxRays
);
503 if (F
->Dimension
== nparam
) {
509 sum
= sum_factors(F
, C
, E
, nvar
, 0, options
);
517 /* Add two constraints corresponding to floor = floor(e/d),
520 * -e + d t + d-1 >= 0
522 * e is assumed to be an affine expression.
524 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
525 struct barvinok_options
*options
)
528 unsigned dim
= P
->Dimension
+1;
529 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
531 Value
*d
= &M
->p
[0][1+nvar
];
532 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
533 value_oppose(*d
, *d
);
534 value_set_si(M
->p
[0][0], 1);
535 value_set_si(M
->p
[1][0], 1);
536 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
537 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
538 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
540 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
541 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
542 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
545 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
550 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
561 /* Compute sum of a step-polynomial over a polytope by grouping
562 * terms containing the same floor-expressions and introducing
563 * new variables for each such expression.
564 * In particular, while there is any floor-expression left,
565 * the step-polynomial is split into a polynomial containing
566 * the expression, which is then converted to a new variable,
567 * and a polynomial not containing the expression.
569 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
570 struct barvinok_options
*options
)
577 while ((floor
= evalue_outer_floor(cur
))) {
580 evalue
*converted_floor
;
582 /* Ignore floors that do not depend on variables. */
583 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
586 converted
= evalue_dup(cur
);
587 converted_floor
= evalue_dup(floor
);
588 evalue_shift_variables(converted
, nvar
, 1);
589 evalue_shift_variables(converted_floor
, nvar
, 1);
590 evalue_replace_floor(converted
, converted_floor
, nvar
);
591 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
592 evalue_free(converted_floor
);
593 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
594 evalue_free(converted
);
596 sum
= evalue_add(t
, sum
);
599 cur
= evalue_dup(cur
);
600 evalue_drop_floor(cur
, floor
);
604 evalue_floor2frac(cur
);
608 if (EVALUE_IS_ZERO(*cur
))
612 unsigned nparam
= P
->Dimension
- nvar
;
613 Polyhedron
*F
= Polyhedron_Factor(P
, nparam
, &T
, options
->MaxRays
);
615 t
= sum_base(P
, cur
, nvar
, options
);
618 cur
= evalue_dup(cur
);
619 t
= sum_product(F
, cur
, T
, nparam
, options
);
626 return evalue_add(t
, sum
);
629 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
630 struct evalue_section_array
*sections
,
631 struct barvinok_options
*options
)
634 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
636 if (options
->summation
== BV_SUM_BERNOULLI
)
637 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
638 else if (options
->summation
== BV_SUM_BOX
)
639 return box_summate(P
, E
, nvar
, options
->MaxRays
);
641 evalue_frac2floor2(E
, 0);
643 return sum_step_polynomial(P
, E
, nvar
, options
);
646 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
649 struct evalue_section_array sections
;
653 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
654 return evalue_dup(e
);
656 assert(value_zero_p(e
->d
));
657 assert(e
->x
.p
->type
== partition
);
659 evalue_section_array_init(§ions
);
662 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
664 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
665 Polyhedron
*next
= D
->next
;
669 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
685 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
688 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
689 options
->MaxRays
= MaxRays
;
690 sum
= barvinok_summate(E
, nvar
, options
);
691 barvinok_options_free(options
);
695 evalue
*esum(evalue
*e
, int nvar
)
698 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
699 sum
= barvinok_summate(e
, nvar
, options
);
700 barvinok_options_free(options
);
704 /* Turn unweighted counting problem into "weighted" counting problem
705 * with weight equal to 1 and call barvinok_summate on this weighted problem.
707 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
708 struct barvinok_options
*options
)
714 if (emptyQ(P
) || emptyQ(C
))
715 return evalue_zero();
717 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
718 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
723 return evalue_zero();
727 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
728 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
729 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
730 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
731 free_evalue_refs(&e
);