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