add isl_union_map_card
[barvinok.git] / evalue_isl.c
blobd6df6ec6f8c832f0d23c29fcef8ffbcf0ef58662
1 #include <isl_set_polylib.h>
2 #include <isl_constraint.h>
3 #include <barvinok/evalue.h>
5 static __isl_give isl_qpolynomial *extract_base(__isl_take isl_dim *dim,
6 const evalue *e)
8 int i;
9 isl_vec *v;
10 isl_div *div;
11 isl_qpolynomial *base, *c;
12 unsigned nparam;
14 if (!dim)
15 return NULL;
17 if (e->x.p->type == polynomial)
18 return isl_qpolynomial_var(dim, isl_dim_param, e->x.p->pos - 1);
20 nparam = isl_dim_size(dim, isl_dim_param);
21 v = isl_vec_alloc(dim->ctx, 2 + nparam);
22 if (!v)
23 goto error;
25 isl_seq_clr(v->el + 2, nparam);
26 evalue_extract_affine(&e->x.p->arr[0], v->el + 2, &v->el[1], &v->el[0]);
28 div = isl_div_alloc(isl_dim_copy(dim));
29 isl_div_set_constant(div, v->el[1]);
30 isl_div_set_denominator(div, v->el[0]);
32 for (i = 0; i < nparam; ++i)
33 isl_div_set_coefficient(div, isl_dim_param, i, v->el[2 + i]);
35 base = isl_qpolynomial_div(div);
37 if (e->x.p->type == fractional) {
38 base = isl_qpolynomial_neg(base);
40 c = isl_qpolynomial_rat_cst(isl_dim_copy(dim), v->el[1], v->el[0]);
41 base = isl_qpolynomial_add(base, c);
43 for (i = 0; i < nparam; ++i) {
44 isl_qpolynomial *t;
45 c = isl_qpolynomial_rat_cst(isl_dim_copy(dim),
46 v->el[2 + i], v->el[0]);
47 t = isl_qpolynomial_var(isl_dim_copy(dim),
48 isl_dim_param, i);
49 t = isl_qpolynomial_mul(c, t);
50 base = isl_qpolynomial_add(base, t);
54 isl_vec_free(v);
55 isl_dim_free(dim);
57 return base;
58 error:
59 isl_dim_free(dim);
60 return NULL;
63 static int type_offset(enode *p)
65 return p->type == fractional ? 1 :
66 p->type == flooring ? 1 : 0;
69 __isl_give isl_qpolynomial *isl_qpolynomial_from_evalue(__isl_take isl_dim *dim,
70 const evalue *e)
72 int i;
73 isl_qpolynomial *qp;
74 isl_qpolynomial *base;
75 int offset;
77 if (EVALUE_IS_NAN(*e))
78 return isl_qpolynomial_infty(dim);
79 if (value_notzero_p(e->d))
80 return isl_qpolynomial_rat_cst(dim, e->x.n, e->d);
82 offset = type_offset(e->x.p);
84 assert(e->x.p->type == polynomial ||
85 e->x.p->type == flooring ||
86 e->x.p->type == fractional);
87 assert(e->x.p->size >= 1 + offset);
89 base = extract_base(isl_dim_copy(dim), e);
90 qp = isl_qpolynomial_from_evalue(isl_dim_copy(dim),
91 &e->x.p->arr[e->x.p->size - 1]);
93 for (i = e->x.p->size - 2; i >= offset; --i) {
94 qp = isl_qpolynomial_mul(qp, isl_qpolynomial_copy(base));
95 qp = isl_qpolynomial_add(qp,
96 isl_qpolynomial_from_evalue(isl_dim_copy(dim),
97 &e->x.p->arr[i]));
100 isl_qpolynomial_free(base);
101 isl_dim_free(dim);
103 return qp;
106 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
107 const evalue *e);
109 static __isl_give isl_pw_qpolynomial *relation2pwqp(__isl_take isl_set *set,
110 const evalue *e)
112 int i;
113 isl_vec *vec;
114 isl_dim *dim;
115 unsigned nparam;
116 isl_pw_qpolynomial *pwqp;
117 struct isl_constraint *c;
118 struct isl_basic_set *bset;
119 struct isl_set *guard;
120 const evalue *fract;
122 if (!set || !e)
123 goto error;
125 if (e->x.p->size == 1) {
126 dim = isl_set_get_dim(set);
127 isl_set_free(set);
128 return isl_pw_qpolynomial_zero(dim);
131 isl_assert(set->ctx, e->x.p->size > 0, goto error);
132 isl_assert(set->ctx, e->x.p->size <= 3, goto error);
133 isl_assert(set->ctx, value_zero_p(e->x.p->arr[0].d), goto error);
134 isl_assert(set->ctx, e->x.p->arr[0].x.p->type == fractional, goto error);
135 fract = &e->x.p->arr[0];
137 nparam = isl_set_dim(set, isl_dim_param);
138 vec = isl_vec_alloc(set->ctx, 2 + nparam + 1);
139 if (!vec)
140 goto error;
142 isl_seq_clr(vec->el + 2, nparam);
143 evalue_extract_affine(&fract->x.p->arr[0],
144 vec->el + 2, &vec->el[1], &vec->el[0]);
146 dim = isl_set_get_dim(set);
147 dim = isl_dim_add(dim, isl_dim_set, 1);
149 bset = isl_basic_set_universe(dim);
150 c = isl_equality_alloc(isl_dim_copy(dim));
151 isl_int_neg(vec->el[0], vec->el[0]);
152 isl_constraint_set_coefficient(c, isl_dim_set, 0, vec->el[0]);
153 isl_constraint_set_constant(c, vec->el[1]);
154 for (i = 0; i < nparam; ++i)
155 isl_constraint_set_coefficient(c, isl_dim_param, i, vec->el[2+i]);
156 bset = isl_basic_set_add_constraint(bset, c);
157 bset = isl_basic_set_project_out(bset, isl_dim_set, 0, 1);
158 guard = isl_set_from_basic_set(bset);
159 isl_vec_free(vec);
161 pwqp = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
162 isl_set_copy(guard)),
163 &e->x.p->arr[1]);
165 if (e->x.p->size == 3) {
166 isl_pw_qpolynomial *pwqpc;
167 guard = isl_set_complement(guard);
168 pwqpc = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
169 isl_set_copy(guard)),
170 &e->x.p->arr[2]);
171 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqpc);
174 isl_set_free(set);
175 isl_set_free(guard);
177 return pwqp;
178 error:
179 isl_set_free(set);
180 return NULL;
183 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
184 const evalue *e)
186 isl_qpolynomial *qp;
188 if (value_zero_p(e->d) && e->x.p->type == relation)
189 return relation2pwqp(set, e);
191 qp = isl_qpolynomial_from_evalue(isl_set_get_dim(set), e);
193 return isl_pw_qpolynomial_alloc(set, qp);
196 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_evalue(__isl_take isl_dim *dim, const evalue *e)
198 int i;
199 isl_pw_qpolynomial *pwqp;
201 if (!dim || !e)
202 goto error;
203 if (EVALUE_IS_ZERO(*e))
204 return isl_pw_qpolynomial_zero(dim);
206 if (value_notzero_p(e->d)) {
207 isl_set *set = isl_set_universe(isl_dim_copy(dim));
208 isl_qpolynomial *qp = isl_qpolynomial_rat_cst(dim, e->x.n, e->d);
209 return isl_pw_qpolynomial_alloc(set, qp);
212 assert(!EVALUE_IS_NAN(*e));
214 assert(e->x.p->type == partition);
216 pwqp = isl_pw_qpolynomial_zero(isl_dim_copy(dim));
218 for (i = 0; i < e->x.p->size/2; ++i) {
219 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2 * i]);
220 isl_set *set = isl_set_new_from_polylib(D, isl_dim_copy(dim));
221 isl_pw_qpolynomial *pwqp_i;
223 pwqp_i = guarded_evalue2pwqp(set, &e->x.p->arr[2 * i + 1]);
225 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqp_i);
228 isl_dim_free(dim);
230 return pwqp;
231 error:
232 isl_dim_free(dim);
233 return NULL;
236 static evalue *evalue_pow(evalue *e, int exp)
238 evalue *pow;
240 if (exp == 1)
241 return e;
243 pow = evalue_dup(e);
244 while (--exp > 0)
245 emul(e, pow);
247 evalue_free(e);
249 return pow;
252 static evalue *div2evalue(__isl_take isl_div *div)
254 int i;
255 isl_vec *vec;
256 unsigned dim;
257 unsigned nparam;
258 evalue *e;
259 evalue *aff;
261 if (!div)
262 return NULL;
264 dim = isl_div_dim(div, isl_dim_set);
265 nparam = isl_div_dim(div, isl_dim_param);
267 vec = isl_vec_alloc(div->ctx, 1 + dim + nparam + 1);
268 if (!vec)
269 goto error;
270 for (i = 0; i < dim; ++i)
271 isl_div_get_coefficient(div, isl_dim_set, i, &vec->el[1 + i]);
272 for (i = 0; i < nparam; ++i)
273 isl_div_get_coefficient(div, isl_dim_param, i,
274 &vec->el[1 + dim + i]);
275 isl_div_get_denominator(div, &vec->el[0]);
276 isl_div_get_constant(div, &vec->el[1 + dim + nparam]);
278 e = isl_alloc_type(div->ctx, evalue);
279 if (!e)
280 goto error;
281 value_init(e->d);
282 value_set_si(e->d, 0);
283 e->x.p = new_enode(flooring, 3, -1);
284 evalue_set_si(&e->x.p->arr[1], 0, 1);
285 evalue_set_si(&e->x.p->arr[2], 1, 1);
286 value_clear(e->x.p->arr[0].d);
287 aff = affine2evalue(vec->el + 1, vec->el[0], dim + nparam);
288 e->x.p->arr[0] = *aff;
289 free(aff);
290 isl_vec_free(vec);
291 isl_div_free(div);
292 return e;
293 error:
294 isl_vec_free(vec);
295 isl_div_free(div);
296 return NULL;
299 static int add_term(__isl_take isl_term *term, void *user)
301 int i;
302 evalue *sum = (evalue *)user;
303 unsigned nparam;
304 unsigned dim;
305 unsigned n_div;
306 isl_ctx *ctx;
307 isl_int n, d;
308 evalue *e;
310 if (!term)
311 return -1;
313 nparam = isl_term_dim(term, isl_dim_param);
314 dim = isl_term_dim(term, isl_dim_set);
315 n_div = isl_term_dim(term, isl_dim_div);
317 ctx = isl_term_get_ctx(term);
318 e = isl_alloc_type(ctx, evalue);
319 if (!e)
320 goto error;
322 isl_int_init(n);
323 isl_int_init(d);
325 isl_term_get_num(term, &n);
326 isl_term_get_den(term, &d);
327 value_init(e->d);
328 evalue_set(e, n, d);
330 for (i = 0; i < dim; ++i) {
331 evalue *pow;
332 int exp = isl_term_get_exp(term, isl_dim_set, i);
334 if (!exp)
335 continue;
337 pow = evalue_pow(evalue_var(i), exp);
338 emul(pow, e);
339 evalue_free(pow);
342 for (i = 0; i < nparam; ++i) {
343 evalue *pow;
344 int exp = isl_term_get_exp(term, isl_dim_param, i);
346 if (!exp)
347 continue;
349 pow = evalue_pow(evalue_var(dim + i), exp);
350 emul(pow, e);
351 evalue_free(pow);
354 for (i = 0; i < n_div; ++i) {
355 evalue *pow;
356 evalue *floor;
357 isl_div *div;
358 int exp = isl_term_get_exp(term, isl_dim_div, i);
360 if (!exp)
361 continue;
363 div = isl_term_get_div(term, i);
364 floor = div2evalue(div);
365 pow = evalue_pow(floor, exp);
366 emul(pow, e);
367 evalue_free(pow);
370 eadd(e, sum);
371 evalue_free(e);
373 isl_int_clear(n);
374 isl_int_clear(d);
376 isl_term_free(term);
378 return 0;
379 error:
380 isl_term_free(term);
381 return -1;
384 evalue *isl_qpolynomial_to_evalue(__isl_keep isl_qpolynomial *qp)
386 evalue *e;
388 e = evalue_zero();
389 if (!e)
390 return NULL;
392 if (isl_qpolynomial_foreach_term(qp, add_term, e) < 0)
393 goto error;
395 return e;
396 error:
397 evalue_free(e);
398 return NULL;
401 static int add_guarded_qp(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
402 void *user)
404 Polyhedron *D;
405 evalue *e = NULL;
406 evalue *f;
407 evalue *sum = (evalue *)user;
408 unsigned dim;
410 e = isl_alloc_type(set->ctx, evalue);
411 if (!e)
412 goto error;
414 D = isl_set_to_polylib(set);
415 if (!D)
416 goto error;
418 f = isl_qpolynomial_to_evalue(qp);
419 if (!e) {
420 Domain_Free(D);
421 goto error;
424 dim = isl_set_dim(set, isl_dim_param) + isl_set_dim(set, isl_dim_set);
425 value_init(e->d);
426 e->x.p = new_enode(partition, 2, D->Dimension);
427 EVALUE_SET_DOMAIN(e->x.p->arr[0], D);
429 value_clear(e->x.p->arr[1].d);
430 e->x.p->arr[1] = *f;
431 free(f);
433 eadd(e, sum);
434 evalue_free(e);
436 isl_set_free(set);
437 isl_qpolynomial_free(qp);
439 return 0;
440 error:
441 free(e);
442 isl_set_free(set);
443 isl_qpolynomial_free(qp);
444 return -1;
447 evalue *isl_pw_qpolynomial_to_evalue(__isl_keep isl_pw_qpolynomial *pwqp)
449 evalue *e;
451 if (!pwqp)
452 return NULL;
453 e = evalue_zero();
455 if (isl_pw_qpolynomial_foreach_piece(pwqp, add_guarded_qp, e))
456 goto error;
458 return e;
459 error:
460 evalue_free(e);
461 return NULL;