edomain.cc

   1 #include <sstream>
   2 #include "fdstream.h"
   3 #include <barvinok/util.h>
   4 #include "edomain.h"
   5 #include "evalue_util.h"
   6
   7 using std::vector;
   8 using std::endl;
   9 using std::ostream;
  10
  11 static void print_term(ostream& os, Value v, int pos, int dim,
  12                         char **names, int *first)
  13 {
  14     if (value_zero_p(v)) {
  15         if (first && *first && pos >= dim)
  16             os << "0";
  17         return;
  18     }
  19
  20     if (first) {
  21         if (!*first && value_pos_p(v))
  22             os << "+";
  23         *first = 0;
  24     }
  25     if (pos < dim) {
  26         if (value_mone_p(v)) {
  27             os << "-";
  28         } else if (!value_one_p(v))
  29             os << VALUE_TO_INT(v);
  30         os << names[pos];
  31     } else
  32         os << VALUE_TO_INT(v);
  33 }
  34
  35 void EDomain_floor::print(ostream& os, char **p) const
  36 {
  37     int first = 1;
  38     os << "[";
  39     os << "(";
  40     for (int i = 0; i < v->Size-2; ++i)
  41         print_term(os, v->p[1+i], i, v->Size-2, p, &first);
  42     os << ")";
  43     os << "/";
  44     print_term(os, v->p[0], v->Size-2, v->Size-2, p, NULL);
  45     os << "]";
  46 }
  47
  48 void EDomain::print_constraints(ostream& os, char **p,
  49                                 barvinok_options *options) const
  50 {
  51     Value tmp;
  52     value_init(tmp);
  53
  54     Matrix *M = Matrix_Alloc(2*floors.size(), 1+D->Dimension+1);
  55     value_set_si(tmp, -1);
  56     for (int i = 0; i < floors.size(); ++i) {
  57         value_set_si(M->p[2*i][0], 1);
  58         Vector_Copy(floors[i]->v->p+1, M->p[2*i]+1, dimension());
  59         value_assign(M->p[2*i][1+D->Dimension], floors[i]->v->p[1+dimension()]);
  60         value_oppose(M->p[2*i][1+dimension()+i], floors[i]->v->p[0]);
  61
  62         Vector_Scale(M->p[2*i]+1, M->p[2*i+1]+1, tmp, D->Dimension+1);
  63         value_addto(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension],
  64                     M->p[2*i+1][1+dimension()+i]);
  65         value_decrement(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension]);
  66         value_set_si(M->p[2*i+1][0], 1);
  67     }
  68     Polyhedron *E = Constraints2Polyhedron(M, options->MaxRays);
  69     Matrix_Free(M);
  70     Polyhedron *SD = DomainSimplify(D, E, options->MaxRays);
  71     Polyhedron_Free(E);
  72
  73     char **names = p;
  74     unsigned dim = dimension();
  75     if (dim < SD->Dimension) {
  76         names = new char * [SD->Dimension];
  77         int i;
  78         for (i = 0; i < dim; ++i)
  79             names[i] = p[i];
  80         for ( ; i < SD->Dimension; ++i) {
  81             std::ostringstream strm;
  82             floors[i-dim]->print(strm, p);
  83             names[i] = strdup(strm.str().c_str());
  84         }
  85     }
  86
  87     for (int i = 0; i < SD->NbConstraints; ++i) {
  88         int first = 1;
  89         int v = First_Non_Zero(SD->Constraint[i]+1, SD->Dimension);
  90         if (v == -1)
  91             continue;
  92         if (i)
  93             os << " && ";
  94         if (value_pos_p(SD->Constraint[i][v+1])) {
  95             print_term(os, SD->Constraint[i][v+1], v, SD->Dimension,
  96                        names, NULL);
  97             if (value_zero_p(SD->Constraint[i][0]))
  98                 os << " = ";
  99             else
 100                 os << " >= ";
 101             for (int j = v+1; j <= SD->Dimension; ++j) {
 102                 value_oppose(tmp, SD->Constraint[i][1+j]);
 103                 print_term(os, tmp, j, SD->Dimension,
 104                            names, &first);
 105             }
 106         } else {
 107             value_oppose(tmp, SD->Constraint[i][1+v]);
 108             print_term(os, tmp, v, SD->Dimension,
 109                        names, NULL);
 110             if (value_zero_p(SD->Constraint[i][0]))
 111                 os << " = ";
 112             else
 113                 os << " <= ";
 114             for (int j = v+1; j <= SD->Dimension; ++j)
 115                 print_term(os, SD->Constraint[i][1+j], j, SD->Dimension,
 116                            names, &first);
 117         }
 118     }
 119     value_clear(tmp);
 120     Domain_Free(SD);
 121
 122     if (dim < D->Dimension) {
 123         for (int i = dim; i < D->Dimension; ++i)
 124             free(names[i]);
 125         delete [] names;
 126     }
 127 }
 128
 129 void EDomain::print(FILE *out, char **p)
 130 {
 131     fdostream os(dup(fileno(out)));
 132     for (int i = 0; i < floors.size(); ++i) {
 133         os << "floor " << i << ": [";
 134         evalue_print(os, floors[i]->e, p);
 135         os << "]" << endl;
 136     }
 137     Polyhedron_Print(out, P_VALUE_FMT, D);
 138 }
 139
 140 static int type_offset(enode *p)
 141 {
 142    return p->type == fractional ? 1 :
 143           p->type == flooring ? 1 : 0;
 144 }
 145
 146 static void add_coeff(Value *cons, int len, evalue *coeff, int pos)
 147 {
 148     Value tmp;
 149
 150     assert(value_notzero_p(coeff->d));
 151
 152     value_init(tmp);
 153
 154     value_lcm(cons[0], coeff->d, &tmp);
 155     value_division(tmp, tmp, cons[0]);
 156     Vector_Scale(cons, cons, tmp, len);
 157     value_division(tmp, cons[0], coeff->d);
 158     value_addmul(cons[pos], tmp, coeff->x.n);
 159
 160     value_clear(tmp);
 161 }
 162
 163 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len);
 164
 165 static void add_fract(evalue *e, Value *cons, int len, int pos)
 166 {
 167     evalue mone;
 168     value_init(mone.d);
 169     evalue_set_si(&mone, -1, 1);
 170
 171     /* contribution of alpha * fract(X) is
 172      *      alpha * X ...
 173      */
 174     assert(e->x.p->size == 3);
 175     evalue argument;
 176     value_init(argument.d);
 177     evalue_copy(&argument, &e->x.p->arr[0]);
 178     emul(&e->x.p->arr[2], &argument);
 179     evalue2constraint_r(NULL, &argument, cons, len);
 180     free_evalue_refs(&argument);
 181
 182     /*      - alpha * floor(X) */
 183     emul(&mone, &e->x.p->arr[2]);
 184     add_coeff(cons, len, &e->x.p->arr[2], pos);
 185     emul(&mone, &e->x.p->arr[2]);
 186
 187     free_evalue_refs(&mone);
 188 }
 189
 190 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len)
 191 {
 192     int r = 0;
 193     if (value_zero_p(E->d) && E->x.p->type == fractional) {
 194         int i;
 195         assert(E->x.p->size == 3);
 196         r = evalue2constraint_r(D, &E->x.p->arr[1], cons, len);
 197         assert(value_notzero_p(E->x.p->arr[2].d));
 198         if (D && (i = D->find_floor(&E->x.p->arr[0])) >= 0) {
 199             add_fract(E, cons, len, 1+D->D->Dimension-D->floors.size()+i);
 200         } else {
 201             if (value_pos_p(E->x.p->arr[2].x.n)) {
 202                 evalue coeff;
 203                 value_init(coeff.d);
 204                 value_init(coeff.x.n);
 205                 value_set_si(coeff.d, 1);
 206                 evalue_denom(&E->x.p->arr[0], &coeff.d);
 207                 value_decrement(coeff.x.n, coeff.d);
 208                 emul(&E->x.p->arr[2], &coeff);
 209                 add_coeff(cons, len, &coeff, len-1);
 210                 free_evalue_refs(&coeff);
 211             }
 212             r = 1;
 213         }
 214     } else if (value_zero_p(E->d)) {
 215         assert(E->x.p->type == polynomial);
 216         assert(E->x.p->size == 2);
 217         r = evalue2constraint_r(D, &E->x.p->arr[0], cons, len) || r;
 218         add_coeff(cons, len, &E->x.p->arr[1], E->x.p->pos);
 219     } else {
 220         add_coeff(cons, len, E, len-1);
 221     }
 222     return r;
 223 }
 224
 225 EDomain_floor::EDomain_floor(const evalue *f, int dim)
 226 {
 227     e = new evalue;
 228     value_init(e->d);
 229     evalue_copy(e, f);
 230     v = Vector_Alloc(1+dim+1);
 231     value_set_si(v->p[0], 1);
 232     evalue2constraint_r(NULL, e, v->p, v->Size);
 233     refcount = 1;
 234 }
 235
 236 void EDomain_floor::eval(Value *values, Value *res) const
 237 {
 238     Inner_Product(v->p+1, values, v->Size-2, res);
 239     value_addto(*res, *res, v->p[v->Size-1]);
 240     value_pdivision(*res, *res, v->p[0]);
 241 }
 242
 243 int evalue2constraint(EDomain *D, evalue *E, Value *cons, int len)
 244 {
 245     Vector_Set(cons, 0, len);
 246     value_set_si(cons[0], 1);
 247     return evalue2constraint_r(D, E, cons, len);
 248 }
 249
 250 bool EDomain::contains(Value *point, int len) const
 251 {
 252     assert(len <= D->Dimension);
 253     if (len == D->Dimension)
 254         return in_domain(D, point);
 255
 256     Vector *all_val = Vector_Alloc(D->Dimension);
 257     Vector_Copy(point, all_val->p, len);
 258     for (int i = len-dimension(); i < floors.size(); ++i)
 259         floors[i]->eval(all_val->p, &all_val->p[dimension()+i]);
 260     bool in = in_domain(D, all_val->p);
 261     Vector_Free(all_val);
 262     return in;
 263 }
 264
 265 Matrix *EDomain::add_ge_constraint(evalue *constraint,
 266                                    vector<EDomain_floor *>& new_floors) const
 267 {
 268     evalue mone;
 269     value_init(mone.d);
 270     evalue_set_si(&mone, -1, 1);
 271     int fract = 0;
 272     for (evalue *e = constraint; value_zero_p(e->d);
 273          e = &e->x.p->arr[type_offset(e->x.p)]) {
 274         int i;
 275         if (e->x.p->type != fractional)
 276             continue;
 277         if (find_floor(&e->x.p->arr[0]) == -1)
 278             ++fract;
 279     }
 280
 281     int rows = D->NbConstraints+2*fract+1;
 282     int cols = 2+D->Dimension+fract;
 283     Matrix *M = Matrix_Alloc(rows, cols);
 284     for (int i = 0; i < D->NbConstraints; ++i) {
 285         Vector_Copy(D->Constraint[i], M->p[i], 1+D->Dimension);
 286         value_assign(M->p[i][1+D->Dimension+fract],
 287                      D->Constraint[i][1+D->Dimension]);
 288     }
 289     value_set_si(M->p[rows-1][0], 1);
 290     fract = 0;
 291     evalue *e;
 292     for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
 293         if (e->x.p->type == fractional) {
 294             int i, pos;
 295
 296             i = find_floor(&e->x.p->arr[0]);
 297             if (i >= 0)
 298                 pos = D->Dimension-floors.size()+i;
 299             else
 300                 pos = D->Dimension+fract;
 301
 302             add_fract(e, M->p[rows-1], cols, 1+pos);
 303
 304             if (pos < D->Dimension)
 305                 continue;
 306
 307             EDomain_floor *new_floor;
 308             new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
 309
 310             /* constraints for the new floor */
 311             int row = D->NbConstraints+2*fract;
 312             Vector_Copy(new_floor->v->p+1, M->p[row]+1, dimension());
 313             value_assign(M->p[row][cols-1], new_floor->v->p[1+dimension()]);
 314             value_oppose(M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
 315             value_set_si(M->p[row][0], 1);
 316             assert(value_eq(M->p[row][cols-1], new_floor->v->p[1+dimension()]));
 317             assert(Vector_Equal(new_floor->v->p+1, M->p[row]+1, dimension()));
 318
 319             Vector_Scale(M->p[row]+1, M->p[row+1]+1, mone.x.n, cols-1);
 320             value_set_si(M->p[row+1][0], 1);
 321             value_addto(M->p[row+1][cols-1], M->p[row+1][cols-1],
 322                         M->p[row+1][1+D->Dimension+fract]);
 323             value_decrement(M->p[row+1][cols-1], M->p[row+1][cols-1]);
 324
 325             new_floors.push_back(new_floor);
 326
 327             ++fract;
 328         } else {
 329             assert(e->x.p->type == polynomial);
 330             assert(e->x.p->size == 2);
 331             add_coeff(M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
 332         }
 333     }
 334     add_coeff(M->p[rows-1], cols, e, cols-1);
 335     value_set_si(M->p[rows-1][0], 1);
 336     free_evalue_refs(&mone);
 337     return M;
 338 }