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     substituted = false;
 235 }
 236
 237 void EDomain_floor::eval(Value *values, Value *res) const
 238 {
 239     Inner_Product(v->p+1, values, v->Size-2, res);
 240     value_addto(*res, *res, v->p[v->Size-1]);
 241     value_pdivision(*res, *res, v->p[0]);
 242 }
 243
 244 int evalue2constraint(EDomain *D, evalue *E, Value *cons, int len)
 245 {
 246     Vector_Set(cons, 0, len);
 247     value_set_si(cons[0], 1);
 248     return evalue2constraint_r(D, E, cons, len);
 249 }
 250
 251 bool EDomain::contains(Value *point, int len) const
 252 {
 253     assert(len <= D->Dimension);
 254     if (len == D->Dimension)
 255         return in_domain(D, point);
 256
 257     Vector *all_val = Vector_Alloc(D->Dimension);
 258     Vector_Copy(point, all_val->p, len);
 259     for (int i = len-dimension(); i < floors.size(); ++i)
 260         floors[i]->eval(all_val->p, &all_val->p[dimension()+i]);
 261     bool in = in_domain(D, all_val->p);
 262     Vector_Free(all_val);
 263     return in;
 264 }
 265
 266 Matrix *EDomain::add_ge_constraint(evalue *constraint,
 267                                    vector<EDomain_floor *>& new_floors) const
 268 {
 269     evalue mone;
 270     value_init(mone.d);
 271     evalue_set_si(&mone, -1, 1);
 272     int fract = 0;
 273     for (evalue *e = constraint; value_zero_p(e->d);
 274          e = &e->x.p->arr[type_offset(e->x.p)]) {
 275         int i;
 276         if (e->x.p->type != fractional)
 277             continue;
 278         if (find_floor(&e->x.p->arr[0]) == -1)
 279             ++fract;
 280     }
 281
 282     int rows = D->NbConstraints+2*fract+1;
 283     int cols = 2+D->Dimension+fract;
 284     Matrix *M = Matrix_Alloc(rows, cols);
 285     for (int i = 0; i < D->NbConstraints; ++i) {
 286         Vector_Copy(D->Constraint[i], M->p[i], 1+D->Dimension);
 287         value_assign(M->p[i][1+D->Dimension+fract],
 288                      D->Constraint[i][1+D->Dimension]);
 289     }
 290     value_set_si(M->p[rows-1][0], 1);
 291     fract = 0;
 292     evalue *e;
 293     for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
 294         if (e->x.p->type == fractional) {
 295             int i, pos;
 296
 297             i = find_floor(&e->x.p->arr[0]);
 298             if (i >= 0)
 299                 pos = D->Dimension-floors.size()+i;
 300             else
 301                 pos = D->Dimension+fract;
 302
 303             add_fract(e, M->p[rows-1], cols, 1+pos);
 304
 305             if (pos < D->Dimension)
 306                 continue;
 307
 308             EDomain_floor *new_floor;
 309             new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
 310
 311             /* constraints for the new floor */
 312             int row = D->NbConstraints+2*fract;
 313             Vector_Copy(new_floor->v->p+1, M->p[row]+1, dimension());
 314             value_assign(M->p[row][cols-1], new_floor->v->p[1+dimension()]);
 315             value_oppose(M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
 316             value_set_si(M->p[row][0], 1);
 317             assert(value_eq(M->p[row][cols-1], new_floor->v->p[1+dimension()]));
 318             assert(Vector_Equal(new_floor->v->p+1, M->p[row]+1, dimension()));
 319
 320             Vector_Scale(M->p[row]+1, M->p[row+1]+1, mone.x.n, cols-1);
 321             value_set_si(M->p[row+1][0], 1);
 322             value_addto(M->p[row+1][cols-1], M->p[row+1][cols-1],
 323                         M->p[row+1][1+D->Dimension+fract]);
 324             value_decrement(M->p[row+1][cols-1], M->p[row+1][cols-1]);
 325
 326             new_floors.push_back(new_floor);
 327
 328             ++fract;
 329         } else {
 330             assert(e->x.p->type == polynomial);
 331             assert(e->x.p->size == 2);
 332             add_coeff(M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
 333         }
 334     }
 335     add_coeff(M->p[rows-1], cols, e, cols-1);
 336     value_set_si(M->p[rows-1][0], 1);
 337     free_evalue_refs(&mone);
 338     return M;
 339 }
 340
 341 /* "align" matrix to have nrows by inserting
 342  * the necessary number of rows and an equal number of columns at the end
 343  * right before the constant row/column
 344  */
 345 static Matrix *align_matrix_initial(Matrix *M, int nrows)
 346 {
 347     int i;
 348     int newrows = nrows - M->NbRows;
 349     Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
 350     for (i = 0; i < M->NbRows-1; ++i) {
 351         Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
 352         value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
 353     }
 354     for (i = 0; i <= newrows; ++i)
 355         value_assign(M2->p[M->NbRows-1+i][M->NbColumns-1+i],
 356                      M->p[M->NbRows-1][M->NbColumns-1]);
 357     return M2;
 358 }
 359
 360 static Matrix *InsertColumns(Matrix *M, int pos, int n)
 361 {
 362     Matrix *R;
 363     int i;
 364
 365     R = Matrix_Alloc(M->NbRows, M->NbColumns+n);
 366     for (i = 0; i < M->NbRows; ++i) {
 367         Vector_Copy(M->p[i], R->p[i], pos);
 368         Vector_Copy(M->p[i]+pos, R->p[i]+pos+n, M->NbColumns-pos);
 369     }
 370     return R;
 371 }
 372
 373 void evalue_substitute(evalue *e, evalue **subs)
 374 {
 375     evalue *v;
 376
 377     if (value_notzero_p(e->d))
 378         return;
 379
 380     enode *p = e->x.p;
 381     for (int i = 0; i < p->size; ++i)
 382         evalue_substitute(&p->arr[i], subs);
 383
 384     if (p->type == polynomial)
 385         v = subs[p->pos-1];
 386     else {
 387         v = new evalue;
 388         value_init(v->d);
 389         value_set_si(v->d, 0);
 390         v->x.p = new_enode(p->type, 3, -1);
 391         value_clear(v->x.p->arr[0].d);
 392         v->x.p->arr[0] = p->arr[0];
 393         evalue_set_si(&v->x.p->arr[1], 0, 1);
 394         evalue_set_si(&v->x.p->arr[2], 1, 1);
 395     }
 396
 397     int offset = type_offset(p);
 398
 399     for (int i = p->size-1; i >= offset+1; i--) {
 400         emul(v, &p->arr[i]);
 401         eadd(&p->arr[i], &p->arr[i-1]);
 402         free_evalue_refs(&(p->arr[i]));
 403     }
 404
 405     if (p->type != polynomial) {
 406         free_evalue_refs(v);
 407         delete v;
 408     }
 409
 410     value_clear(e->d);
 411     *e = p->arr[offset];
 412     free(p);
 413 }
 414
 415 void EDomain_floor::substitute(evalue **subs, Matrix *T)
 416 {
 417     /* This is a hack.  The EDomain_floor elements are possibly shared
 418      * by many EDomain structures and we want to perform the substitution
 419      * only once.
 420      */
 421     if (substituted)
 422         return;
 423     substituted = true;
 424
 425     evalue_substitute(e, subs);
 426
 427     assert(T->NbRows == v->Size-1);
 428     Vector *tmp = Vector_Alloc(1+T->NbColumns);
 429     Vector_Matrix_Product(v->p+1, T, tmp->p+1);
 430     value_multiply(tmp->p[0], v->p[0], T->p[T->NbRows-1][T->NbColumns-1]);
 431     Vector_Free(v);
 432     v = tmp;
 433 }
 434
 435 /* T is a homogeneous matrix that maps the original variables to the new variables.
 436  * this has constraints in the new variables and this method
 437  * transforms this to constraints in the original variables.
 438  */
 439 void EDomain::substitute(evalue **subs, Matrix *T, Matrix *Eq, unsigned MaxRays)
 440 {
 441     int nexist = floors.size();
 442     Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
 443     Polyhedron *new_D = DomainPreimage(D, M, MaxRays);
 444     Polyhedron_Free(D);
 445     D = new_D;
 446     Matrix_Free(M);
 447
 448     M = nexist ? InsertColumns(Eq, 1+dimension(), nexist) : Eq;
 449     new_D = DomainAddConstraints(D, M, MaxRays);
 450     Polyhedron_Free(D);
 451     D = new_D;
 452     if (nexist)
 453         Matrix_Free(M);
 454     for (int i = 0; i < floors.size(); ++i)
 455         floors[i]->substitute(subs, T);
 456 }