evalue_isl.c

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