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                 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);
  99         }
 100
 101         offset = type_offset(e->x.p);
 102
 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);
 107
 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]);
 111
 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]));
 117         }
 118
 119         isl_qpolynomial_free(base);
 120         isl_space_free(dim);
 121
 122         return qp;
 123 }
 124
 125 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
 126         const evalue *e);
 127
 128 static __isl_give isl_pw_qpolynomial *relation2pwqp(__isl_take isl_set *set,
 129         const evalue *e)
 130 {
 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;
 143
 144         if (!set || !e)
 145                 goto error;
 146
 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);
 151         }
 152
 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];
 159
 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;
 165
 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]);
 170
 171         dim = isl_set_get_space(set);
 172
 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);
 180         }
 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);
 187
 188         pwqp = guarded_evalue2pwqp(isl_set_intersect(isl_set_copy(set),
 189                                                      isl_set_copy(guard)),
 190                                    &e->x.p->arr[1]);
 191
 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);
 199         }
 200
 201         isl_set_free(set);
 202         isl_set_free(guard);
 203
 204         return pwqp;
 205 error:
 206         isl_set_free(set);
 207         return NULL;
 208 }
 209
 210 static __isl_give isl_pw_qpolynomial *guarded_evalue2pwqp(__isl_take isl_set *set,
 211         const evalue *e)
 212 {
 213         isl_qpolynomial *qp;
 214
 215         if (value_zero_p(e->d) && e->x.p->type == relation)
 216                 return relation2pwqp(set, e);
 217
 218         qp = isl_qpolynomial_from_evalue(isl_set_get_space(set), e);
 219
 220         return isl_pw_qpolynomial_alloc(set, qp);
 221 }
 222
 223 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_evalue(__isl_take isl_space *dim, const evalue *e)
 224 {
 225         int i;
 226         isl_space *pw_space;
 227         isl_pw_qpolynomial *pwqp;
 228
 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);
 235         }
 236
 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);
 243         }
 244
 245         assert(!EVALUE_IS_NAN(*e));
 246
 247         assert(e->x.p->type == partition);
 248
 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);
 253
 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;
 258
 259                 pwqp_i = guarded_evalue2pwqp(set, &e->x.p->arr[2 * i + 1]);
 260
 261                 pwqp = isl_pw_qpolynomial_add_disjoint(pwqp, pwqp_i);
 262         }
 263
 264         isl_space_free(dim);
 265
 266         return pwqp;
 267 error:
 268         isl_space_free(dim);
 269         return NULL;
 270 }
 271
 272 static evalue *evalue_pow(evalue *e, int exp)
 273 {
 274         evalue *pow;
 275
 276         if (exp == 1)
 277                 return e;
 278
 279         pow = evalue_dup(e);
 280         while (--exp > 0)
 281                 emul(e, pow);
 282
 283         evalue_free(e);
 284
 285         return pow;
 286 }
 287
 288 static evalue *div2evalue(__isl_take isl_aff *div)
 289 {
 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;
 298
 299         if (!div)
 300                 return NULL;
 301
 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);
 305
 306         dim = isl_aff_dim(div, isl_dim_in);
 307         nparam = isl_aff_dim(div, isl_dim_param);
 308
 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);
 320         }
 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);
 325         }
 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);
 329
 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;
 349 }
 350
 351 static int add_term(__isl_take isl_term *term, void *user)
 352 {
 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;
 362
 363         if (!term)
 364                 return -1;
 365
 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);
 369
 370         ctx = isl_term_get_ctx(term);
 371         e = isl_alloc_type(ctx, evalue);
 372         if (!e)
 373                 goto error;
 374
 375         value_init(n);
 376         value_init(d);
 377
 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);
 386
 387         for (i = 0; i < dim; ++i) {
 388                 evalue *pow;
 389                 int exp = isl_term_get_exp(term, isl_dim_set, i);
 390
 391                 if (!exp)
 392                         continue;
 393
 394                 pow = evalue_pow(evalue_var(i), exp);
 395                 emul(pow, e);
 396                 evalue_free(pow);
 397         }
 398
 399         for (i = 0; i < nparam; ++i) {
 400                 evalue *pow;
 401                 int exp = isl_term_get_exp(term, isl_dim_param, i);
 402
 403                 if (!exp)
 404                         continue;
 405
 406                 pow = evalue_pow(evalue_var(dim + i), exp);
 407                 emul(pow, e);
 408                 evalue_free(pow);
 409         }
 410
 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);
 416
 417                 if (!exp)
 418                         continue;
 419
 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);
 427         }
 428
 429         eadd(e, sum);
 430         evalue_free(e);
 431
 432         value_clear(n);
 433         value_clear(d);
 434
 435         isl_term_free(term);
 436
 437         return 0;
 438 error2:
 439         evalue_free(e);
 440         value_clear(n);
 441         value_clear(d);
 442 error:
 443         isl_term_free(term);
 444         return -1;
 445 }
 446
 447 evalue *isl_qpolynomial_to_evalue(__isl_keep isl_qpolynomial *qp)
 448 {
 449         evalue *e;
 450
 451         e = evalue_zero();
 452         if (!e)
 453                 return NULL;
 454
 455         if (isl_qpolynomial_foreach_term(qp, add_term, e) < 0)
 456                 goto error;
 457
 458         return e;
 459 error:
 460         evalue_free(e);
 461         return NULL;
 462 }
 463
 464 static int add_guarded_qp(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
 465         void *user)
 466 {
 467         Polyhedron *D;
 468         evalue *e = NULL;
 469         evalue *f;
 470         evalue *sum = (evalue *)user;
 471         unsigned dim;
 472
 473         e = isl_alloc_type(isl_set_get_ctx(set), evalue);
 474         if (!e)
 475                 goto error;
 476
 477         D = isl_set_to_polylib(set);
 478         if (!D)
 479                 goto error;
 480
 481         f = isl_qpolynomial_to_evalue(qp);
 482         if (!e) {
 483                 Domain_Free(D);
 484                 goto error;
 485         }
 486
 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);
 491
 492         value_clear(e->x.p->arr[1].d);
 493         e->x.p->arr[1] = *f;
 494         free(f);
 495
 496         eadd(e, sum);
 497         evalue_free(e);
 498
 499         isl_set_free(set);
 500         isl_qpolynomial_free(qp);
 501
 502         return 0;
 503 error:
 504         free(e);
 505         isl_set_free(set);
 506         isl_qpolynomial_free(qp);
 507         return -1;
 508 }
 509
 510 evalue *isl_pw_qpolynomial_to_evalue(__isl_keep isl_pw_qpolynomial *pwqp)
 511 {
 512         evalue *e;
 513
 514         if (!pwqp)
 515                 return NULL;
 516         e = evalue_zero();
 517
 518         if (isl_pw_qpolynomial_foreach_piece(pwqp, add_guarded_qp, e))
 519                 goto error;
 520
 521         return e;
 522 error:
 523         evalue_free(e);
 524         return NULL;
 525 }