evalue_isl.c

   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>
   6
   7 static __isl_give isl_qpolynomial *extract_base(__isl_take isl_space *dim,
   8         const evalue *e)
   9 {
  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;
  18
  19         if (!dim)
  20                 return NULL;
  21
  22         if (e->x.p->type == polynomial)
  23                 return isl_qpolynomial_var_on_domain(dim, isl_dim_param, e->x.p->pos - 1);
  24
  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;
  30
  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]);
  34
  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);
  39
  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);
  43         }
  44
  45         val = isl_val_int_from_gmp(ctx, v->p[0]);
  46         aff = isl_aff_scale_down_val(aff, val);
  47
  48         aff = isl_aff_floor(aff);
  49         base = isl_qpolynomial_from_aff(aff);
  50
  51         if (e->x.p->type == fractional) {
  52                 base = isl_qpolynomial_neg(base);
  53
  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);
  57
  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);
  67                 }
  68         }
  69
  70         Vector_Free(v);
  71         isl_space_free(dim);
  72
  73         return base;
  74 error:
  75         isl_space_free(dim);
  76         return NULL;
  77 }
  78
  79 static int type_offset(enode *p)
  80 {
  81         return p->type == fractional ? 1 :
  82                p->type == flooring ? 1 : 0;
  83 }
  84
  85 __isl_give isl_qpolynomial *isl_qpolynomial_from_evalue(__isl_take isl_space *dim,
  86         const evalue *e)
  87 {
  88         int i;
  89         isl_qpolynomial *qp;
  90         isl_qpolynomial *base;
  91         int offset;
  92
  93         if (EVALUE_IS_NAN(*e))
  94                 return isl_qpolynomial_infty_on_domain(dim);
  95         if (value_notzero_p(e->d))
  96                 return isl_qpolynomial_rat_cst_on_domain(dim, e->x.n, e->d);
  97
  98         offset = type_offset(e->x.p);
  99
 100         assert(e->x.p->type == polynomial ||
 101                e->x.p->type == flooring ||
 102                e->x.p->type == fractional);
 103         assert(e->x.p->size >= 1 + offset);
 104
 105         base = extract_base(isl_space_copy(dim), e);
 106         qp = isl_qpolynomial_from_evalue(isl_space_copy(dim),
 107                                          &e->x.p->arr[e->x.p->size - 1]);
 108
 109         for (i = e->x.p->size - 2; i >= offset; --i) {
 110                 qp = isl_qpolynomial_mul(qp, isl_qpolynomial_copy(base));
 111                 qp = isl_qpolynomial_add(qp,
 112                                     isl_qpolynomial_from_evalue(isl_space_copy(dim),
 113                                     &e->x.p->arr[i]));
 114         }
 115
 116         isl_qpolynomial_free(base);
 117         isl_space_free(dim);
 118
 119         return qp;
 120 }
 121
 122 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
 123         const evalue *e);
 124
 125 static __isl_give isl_pw_qpolynomial *relation2pwqp(__isl_take isl_set *set,
 126         const evalue *e)
 127 {
 128         int i;
 129         Vector *vec;
 130         isl_space *dim;
 131         isl_ctx *ctx;
 132         isl_val *v;
 133         unsigned nparam;
 134         isl_pw_qpolynomial *pwqp;
 135         struct isl_constraint *c;
 136         struct isl_basic_set *bset;
 137         struct isl_set *guard;
 138         const evalue *fract;
 139
 140         if (!set || !e)
 141                 goto error;
 142
 143         if (e->x.p->size == 1) {
 144                 dim = isl_set_get_space(set);
 145                 isl_set_free(set);
 146                 return isl_pw_qpolynomial_zero(dim);
 147         }
 148
 149         ctx = isl_set_get_ctx(set);
 150         isl_assert(ctx, e->x.p->size > 0, goto error);
 151         isl_assert(ctx, e->x.p->size <= 3, goto error);
 152         isl_assert(ctx, value_zero_p(e->x.p->arr[0].d), goto error);
 153         isl_assert(ctx, e->x.p->arr[0].x.p->type == fractional, goto error);
 154         fract = &e->x.p->arr[0];
 155
 156         nparam = isl_set_dim(set, isl_dim_param);
 157         vec = Vector_Alloc(2 + nparam + 1);
 158         if (!vec)
 159                 goto error;
 160
 161         for (i = 0; i < nparam; ++i)
 162                 value_set_si(vec->p[2 + i], 0);
 163         evalue_extract_affine(&fract->x.p->arr[0],
 164                                 vec->p + 2, &vec->p[1], &vec->p[0]);
 165
 166         dim = isl_set_get_space(set);
 167         dim = isl_space_add_dims(dim, isl_dim_set, 1);
 168
 169         bset = isl_basic_set_universe(isl_space_copy(dim));
 170         c = isl_equality_alloc(isl_local_space_from_space(dim));
 171         v = isl_val_int_from_gmp(ctx, vec->p[0]);
 172         v = isl_val_neg(v);
 173         c = isl_constraint_set_coefficient_val(c, isl_dim_set, 0, v);
 174         v = isl_val_int_from_gmp(ctx, vec->p[1]);
 175         c = isl_constraint_set_constant_val(c, v);
 176         for (i = 0; i < nparam; ++i) {
 177                 v = isl_val_int_from_gmp(ctx, vec->p[2 + i]);
 178                 c = isl_constraint_set_coefficient_val(c, isl_dim_param, i, v);
 179         }
 180         bset = isl_basic_set_add_constraint(bset, c);
 181         bset = isl_basic_set_params(bset);
 182         guard = isl_set_from_basic_set(bset);
 183         Vector_Free(vec);
 184
 185         pwqp = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
 186                                                      isl_set_copy(guard)),
 187                                    &e->x.p->arr[1]);
 188
 189         if (e->x.p->size == 3) {
 190                 isl_pw_qpolynomial *pwqpc;
 191                 guard = isl_set_complement(guard);
 192                 pwqpc = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
 193                                                               isl_set_copy(guard)),
 194                                             &e->x.p->arr[2]);
 195                 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqpc);
 196         }
 197
 198         isl_set_free(set);
 199         isl_set_free(guard);
 200
 201         return pwqp;
 202 error:
 203         isl_set_free(set);
 204         return NULL;
 205 }
 206
 207 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
 208         const evalue *e)
 209 {
 210         isl_qpolynomial *qp;
 211
 212         if (value_zero_p(e->d) && e->x.p->type == relation)
 213                 return relation2pwqp(set, e);
 214
 215         qp = isl_qpolynomial_from_evalue(isl_set_get_space(set), e);
 216
 217         return isl_pw_qpolynomial_alloc(set, qp);
 218 }
 219
 220 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_evalue(__isl_take isl_space *dim, const evalue *e)
 221 {
 222         int i;
 223         isl_space *pw_space;
 224         isl_pw_qpolynomial *pwqp;
 225
 226         if (!dim || !e)
 227                 goto error;
 228         if (EVALUE_IS_ZERO(*e)) {
 229                 dim = isl_space_from_domain(dim);
 230                 dim = isl_space_add_dims(dim, isl_dim_out, 1);
 231                 return isl_pw_qpolynomial_zero(dim);
 232         }
 233
 234         if (value_notzero_p(e->d)) {
 235                 isl_set *set = isl_set_universe(isl_space_copy(dim));
 236                 isl_qpolynomial *qp = isl_qpolynomial_rat_cst_on_domain(dim, e->x.n, e->d);
 237                 return isl_pw_qpolynomial_alloc(set, qp);
 238         }
 239
 240         assert(!EVALUE_IS_NAN(*e));
 241
 242         assert(e->x.p->type == partition);
 243
 244         pw_space = isl_space_copy(dim);
 245         pw_space = isl_space_from_domain(pw_space);
 246         pw_space = isl_space_add_dims(pw_space, isl_dim_out, 1);
 247         pwqp = isl_pw_qpolynomial_zero(pw_space);
 248
 249         for (i = 0; i < e->x.p->size/2; ++i) {
 250                 Polyhedron *D = EVALUE_DOMAIN(e->x.p->arr[2 * i]);
 251                 isl_set *set = isl_set_new_from_polylib(D, isl_space_copy(dim));
 252                 isl_pw_qpolynomial *pwqp_i;
 253
 254                 pwqp_i = guarded_evalue2pwqp(set, &e->x.p->arr[2 * i + 1]);
 255
 256                 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqp_i);
 257         }
 258
 259         isl_space_free(dim);
 260
 261         return pwqp;
 262 error:
 263         isl_space_free(dim);
 264         return NULL;
 265 }
 266
 267 static evalue *evalue_pow(evalue *e, int exp)
 268 {
 269         evalue *pow;
 270
 271         if (exp == 1)
 272                 return e;
 273
 274         pow = evalue_dup(e);
 275         while (--exp > 0)
 276                 emul(e, pow);
 277
 278         evalue_free(e);
 279
 280         return pow;
 281 }
 282
 283 static evalue *div2evalue(__isl_take isl_aff *div)
 284 {
 285         int i;
 286         isl_ctx *ctx;
 287         Vector *vec = NULL;
 288         unsigned dim;
 289         unsigned nparam;
 290         evalue *e;
 291         evalue *aff;
 292
 293         if (!div)
 294                 return NULL;
 295
 296         if (isl_aff_dim(div, isl_dim_div) != 0)
 297                 isl_die(isl_aff_get_ctx(div), isl_error_unsupported,
 298                         "cannot handle nested divs", goto error);
 299
 300         dim = isl_aff_dim(div, isl_dim_in);
 301         nparam = isl_aff_dim(div, isl_dim_param);
 302
 303         ctx = isl_aff_get_ctx(div);
 304         vec = Vector_Alloc(1 + dim + nparam + 1);
 305         if (!vec)
 306                 goto error;
 307         for (i = 0; i < dim; ++i)
 308                 isl_aff_get_coefficient(div, isl_dim_in, i, &vec->p[1 + i]);
 309         for (i = 0; i < nparam; ++i)
 310                 isl_aff_get_coefficient(div, isl_dim_param, i,
 311                                         &vec->p[1 + dim + i]);
 312         isl_aff_get_denominator(div, &vec->p[0]);
 313         isl_aff_get_constant(div, &vec->p[1 + dim + nparam]);
 314
 315         e = isl_alloc_type(ctx, evalue);
 316         if (!e)
 317                 goto error;
 318         value_init(e->d);
 319         value_set_si(e->d, 0);
 320         e->x.p = new_enode(flooring, 3, -1);
 321         evalue_set_si(&e->x.p->arr[1], 0, 1);
 322         evalue_set_si(&e->x.p->arr[2], 1, 1);
 323         value_clear(e->x.p->arr[0].d);
 324         aff = affine2evalue(vec->p + 1, vec->p[0], dim + nparam);
 325         e->x.p->arr[0] = *aff;
 326         free(aff);
 327         Vector_Free(vec);
 328         isl_aff_free(div);
 329         return e;
 330 error:
 331         Vector_Free(vec);
 332         isl_aff_free(div);
 333         return NULL;
 334 }
 335
 336 static int add_term(__isl_take isl_term *term, void *user)
 337 {
 338         int i;
 339         evalue *sum = (evalue *)user;
 340         unsigned nparam;
 341         unsigned dim;
 342         unsigned n_div;
 343         isl_ctx *ctx;
 344         isl_val *v;
 345         Value n, d;
 346         evalue *e;
 347
 348         if (!term)
 349                 return -1;
 350
 351         nparam = isl_term_dim(term, isl_dim_param);
 352         dim = isl_term_dim(term, isl_dim_set);
 353         n_div = isl_term_dim(term, isl_dim_div);
 354
 355         ctx = isl_term_get_ctx(term);
 356         e = isl_alloc_type(ctx, evalue);
 357         if (!e)
 358                 goto error;
 359
 360         value_init(n);
 361         value_init(d);
 362
 363         v = isl_term_get_coefficient_val(term);
 364         isl_val_get_num_gmp(v, n);
 365         isl_val_get_den_gmp(v, d);
 366         isl_val_free(v);
 367         if (!v)
 368                 goto error2;
 369         value_init(e->d);
 370         evalue_set(e, n, d);
 371
 372         for (i = 0; i < dim; ++i) {
 373                 evalue *pow;
 374                 int exp = isl_term_get_exp(term, isl_dim_set, i);
 375
 376                 if (!exp)
 377                         continue;
 378
 379                 pow = evalue_pow(evalue_var(i), exp);
 380                 emul(pow, e);
 381                 evalue_free(pow);
 382         }
 383
 384         for (i = 0; i < nparam; ++i) {
 385                 evalue *pow;
 386                 int exp = isl_term_get_exp(term, isl_dim_param, i);
 387
 388                 if (!exp)
 389                         continue;
 390
 391                 pow = evalue_pow(evalue_var(dim + i), exp);
 392                 emul(pow, e);
 393                 evalue_free(pow);
 394         }
 395
 396         for (i = 0; i < n_div; ++i) {
 397                 evalue *pow;
 398                 evalue *floor;
 399                 isl_aff *div;
 400                 int exp = isl_term_get_exp(term, isl_dim_div, i);
 401
 402                 if (!exp)
 403                         continue;
 404
 405                 div = isl_term_get_div(term, i);
 406                 floor = div2evalue(div);
 407                 if (!floor)
 408                         goto error2;
 409                 pow = evalue_pow(floor, exp);
 410                 emul(pow, e);
 411                 evalue_free(pow);
 412         }
 413
 414         eadd(e, sum);
 415         evalue_free(e);
 416
 417         value_clear(n);
 418         value_clear(d);
 419
 420         isl_term_free(term);
 421
 422         return 0;
 423 error2:
 424         evalue_free(e);
 425         value_clear(n);
 426         value_clear(d);
 427 error:
 428         isl_term_free(term);
 429         return -1;
 430 }
 431
 432 evalue *isl_qpolynomial_to_evalue(__isl_keep isl_qpolynomial *qp)
 433 {
 434         evalue *e;
 435
 436         e = evalue_zero();
 437         if (!e)
 438                 return NULL;
 439
 440         if (isl_qpolynomial_foreach_term(qp, add_term, e) < 0)
 441                 goto error;
 442
 443         return e;
 444 error:
 445         evalue_free(e);
 446         return NULL;
 447 }
 448
 449 static int add_guarded_qp(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
 450         void *user)
 451 {
 452         Polyhedron *D;
 453         evalue *e = NULL;
 454         evalue *f;
 455         evalue *sum = (evalue *)user;
 456         unsigned dim;
 457
 458         e = isl_alloc_type(isl_set_get_ctx(set), evalue);
 459         if (!e)
 460                 goto error;
 461
 462         D = isl_set_to_polylib(set);
 463         if (!D)
 464                 goto error;
 465
 466         f = isl_qpolynomial_to_evalue(qp);
 467         if (!e) {
 468                 Domain_Free(D);
 469                 goto error;
 470         }
 471
 472         dim = isl_set_dim(set, isl_dim_param) + isl_set_dim(set, isl_dim_set);
 473         value_init(e->d);
 474         e->x.p = new_enode(partition, 2, D->Dimension);
 475         EVALUE_SET_DOMAIN(e->x.p->arr[0], D);
 476
 477         value_clear(e->x.p->arr[1].d);
 478         e->x.p->arr[1] = *f;
 479         free(f);
 480
 481         eadd(e, sum);
 482         evalue_free(e);
 483
 484         isl_set_free(set);
 485         isl_qpolynomial_free(qp);
 486
 487         return 0;
 488 error:
 489         free(e);
 490         isl_set_free(set);
 491         isl_qpolynomial_free(qp);
 492         return -1;
 493 }
 494
 495 evalue *isl_pw_qpolynomial_to_evalue(__isl_keep isl_pw_qpolynomial *pwqp)
 496 {
 497         evalue *e;
 498
 499         if (!pwqp)
 500                 return NULL;
 501         e = evalue_zero();
 502
 503         if (isl_pw_qpolynomial_foreach_piece(pwqp, add_guarded_qp, e))
 504                 goto error;
 505
 506         return e;
 507 error:
 508         evalue_free(e);
 509         return NULL;
 510 }