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
))) {
156 evalue
*converted
= evalue_dup(cur
);
157 evalue
*converted_floor
= evalue_dup(floor
);
159 evalue_shift_variables(converted
, nvar
, 1);
160 evalue_shift_variables(converted_floor
, nvar
, 1);
161 evalue_replace_floor(converted
, converted_floor
, nvar
);
162 CP
= add_floor_var(P
, nvar
, converted_floor
, options
);
163 evalue_free(converted_floor
);
164 t
= sum_step_polynomial(CP
, converted
, nvar
+1, options
);
165 evalue_free(converted
);
167 sum
= evalue_add(t
, sum
);
170 cur
= evalue_dup(cur
);
171 evalue_drop_floor(cur
, floor
);
175 if (EVALUE_IS_ZERO(*cur
))
177 else if (options
->summation
== BV_SUM_EULER
)
178 t
= euler_summate(P
, cur
, nvar
, options
);
179 else if (options
->summation
== BV_SUM_LAURENT
)
180 t
= laurent_summate(P
, cur
, nvar
, options
);
187 return evalue_add(t
, sum
);
190 evalue
*barvinok_sum_over_polytope(Polyhedron
*P
, evalue
*E
, unsigned nvar
,
191 struct evalue_section_array
*sections
,
192 struct barvinok_options
*options
)
195 return sum_over_polytope_with_equalities(P
, E
, nvar
, sections
, options
);
197 if (options
->summation
== BV_SUM_BERNOULLI
)
198 return bernoulli_summate(P
, E
, nvar
, sections
, options
);
199 else if (options
->summation
== BV_SUM_BARVINOK
)
200 return box_summate(P
, E
, nvar
, options
->MaxRays
);
202 evalue_frac2floor2(E
, 0);
204 return sum_step_polynomial(P
, E
, nvar
, options
);
207 evalue
*barvinok_summate(evalue
*e
, int nvar
, struct barvinok_options
*options
)
210 struct evalue_section_array sections
;
214 if (nvar
== 0 || EVALUE_IS_ZERO(*e
))
215 return evalue_dup(e
);
217 assert(value_zero_p(e
->d
));
218 assert(e
->x
.p
->type
== partition
);
220 evalue_section_array_init(§ions
);
223 for (i
= 0; i
< e
->x
.p
->size
/2; ++i
) {
225 for (D
= EVALUE_DOMAIN(e
->x
.p
->arr
[2*i
]); D
; D
= D
->next
) {
226 Polyhedron
*next
= D
->next
;
230 tmp
= barvinok_sum_over_polytope(D
, &e
->x
.p
->arr
[2*i
+1], nvar
,
246 evalue
*evalue_sum(evalue
*E
, int nvar
, unsigned MaxRays
)
249 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
250 options
->MaxRays
= MaxRays
;
251 sum
= barvinok_summate(E
, nvar
, options
);
252 barvinok_options_free(options
);
256 evalue
*esum(evalue
*e
, int nvar
)
259 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
260 sum
= barvinok_summate(e
, nvar
, options
);
261 barvinok_options_free(options
);
265 /* Turn unweighted counting problem into "weighted" counting problem
266 * with weight equal to 1 and call barvinok_summate on this weighted problem.
268 evalue
*barvinok_summate_unweighted(Polyhedron
*P
, Polyhedron
*C
,
269 struct barvinok_options
*options
)
275 if (emptyQ(P
) || emptyQ(C
))
276 return evalue_zero();
278 CA
= align_context(C
, P
->Dimension
, options
->MaxRays
);
279 D
= DomainIntersection(P
, CA
, options
->MaxRays
);
284 return evalue_zero();
288 e
.x
.p
= new_enode(partition
, 2, P
->Dimension
);
289 EVALUE_SET_DOMAIN(e
.x
.p
->arr
[0], D
);
290 evalue_set_si(&e
.x
.p
->arr
[1], 1, 1);
291 sum
= barvinok_summate(&e
, P
->Dimension
- C
->Dimension
, options
);
292 free_evalue_refs(&e
);