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