assume NTL has been compiled in ISO mode
[barvinok.git] / evalue_isl.c
blobc710f803e556d813da96d20978c66af1b7f31631
1 #include <isl/aff.h>
2 #include <isl_set_polylib.h>
3 #include <isl/constraint.h>
4 #include <isl/val_gmp.h>
5 #include <barvinok/evalue.h>
7 static __isl_give isl_qpolynomial *extract_base(__isl_take isl_space *dim,
8 const evalue *e)
10 int i;
11 Vector *v;
12 isl_ctx *ctx;
13 isl_local_space *ls;
14 isl_aff *aff;
15 isl_val *val;
16 isl_qpolynomial *base, *c;
17 unsigned nparam;
19 if (!dim)
20 return NULL;
22 if (e->x.p->type == polynomial)
23 return isl_qpolynomial_var_on_domain(dim, isl_dim_param, e->x.p->pos - 1);
25 ctx = isl_space_get_ctx(dim);
26 nparam = isl_space_dim(dim, isl_dim_param);
27 v = Vector_Alloc(2 + nparam);
28 if (!v)
29 goto error;
31 for (i = 0; i < nparam; ++i)
32 value_set_si(v->p[2 + i], 0);
33 evalue_extract_affine(&e->x.p->arr[0], v->p + 2, &v->p[1], &v->p[0]);
35 ls = isl_local_space_from_space(isl_space_copy(dim));
36 aff = isl_aff_zero_on_domain(ls);
37 val = isl_val_int_from_gmp(ctx, v->p[1]);
38 aff = isl_aff_set_constant_val(aff, val);
40 for (i = 0; i < nparam; ++i) {
41 val = isl_val_int_from_gmp(ctx, v->p[2 + i]);
42 aff = isl_aff_set_coefficient_val(aff, isl_dim_param, i, val);
45 val = isl_val_int_from_gmp(ctx, v->p[0]);
46 aff = isl_aff_scale_down_val(aff, val);
48 aff = isl_aff_floor(aff);
49 base = isl_qpolynomial_from_aff(aff);
51 if (e->x.p->type == fractional) {
52 base = isl_qpolynomial_neg(base);
54 val = isl_val_from_gmp(ctx, v->p[1], v->p[0]);
55 c = isl_qpolynomial_val_on_domain(isl_space_copy(dim), val);
56 base = isl_qpolynomial_add(base, c);
58 for (i = 0; i < nparam; ++i) {
59 isl_qpolynomial *t;
60 val = isl_val_from_gmp(ctx, v->p[2 + i], v->p[0]);
61 c = isl_qpolynomial_val_on_domain(isl_space_copy(dim),
62 val);
63 t = isl_qpolynomial_var_on_domain(isl_space_copy(dim),
64 isl_dim_param, i);
65 t = isl_qpolynomial_mul(c, t);
66 base = isl_qpolynomial_add(base, t);
70 Vector_Free(v);
71 isl_space_free(dim);
73 return base;
74 error:
75 isl_space_free(dim);
76 return NULL;
79 static int type_offset(enode *p)
81 return p->type == fractional ? 1 :
82 p->type == flooring ? 1 : 0;
85 __isl_give isl_qpolynomial *isl_qpolynomial_from_evalue(__isl_take isl_space *dim,
86 const evalue *e)
88 int i;
89 isl_qpolynomial *qp;
90 isl_qpolynomial *base;
91 int offset;
93 if (EVALUE_IS_NAN(*e))
94 return isl_qpolynomial_infty_on_domain(dim);
95 if (value_notzero_p(e->d)) {
96 isl_ctx *ctx = isl_space_get_ctx(dim);
97 isl_val *val = isl_val_from_gmp(ctx, e->x.n, e->d);
98 return isl_qpolynomial_val_on_domain(dim, val);
101 offset = type_offset(e->x.p);
103 assert(e->x.p->type == polynomial ||
104 e->x.p->type == flooring ||
105 e->x.p->type == fractional);
106 assert(e->x.p->size >= 1 + offset);
108 base = extract_base(isl_space_copy(dim), e);
109 qp = isl_qpolynomial_from_evalue(isl_space_copy(dim),
110 &e->x.p->arr[e->x.p->size - 1]);
112 for (i = e->x.p->size - 2; i >= offset; --i) {
113 qp = isl_qpolynomial_mul(qp, isl_qpolynomial_copy(base));
114 qp = isl_qpolynomial_add(qp,
115 isl_qpolynomial_from_evalue(isl_space_copy(dim),
116 &e->x.p->arr[i]));
119 isl_qpolynomial_free(base);
120 isl_space_free(dim);
122 return qp;
125 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
126 const evalue *e);
128 static __isl_give isl_pw_qpolynomial *relation2pwqp(__isl_take isl_set *set,
129 const evalue *e)
131 int i;
132 Vector *vec;
133 isl_space *dim;
134 isl_ctx *ctx;
135 isl_val *v;
136 unsigned nparam;
137 isl_pw_qpolynomial *pwqp;
138 isl_aff *aff;
139 struct isl_constraint *c;
140 struct isl_basic_set *bset;
141 struct isl_set *guard;
142 const evalue *fract;
144 if (!set || !e)
145 goto error;
147 if (e->x.p->size == 1) {
148 dim = isl_set_get_space(set);
149 isl_set_free(set);
150 return isl_pw_qpolynomial_zero(dim);
153 ctx = isl_set_get_ctx(set);
154 isl_assert(ctx, e->x.p->size > 0, goto error);
155 isl_assert(ctx, e->x.p->size <= 3, goto error);
156 isl_assert(ctx, value_zero_p(e->x.p->arr[0].d), goto error);
157 isl_assert(ctx, e->x.p->arr[0].x.p->type == fractional, goto error);
158 fract = &e->x.p->arr[0];
160 nparam = isl_set_dim(set, isl_dim_param);
161 assert(isl_set_dim(set, isl_dim_set) == 0);
162 vec = Vector_Alloc(2 + nparam + 1);
163 if (!vec)
164 goto error;
166 for (i = 0; i < nparam; ++i)
167 value_set_si(vec->p[2 + i], 0);
168 evalue_extract_affine(&fract->x.p->arr[0],
169 vec->p + 2, &vec->p[1], &vec->p[0]);
171 dim = isl_set_get_space(set);
173 bset = isl_basic_set_universe(isl_space_copy(dim));
174 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
175 v = isl_val_int_from_gmp(ctx, vec->p[1]);
176 aff = isl_aff_set_constant_val(aff, v);
177 for (i = 0; i < nparam; ++i) {
178 v = isl_val_int_from_gmp(ctx, vec->p[2 + i]);
179 aff = isl_aff_set_coefficient_val(aff, isl_dim_param, i, v);
181 v = isl_val_int_from_gmp(ctx, vec->p[0]);
182 aff = isl_aff_mod_val(aff, v);
183 c = isl_equality_from_aff(aff);
184 bset = isl_basic_set_add_constraint(bset, c);
185 guard = isl_set_from_basic_set(bset);
186 Vector_Free(vec);
188 pwqp = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
189 isl_set_copy(guard)),
190 &e->x.p->arr[1]);
192 if (e->x.p->size == 3) {
193 isl_pw_qpolynomial *pwqpc;
194 guard = isl_set_complement(guard);
195 pwqpc = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
196 isl_set_copy(guard)),
197 &e->x.p->arr[2]);
198 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqpc);
201 isl_set_free(set);
202 isl_set_free(guard);
204 return pwqp;
205 error:
206 isl_set_free(set);
207 return NULL;
210 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
211 const evalue *e)
213 isl_qpolynomial *qp;
215 if (value_zero_p(e->d) && e->x.p->type == relation)
216 return relation2pwqp(set, e);
218 qp = isl_qpolynomial_from_evalue(isl_set_get_space(set), e);
220 return isl_pw_qpolynomial_alloc(set, qp);
223 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_evalue(__isl_take isl_space *dim, const evalue *e)
225 int i;
226 isl_space *pw_space;
227 isl_pw_qpolynomial *pwqp;
229 if (!dim || !e)
230 goto error;
231 if (EVALUE_IS_ZERO(*e)) {
232 dim = isl_space_from_domain(dim);
233 dim = isl_space_add_dims(dim, isl_dim_out, 1);
234 return isl_pw_qpolynomial_zero(dim);
237 if (value_notzero_p(e->d)) {
238 isl_ctx *ctx = isl_space_get_ctx(dim);
239 isl_set *set = isl_set_universe(isl_space_copy(dim));
240 isl_val *val = isl_val_from_gmp(ctx, e->x.n, e->d);
241 isl_qpolynomial *qp = isl_qpolynomial_val_on_domain(dim, val);
242 return isl_pw_qpolynomial_alloc(set, qp);
245 assert(!EVALUE_IS_NAN(*e));
247 assert(e->x.p->type == partition);
249 pw_space = isl_space_copy(dim);
250 pw_space = isl_space_from_domain(pw_space);
251 pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
252 pwqp = isl_pw_qpolynomial_zero(pw_space);
254 for (i = 0; i < e->x.p->size/2; ++i) {
255 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2 * i]);
256 isl_set *set = isl_set_new_from_polylib(D, isl_space_copy(dim));
257 isl_pw_qpolynomial *pwqp_i;
259 pwqp_i = guarded_evalue2pwqp(set, &e->x.p->arr[2 * i + 1]);
261 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqp_i);
264 isl_space_free(dim);
266 return pwqp;
267 error:
268 isl_space_free(dim);
269 return NULL;
272 static evalue *evalue_pow(evalue *e, int exp)
274 evalue *pow;
276 if (exp == 1)
277 return e;
279 pow = evalue_dup(e);
280 while (--exp > 0)
281 emul(e, pow);
283 evalue_free(e);
285 return pow;
288 static evalue *div2evalue(__isl_take isl_aff *div)
290 int i;
291 isl_ctx *ctx;
292 isl_val *v;
293 Vector *vec = NULL;
294 unsigned dim;
295 unsigned nparam;
296 evalue *e;
297 evalue *aff;
299 if (!div)
300 return NULL;
302 if (isl_aff_dim(div, isl_dim_div) != 0)
303 isl_die(isl_aff_get_ctx(div), isl_error_unsupported,
304 "cannot handle nested divs", goto error);
306 dim = isl_aff_dim(div, isl_dim_in);
307 nparam = isl_aff_dim(div, isl_dim_param);
309 ctx = isl_aff_get_ctx(div);
310 vec = Vector_Alloc(1 + dim + nparam + 1);
311 if (!vec)
312 goto error;
313 v = isl_aff_get_denominator_val(div);
314 isl_val_get_num_gmp(v, vec->p[0]);
315 div = isl_aff_scale_val(div, v);
316 for (i = 0; i < dim; ++i) {
317 v = isl_aff_get_coefficient_val(div, isl_dim_in, i);
318 isl_val_get_num_gmp(v, vec->p[1 + i]);
319 isl_val_free(v);
321 for (i = 0; i < nparam; ++i) {
322 v = isl_aff_get_coefficient_val(div, isl_dim_param, i);
323 isl_val_get_num_gmp(v, vec->p[1 + dim + i]);
324 isl_val_free(v);
326 v = isl_aff_get_constant_val(div);
327 isl_val_get_num_gmp(v, vec->p[1 + dim + nparam]);
328 isl_val_free(v);
330 e = isl_alloc_type(ctx, evalue);
331 if (!e)
332 goto error;
333 value_init(e->d);
334 value_set_si(e->d, 0);
335 e->x.p = new_enode(flooring, 3, -1);
336 evalue_set_si(&e->x.p->arr[1], 0, 1);
337 evalue_set_si(&e->x.p->arr[2], 1, 1);
338 value_clear(e->x.p->arr[0].d);
339 aff = affine2evalue(vec->p + 1, vec->p[0], dim + nparam);
340 e->x.p->arr[0] = *aff;
341 free(aff);
342 Vector_Free(vec);
343 isl_aff_free(div);
344 return e;
345 error:
346 Vector_Free(vec);
347 isl_aff_free(div);
348 return NULL;
351 static isl_stat add_term(__isl_take isl_term *term, void *user)
353 int i;
354 evalue *sum = (evalue *)user;
355 unsigned nparam;
356 unsigned dim;
357 unsigned n_div;
358 isl_ctx *ctx;
359 isl_val *v;
360 Value n, d;
361 evalue *e;
363 if (!term)
364 return isl_stat_error;
366 nparam = isl_term_dim(term, isl_dim_param);
367 dim = isl_term_dim(term, isl_dim_set);
368 n_div = isl_term_dim(term, isl_dim_div);
370 ctx = isl_term_get_ctx(term);
371 e = isl_alloc_type(ctx, evalue);
372 if (!e)
373 goto error;
375 value_init(n);
376 value_init(d);
378 v = isl_term_get_coefficient_val(term);
379 isl_val_get_num_gmp(v, n);
380 isl_val_get_den_gmp(v, d);
381 isl_val_free(v);
382 if (!v)
383 goto error2;
384 value_init(e->d);
385 evalue_set(e, n, d);
387 for (i = 0; i < dim; ++i) {
388 evalue *pow;
389 int exp = isl_term_get_exp(term, isl_dim_set, i);
391 if (!exp)
392 continue;
394 pow = evalue_pow(evalue_var(i), exp);
395 emul(pow, e);
396 evalue_free(pow);
399 for (i = 0; i < nparam; ++i) {
400 evalue *pow;
401 int exp = isl_term_get_exp(term, isl_dim_param, i);
403 if (!exp)
404 continue;
406 pow = evalue_pow(evalue_var(dim + i), exp);
407 emul(pow, e);
408 evalue_free(pow);
411 for (i = 0; i < n_div; ++i) {
412 evalue *pow;
413 evalue *floor;
414 isl_aff *div;
415 int exp = isl_term_get_exp(term, isl_dim_div, i);
417 if (!exp)
418 continue;
420 div = isl_term_get_div(term, i);
421 floor = div2evalue(div);
422 if (!floor)
423 goto error2;
424 pow = evalue_pow(floor, exp);
425 emul(pow, e);
426 evalue_free(pow);
429 eadd(e, sum);
430 evalue_free(e);
432 value_clear(n);
433 value_clear(d);
435 isl_term_free(term);
437 return isl_stat_ok;
438 error2:
439 evalue_free(e);
440 value_clear(n);
441 value_clear(d);
442 error:
443 isl_term_free(term);
444 return isl_stat_error;
447 evalue *isl_qpolynomial_to_evalue(__isl_keep isl_qpolynomial *qp)
449 evalue *e;
451 e = evalue_zero();
452 if (!e)
453 return NULL;
455 if (isl_qpolynomial_foreach_term(qp, add_term, e) < 0)
456 goto error;
458 return e;
459 error:
460 evalue_free(e);
461 return NULL;
464 static isl_stat add_guarded_qp(__isl_take isl_set *set,
465 __isl_take isl_qpolynomial *qp, void *user)
467 Polyhedron *D;
468 evalue *e = NULL;
469 evalue *f;
470 evalue *sum = (evalue *)user;
471 unsigned dim;
473 e = isl_alloc_type(isl_set_get_ctx(set), evalue);
474 if (!e)
475 goto error;
477 D = isl_set_to_polylib(set);
478 if (!D)
479 goto error;
481 f = isl_qpolynomial_to_evalue(qp);
482 if (!e) {
483 Domain_Free(D);
484 goto error;
487 dim = isl_set_dim(set, isl_dim_param) + isl_set_dim(set, isl_dim_set);
488 value_init(e->d);
489 e->x.p = new_enode(partition, 2, D->Dimension);
490 EVALUE_SET_DOMAIN(e->x.p->arr[0], D);
492 value_clear(e->x.p->arr[1].d);
493 e->x.p->arr[1] = *f;
494 free(f);
496 eadd(e, sum);
497 evalue_free(e);
499 isl_set_free(set);
500 isl_qpolynomial_free(qp);
502 return isl_stat_ok;
503 error:
504 free(e);
505 isl_set_free(set);
506 isl_qpolynomial_free(qp);
507 return isl_stat_error;
510 evalue *isl_pw_qpolynomial_to_evalue(__isl_keep isl_pw_qpolynomial *pwqp)
512 evalue *e;
514 if (!pwqp)
515 return NULL;
516 e = evalue_zero();
518 if (isl_pw_qpolynomial_foreach_piece(pwqp, add_guarded_qp, e))
519 goto error;
521 return e;
522 error:
523 evalue_free(e);
524 return NULL;