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 *evalue2isl(__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 }