evalue_isl.c

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