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