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