evalue_isl.c

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