evalue_isl.c

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