1 #include <barvinok/options.h>
6 #include "section_array.h"
8 extern evalue
*evalue_outer_floor(evalue
*e
);
9 extern int evalue_replace_floor(evalue
*e
, const evalue
*floor
, int var
);
10 extern void evalue_drop_floor(evalue
*e
, const evalue
*floor
);
12 #define ALLOC(type) (type*)malloc(sizeof(type))
13 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
15 static evalue
*sum_over_polytope_with_equalities(Polyhedron
*P
, evalue
*E
,
17 struct evalue_section_array
*sections
,
18 struct barvinok_options
*options
)
20 unsigned dim
= P
->Dimension
;
21 unsigned new_dim
, new_nparam
;
22 Matrix
*T
= NULL
, *CP
= NULL
;
32 remove_all_equalities(&P
, NULL
, &CP
, &T
, dim
-nvar
, options
->MaxRays
);
39 new_nparam
= CP
? CP
->NbColumns
-1 : dim
- nvar
;
40 new_dim
= T
? T
->NbColumns
-1 : nvar
+ new_nparam
;
42 /* We can avoid these substitutions if E is a constant */
43 subs
= ALLOCN(evalue
*, dim
);
44 for (j
= 0; j
< nvar
; ++j
) {
46 subs
[j
] = affine2evalue(T
->p
[j
], T
->p
[nvar
+new_nparam
][new_dim
],
49 subs
[j
] = evalue_var(j
);
51 for (j
= 0; j
< dim
-nvar
; ++j
) {
53 subs
[nvar
+j
] = affine2evalue(CP
->p
[j
], CP
->p
[dim
-nvar
][new_nparam
],
56 subs
[nvar
+j
] = evalue_var(j
);
57 evalue_shift_variables(subs
[nvar
+j
], 0, new_dim
-new_nparam
);
61 evalue_substitute(E
, subs
);
64 for (j
= 0; j
< dim
; ++j
)
68 if (new_dim
-new_nparam
> 0) {
69 sum
= barvinok_sum_over_polytope(P
, E
, new_dim
-new_nparam
,
76 sum
->x
.p
= new_enode(partition
, 2, new_dim
);
77 EVALUE_SET_DOMAIN(sum
->x
.p
->arr
[0], P
);
78 value_clear(sum
->x
.p
->arr
[1].d
);
79 sum
->x
.p
->arr
[1] = *E
;
84 evalue_backsubstitute(sum
, CP
, options
->MaxRays
);
94 /* Add two constraints corresponding to floor = floor(e/d),
99 * e is assumed to be an affine expression.
101 Polyhedron
*add_floor_var(Polyhedron
*P
, unsigned nvar
, const evalue
*floor
,
102 struct barvinok_options
*options
)
105 unsigned dim
= P
->Dimension
+1;
106 Matrix
*M
= Matrix_Alloc(P
->NbConstraints
+2, 2+dim
);
108 Value
*d
= &M
->p
[0][1+nvar
];
109 evalue_extract_affine(floor
, M
->p
[0]+1, M
->p
[0]+1+dim
, d
);
110 value_oppose(*d
, *d
);
111 value_set_si(M
->p
[0][0], 1);
112 value_set_si(M
->p
[1][0], 1);
113 Vector_Oppose(M
->p
[0]+1, M
->p
[1]+1, M
->NbColumns
-1);
114 value_subtract(M
->p
[1][1+dim
], M
->p
[1][1+dim
], *d
);
115 value_decrement(M
->p
[1][1+dim
], M
->p
[1][1+dim
]);
117 for (i
= 0; i
< P
->NbConstraints
; ++i
) {
118 Vector_Copy(P
->Constraint
[i
], M
->p
[i
+2], 1+nvar
);
119 Vector_Copy(P
->Constraint
[i
]+1+nvar
, M
->p
[i
+2]+1+nvar
+1, dim
-nvar
-1+1);
122 CP
= Constraints2Polyhedron(M
, options
->MaxRays
);
127 static evalue
*evalue_add(evalue
*a
, evalue
*b
)
138 /* Compute sum of a step-polynomial over a polytope by grouping
139 * terms containing the same floor-expressions and introducing
140 * new variables for each such expression.
141 * In particular, while there is any floor-expression left,
142 * the step-polynomial is split into a polynomial containing
143 * the expression, which is then converted to a new variable,
144 * and a polynomial not containing the expression.
146 static evalue
*sum_step_polynomial(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
147 struct barvinok_options
*options
)
154 while ((floor
= evalue_outer_floor(cur
))) {
157 evalue
*converted_floor
;
159 /* Ignore floors that do not depend on variables. */
160 if (value_notzero_p(floor
->d
) || floor
->x
.p
->pos
>= nvar
+1)
163 converted
= evalue_dup(cur
);
164 converted_floor
= evalue_dup(floor
);
165 evalue_shift_variables(converted
, nvar
, 1);
166 evalue_shift_variables(converted_floor
, nvar
, 1);
167 evalue_replace_floor(converted
, converted_floor
, nvar
);
168 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
169 evalue_free(converted_floor
);
170 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
171 evalue_free(converted
);
173 sum
= evalue_add(t
, sum
);
176 cur
= evalue_dup(cur
);
177 evalue_drop_floor(cur
, floor
);
181 evalue_floor2frac(cur
);
185 if (EVALUE_IS_ZERO(*cur
))
187 else if (options
->summation
== BV_SUM_EULER
)
188 t
= euler_summate(P
, cur
, nvar
, options
);
189 else if (options
->summation
== BV_SUM_LAURENT
)
190 t
= laurent_summate(P
, cur
, nvar
, options
);
197 return evalue_add(t
, sum
);
200 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
201 struct evalue_section_array
*sections
,
202 struct barvinok_options
*options
)
205 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
207 if (options
->summation
== BV_SUM_BERNOULLI
)
208 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
209 else if (options
->summation
== BV_SUM_BARVINOK
)
210 return box_summate(P
, E
, nvar
, options
->MaxRays
);
212 evalue_frac2floor2(E
, 0);
214 return sum_step_polynomial(P
, E
, nvar
, options
);
217 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
220 struct evalue_section_array sections
;
224 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
225 return evalue_dup(e
);
227 assert(value_zero_p(e
->d
));
228 assert(e
->x
.p
->type
== partition
);
230 evalue_section_array_init(§ions
);
233 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
235 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
236 Polyhedron
*next
= D
->next
;
240 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
256 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
259 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
260 options
->MaxRays
= MaxRays
;
261 sum
= barvinok_summate(E
, nvar
, options
);
262 barvinok_options_free(options
);
266 evalue
*esum(evalue
*e
, int nvar
)
269 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
270 sum
= barvinok_summate(e
, nvar
, options
);
271 barvinok_options_free(options
);
275 /* Turn unweighted counting problem into "weighted" counting problem
276 * with weight equal to 1 and call barvinok_summate on this weighted problem.
278 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
279 struct barvinok_options
*options
)
285 if (emptyQ(P
) || emptyQ(C
))
286 return evalue_zero();
288 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
289 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
294 return evalue_zero();
298 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
299 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
300 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
301 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
302 free_evalue_refs(&e
);