edomain.cc

   1 #include <sstream>
   2 #include "fdstream.h"
   3 #include <barvinok/util.h>
   4 #include <barvinok/sample.h>
   5 #include <barvinok/barvinok.h>
   6 #include "edomain.h"
   7 #include "evalue_util.h"
   8
   9 using std::vector;
  10 using std::endl;
  11 using std::ostream;
  12
  13 EDomain *EDomain::new_from_ge_constraint(ge_constraint *ge, int sign,
  14                                          barvinok_options *options)
  15 {
  16     if (ge->simplified && sign == 0)
  17         return NULL;
  18
  19     EDomain *ED;
  20     Matrix *M2;
  21     M2 = Matrix_Copy(ge->M);
  22     if (sign == 1) {
  23         if (!ge->simplified)
  24             value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
  25                             M2->p[M2->NbRows-1][M2->NbColumns-1]);
  26     } else if (sign == -1) {
  27         Value mone;
  28         value_init(mone);
  29         value_set_si(mone, -1);
  30         Vector_Scale(M2->p[M2->NbRows-1]+1, M2->p[M2->NbRows-1]+1,
  31                      mone, M2->NbColumns-1);
  32         value_decrement(M2->p[M2->NbRows-1][M2->NbColumns-1],
  33                         M2->p[M2->NbRows-1][M2->NbColumns-1]);
  34         value_clear(mone);
  35     } else {
  36         value_set_si(M2->p[M2->NbRows-1][0], 0);
  37     }
  38     Polyhedron *D2 = Constraints2Polyhedron(M2, options->MaxRays);
  39     ED = new EDomain(D2, ge->D, ge->new_floors);
  40     Polyhedron_Free(D2);
  41     Matrix_Free(M2);
  42     return ED;
  43 }
  44
  45 static void print_term(ostream& os, Value v, int pos, int dim,
  46                         char **names, int *first)
  47 {
  48     if (value_zero_p(v)) {
  49         if (first && *first && pos >= dim)
  50             os << "0";
  51         return;
  52     }
  53
  54     if (first) {
  55         if (!*first && value_pos_p(v))
  56             os << "+";
  57         *first = 0;
  58     }
  59     if (pos < dim) {
  60         if (value_mone_p(v)) {
  61             os << "-";
  62         } else if (!value_one_p(v))
  63             os << VALUE_TO_INT(v);
  64         os << names[pos];
  65     } else
  66         os << VALUE_TO_INT(v);
  67 }
  68
  69 void EDomain_floor::print(ostream& os, char **p) const
  70 {
  71     int first = 1;
  72     os << "[";
  73     os << "(";
  74     for (int i = 0; i < v->Size-2; ++i)
  75         print_term(os, v->p[1+i], i, v->Size-2, p, &first);
  76     os << ")";
  77     os << "/";
  78     print_term(os, v->p[0], v->Size-2, v->Size-2, p, NULL);
  79     os << "]";
  80 }
  81
  82 void EDomain::print_constraints(ostream& os, char **p,
  83                                 barvinok_options *options) const
  84 {
  85     Value tmp;
  86     value_init(tmp);
  87
  88     Matrix *M = Matrix_Alloc(2*floors.size(), 1+D->Dimension+1);
  89     value_set_si(tmp, -1);
  90     for (int i = 0; i < floors.size(); ++i) {
  91         value_set_si(M->p[2*i][0], 1);
  92         Vector_Copy(floors[i]->v->p+1, M->p[2*i]+1, dimension());
  93         value_assign(M->p[2*i][1+D->Dimension], floors[i]->v->p[1+dimension()]);
  94         value_oppose(M->p[2*i][1+dimension()+i], floors[i]->v->p[0]);
  95
  96         Vector_Scale(M->p[2*i]+1, M->p[2*i+1]+1, tmp, D->Dimension+1);
  97         value_addto(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension],
  98                     M->p[2*i+1][1+dimension()+i]);
  99         value_decrement(M->p[2*i+1][1+D->Dimension], M->p[2*i+1][1+D->Dimension]);
 100         value_set_si(M->p[2*i+1][0], 1);
 101     }
 102     Polyhedron *E = Constraints2Polyhedron(M, options->MaxRays);
 103     Matrix_Free(M);
 104     Polyhedron *SD = DomainSimplify(D, E, options->MaxRays);
 105     Polyhedron_Free(E);
 106
 107     char **names = p;
 108     unsigned dim = dimension();
 109     if (dim < SD->Dimension) {
 110         names = new char * [SD->Dimension];
 111         int i;
 112         for (i = 0; i < dim; ++i)
 113             names[i] = p[i];
 114         for ( ; i < SD->Dimension; ++i) {
 115             std::ostringstream strm;
 116             floors[i-dim]->print(strm, p);
 117             names[i] = strdup(strm.str().c_str());
 118         }
 119     }
 120
 121     for (int i = 0; i < SD->NbConstraints; ++i) {
 122         int first = 1;
 123         int v = First_Non_Zero(SD->Constraint[i]+1, SD->Dimension);
 124         if (v == -1)
 125             continue;
 126         if (i)
 127             os << " && ";
 128         if (value_pos_p(SD->Constraint[i][v+1])) {
 129             print_term(os, SD->Constraint[i][v+1], v, SD->Dimension,
 130                        names, NULL);
 131             if (value_zero_p(SD->Constraint[i][0]))
 132                 os << " = ";
 133             else
 134                 os << " >= ";
 135             for (int j = v+1; j <= SD->Dimension; ++j) {
 136                 value_oppose(tmp, SD->Constraint[i][1+j]);
 137                 print_term(os, tmp, j, SD->Dimension,
 138                            names, &first);
 139             }
 140         } else {
 141             value_oppose(tmp, SD->Constraint[i][1+v]);
 142             print_term(os, tmp, v, SD->Dimension,
 143                        names, NULL);
 144             if (value_zero_p(SD->Constraint[i][0]))
 145                 os << " = ";
 146             else
 147                 os << " <= ";
 148             for (int j = v+1; j <= SD->Dimension; ++j)
 149                 print_term(os, SD->Constraint[i][1+j], j, SD->Dimension,
 150                            names, &first);
 151         }
 152     }
 153     value_clear(tmp);
 154     Domain_Free(SD);
 155
 156     if (dim < D->Dimension) {
 157         for (int i = dim; i < D->Dimension; ++i)
 158             free(names[i]);
 159         delete [] names;
 160     }
 161 }
 162
 163 void EDomain::print(FILE *out, char **p)
 164 {
 165     fdostream os(dup(fileno(out)));
 166     for (int i = 0; i < floors.size(); ++i) {
 167         os << "floor " << i << ": [";
 168         evalue_print(os, floors[i]->e, p);
 169         os << "]" << endl;
 170     }
 171     Polyhedron_Print(out, P_VALUE_FMT, D);
 172 }
 173
 174 static int type_offset(enode *p)
 175 {
 176    return p->type == fractional ? 1 :
 177           p->type == flooring ? 1 : 0;
 178 }
 179
 180 static void add_coeff(Value *cons, int len, evalue *coeff, int pos)
 181 {
 182     Value tmp;
 183
 184     assert(value_notzero_p(coeff->d));
 185
 186     value_init(tmp);
 187
 188     value_lcm(cons[0], coeff->d, &tmp);
 189     value_division(tmp, tmp, cons[0]);
 190     Vector_Scale(cons, cons, tmp, len);
 191     value_division(tmp, cons[0], coeff->d);
 192     value_addmul(cons[pos], tmp, coeff->x.n);
 193
 194     value_clear(tmp);
 195 }
 196
 197 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len);
 198
 199 static void add_fract(evalue *e, Value *cons, int len, int pos)
 200 {
 201     evalue mone;
 202     value_init(mone.d);
 203     evalue_set_si(&mone, -1, 1);
 204
 205     /* contribution of alpha * fract(X) is
 206      *      alpha * X ...
 207      */
 208     assert(e->x.p->size == 3);
 209     evalue argument;
 210     value_init(argument.d);
 211     evalue_copy(&argument, &e->x.p->arr[0]);
 212     emul(&e->x.p->arr[2], &argument);
 213     evalue2constraint_r(NULL, &argument, cons, len);
 214     free_evalue_refs(&argument);
 215
 216     /*      - alpha * floor(X) */
 217     emul(&mone, &e->x.p->arr[2]);
 218     add_coeff(cons, len, &e->x.p->arr[2], pos);
 219     emul(&mone, &e->x.p->arr[2]);
 220
 221     free_evalue_refs(&mone);
 222 }
 223
 224 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len)
 225 {
 226     int r = 0;
 227     if (value_zero_p(E->d) && E->x.p->type == fractional) {
 228         int i;
 229         assert(E->x.p->size == 3);
 230         r = evalue2constraint_r(D, &E->x.p->arr[1], cons, len);
 231         assert(value_notzero_p(E->x.p->arr[2].d));
 232         if (D && (i = D->find_floor(&E->x.p->arr[0])) >= 0) {
 233             add_fract(E, cons, len, 1+D->D->Dimension-D->floors.size()+i);
 234         } else {
 235             if (value_pos_p(E->x.p->arr[2].x.n)) {
 236                 evalue coeff;
 237                 value_init(coeff.d);
 238                 value_init(coeff.x.n);
 239                 value_set_si(coeff.d, 1);
 240                 evalue_denom(&E->x.p->arr[0], &coeff.d);
 241                 value_decrement(coeff.x.n, coeff.d);
 242                 emul(&E->x.p->arr[2], &coeff);
 243                 add_coeff(cons, len, &coeff, len-1);
 244                 free_evalue_refs(&coeff);
 245             }
 246             r = 1;
 247         }
 248     } else if (value_zero_p(E->d)) {
 249         assert(E->x.p->type == polynomial);
 250         assert(E->x.p->size == 2);
 251         r = evalue2constraint_r(D, &E->x.p->arr[0], cons, len) || r;
 252         add_coeff(cons, len, &E->x.p->arr[1], E->x.p->pos);
 253     } else {
 254         add_coeff(cons, len, E, len-1);
 255     }
 256     return r;
 257 }
 258
 259 EDomain_floor::EDomain_floor(const evalue *f, int dim)
 260 {
 261     e = new evalue;
 262     value_init(e->d);
 263     evalue_copy(e, f);
 264     v = Vector_Alloc(1+dim+1);
 265     value_set_si(v->p[0], 1);
 266     evalue2constraint_r(NULL, e, v->p, v->Size);
 267     refcount = 1;
 268     substituted = false;
 269 }
 270
 271 void EDomain_floor::eval(Value *values, Value *res) const
 272 {
 273     Inner_Product(v->p+1, values, v->Size-2, res);
 274     value_addto(*res, *res, v->p[v->Size-1]);
 275     value_pdivision(*res, *res, v->p[0]);
 276 }
 277
 278 int evalue2constraint(EDomain *D, evalue *E, Value *cons, int len)
 279 {
 280     Vector_Set(cons, 0, len);
 281     value_set_si(cons[0], 1);
 282     return evalue2constraint_r(D, E, cons, len);
 283 }
 284
 285 bool EDomain::contains(Value *point, int len) const
 286 {
 287     assert(len <= D->Dimension);
 288     if (len == D->Dimension)
 289         return in_domain(D, point);
 290
 291     Vector *all_val = Vector_Alloc(D->Dimension);
 292     Vector_Copy(point, all_val->p, len);
 293     for (int i = len-dimension(); i < floors.size(); ++i)
 294         floors[i]->eval(all_val->p, &all_val->p[dimension()+i]);
 295     bool in = in_domain(D, all_val->p);
 296     Vector_Free(all_val);
 297     return in;
 298 }
 299
 300 ge_constraint *EDomain::compute_ge_constraint(evalue *constraint) const
 301 {
 302     ge_constraint *ge = new ge_constraint(this);
 303     evalue mone;
 304     value_init(mone.d);
 305     evalue_set_si(&mone, -1, 1);
 306     int fract = 0;
 307     for (evalue *e = constraint; value_zero_p(e->d);
 308          e = &e->x.p->arr[type_offset(e->x.p)]) {
 309         int i;
 310         if (e->x.p->type != fractional)
 311             continue;
 312         if (find_floor(&e->x.p->arr[0]) == -1)
 313             ++fract;
 314     }
 315
 316     int rows = D->NbConstraints+2*fract+1;
 317     int cols = 2+D->Dimension+fract;
 318     ge->M = Matrix_Alloc(rows, cols);
 319     for (int i = 0; i < D->NbConstraints; ++i) {
 320         Vector_Copy(D->Constraint[i], ge->M->p[i], 1+D->Dimension);
 321         value_assign(ge->M->p[i][1+D->Dimension+fract],
 322                      D->Constraint[i][1+D->Dimension]);
 323     }
 324     value_set_si(ge->M->p[rows-1][0], 1);
 325     fract = 0;
 326     evalue *e;
 327     for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
 328         if (e->x.p->type == fractional) {
 329             int i, pos;
 330
 331             i = find_floor(&e->x.p->arr[0]);
 332             if (i >= 0)
 333                 pos = D->Dimension-floors.size()+i;
 334             else
 335                 pos = D->Dimension+fract;
 336
 337             add_fract(e, ge->M->p[rows-1], cols, 1+pos);
 338
 339             if (pos < D->Dimension)
 340                 continue;
 341
 342             EDomain_floor *new_floor;
 343             new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
 344
 345             /* constraints for the new floor */
 346             int row = D->NbConstraints+2*fract;
 347             Vector_Copy(new_floor->v->p+1, ge->M->p[row]+1, dimension());
 348             value_assign(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]);
 349             value_oppose(ge->M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
 350             value_set_si(ge->M->p[row][0], 1);
 351             assert(value_eq(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]));
 352             assert(Vector_Equal(new_floor->v->p+1, ge->M->p[row]+1, dimension()));
 353
 354             Vector_Scale(ge->M->p[row]+1, ge->M->p[row+1]+1, mone.x.n, cols-1);
 355             value_set_si(ge->M->p[row+1][0], 1);
 356             value_addto(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1],
 357                         ge->M->p[row+1][1+D->Dimension+fract]);
 358             value_decrement(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1]);
 359
 360             ge->new_floors.push_back(new_floor);
 361
 362             ++fract;
 363         } else {
 364             assert(e->x.p->type == polynomial);
 365             assert(e->x.p->size == 2);
 366             add_coeff(ge->M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
 367         }
 368     }
 369     add_coeff(ge->M->p[rows-1], cols, e, cols-1);
 370     ge->simplified = ConstraintSimplify(ge->M->p[rows-1], ge->M->p[rows-1],
 371                                          cols, &ge->M->p[rows-1][0]);
 372     value_set_si(ge->M->p[rows-1][0], 1);
 373     free_evalue_refs(&mone);
 374     return ge;
 375 }
 376
 377 /* "align" matrix to have nrows by inserting
 378  * the necessary number of rows and an equal number of columns at the end
 379  * right before the constant row/column
 380  */
 381 static Matrix *align_matrix_initial(Matrix *M, int nrows)
 382 {
 383     int i;
 384     int newrows = nrows - M->NbRows;
 385     Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
 386     for (i = 0; i < M->NbRows-1; ++i) {
 387         Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
 388         value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
 389     }
 390     for (i = 0; i <= newrows; ++i)
 391         value_assign(M2->p[M->NbRows-1+i][M->NbColumns-1+i],
 392                      M->p[M->NbRows-1][M->NbColumns-1]);
 393     return M2;
 394 }
 395
 396 static Matrix *InsertColumns(Matrix *M, int pos, int n)
 397 {
 398     Matrix *R;
 399     int i;
 400
 401     R = Matrix_Alloc(M->NbRows, M->NbColumns+n);
 402     for (i = 0; i < M->NbRows; ++i) {
 403         Vector_Copy(M->p[i], R->p[i], pos);
 404         Vector_Copy(M->p[i]+pos, R->p[i]+pos+n, M->NbColumns-pos);
 405     }
 406     return R;
 407 }
 408
 409 void evalue_substitute(evalue *e, evalue **subs)
 410 {
 411     evalue *v;
 412
 413     if (value_notzero_p(e->d))
 414         return;
 415
 416     enode *p = e->x.p;
 417     for (int i = 0; i < p->size; ++i)
 418         evalue_substitute(&p->arr[i], subs);
 419
 420     if (p->type == polynomial)
 421         v = subs[p->pos-1];
 422     else {
 423         v = new evalue;
 424         value_init(v->d);
 425         value_set_si(v->d, 0);
 426         v->x.p = new_enode(p->type, 3, -1);
 427         value_clear(v->x.p->arr[0].d);
 428         v->x.p->arr[0] = p->arr[0];
 429         evalue_set_si(&v->x.p->arr[1], 0, 1);
 430         evalue_set_si(&v->x.p->arr[2], 1, 1);
 431     }
 432
 433     int offset = type_offset(p);
 434
 435     for (int i = p->size-1; i >= offset+1; i--) {
 436         emul(v, &p->arr[i]);
 437         eadd(&p->arr[i], &p->arr[i-1]);
 438         free_evalue_refs(&(p->arr[i]));
 439     }
 440
 441     if (p->type != polynomial) {
 442         free_evalue_refs(v);
 443         delete v;
 444     }
 445
 446     value_clear(e->d);
 447     *e = p->arr[offset];
 448     free(p);
 449 }
 450
 451 void EDomain_floor::substitute(evalue **subs, Matrix *T)
 452 {
 453     /* This is a hack.  The EDomain_floor elements are possibly shared
 454      * by many EDomain structures and we want to perform the substitution
 455      * only once.
 456      */
 457     if (substituted)
 458         return;
 459     substituted = true;
 460
 461     evalue_substitute(e, subs);
 462
 463     assert(T->NbRows == v->Size-1);
 464     Vector *tmp = Vector_Alloc(1+T->NbColumns);
 465     Vector_Matrix_Product(v->p+1, T, tmp->p+1);
 466     value_multiply(tmp->p[0], v->p[0], T->p[T->NbRows-1][T->NbColumns-1]);
 467     Vector_Free(v);
 468     v = tmp;
 469 }
 470
 471 /* T is a homogeneous matrix that maps the original variables to the new variables.
 472  * this has constraints in the new variables and this method
 473  * transforms this to constraints in the original variables.
 474  */
 475 void EDomain::substitute(evalue **subs, Matrix *T, Matrix *Eq, unsigned MaxRays)
 476 {
 477     int nexist = floors.size();
 478     Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
 479     Polyhedron *new_D = DomainPreimage(D, M, MaxRays);
 480     Polyhedron_Free(D);
 481     D = new_D;
 482     Matrix_Free(M);
 483
 484     M = nexist ? InsertColumns(Eq, 1+dimension(), nexist) : Eq;
 485     new_D = DomainAddConstraints(D, M, MaxRays);
 486     Polyhedron_Free(D);
 487     D = new_D;
 488     if (nexist)
 489         Matrix_Free(M);
 490     for (int i = 0; i < floors.size(); ++i)
 491         floors[i]->substitute(subs, T);
 492 }
 493
 494 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
 495 {
 496     /* Matrix "view" of equalities */
 497     Matrix M;
 498     M.NbRows = (*P)->NbEq;
 499     M.NbColumns = (*P)->Dimension+2;
 500     M.p_Init = (*P)->p_Init;
 501     M.p = (*P)->Constraint;
 502
 503     Matrix *T = compress_variables(&M, nparam);
 504
 505     if (!T) {
 506         *P = NULL;
 507         return NULL;
 508     }
 509     if (isIdentity(T)) {
 510         Matrix_Free(T);
 511         T = NULL;
 512     } else
 513         *P = Polyhedron_Preimage(*P, T, MaxRays);
 514
 515     return T;
 516 }
 517
 518 bool EDomain::not_empty(lexmin_options *options)
 519 {
 520     Polyhedron *P = D;
 521     Polyhedron *Porig = P;
 522
 523     POL_ENSURE_VERTICES(P);
 524     if (emptyQ2(P))
 525         return false;
 526
 527     for (int i = 0; i < P->NbRays; ++i)
 528         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 529             sample = Vector_Alloc(P->Dimension + 1);
 530             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 531             return true;
 532         }
 533
 534     if (options->emptiness_check == BV_LEXMIN_EMPTINESS_CHECK_COUNT) {
 535         bool notzero;
 536         Value cb;
 537         value_init(cb);
 538         barvinok_count_with_options(P, &cb, options->barvinok);
 539         notzero = value_notzero_p(cb);
 540         value_clear(cb);
 541         return notzero;
 542     }
 543
 544     Matrix *T = NULL;
 545     while (P && !emptyQ2(P) && P->NbEq > 0) {
 546         Polyhedron *Q = P;
 547         Matrix *T2 = remove_equalities(&P, 0, options->barvinok->MaxRays);
 548         if (!T)
 549             T = T2;
 550         else {
 551             if (T2) {
 552                 Matrix *T3 = Matrix_Alloc(T->NbRows, T2->NbColumns);
 553                 Matrix_Product(T, T2, T3);
 554                 Matrix_Free(T);
 555                 Matrix_Free(T2);
 556                 T = T3;
 557             }
 558             if (Q != Porig)
 559                 Polyhedron_Free(Q);
 560         }
 561     }
 562     if (P)
 563         sample = Polyhedron_Sample(P, options->barvinok);
 564     if (sample) {
 565         if (T) {
 566             Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
 567             Matrix_Vector_Product(T, sample->p, P_sample->p);
 568             Vector_Free(sample);
 569             sample = P_sample;
 570         }
 571     }
 572     if (T) {
 573         Polyhedron_Free(P);
 574         Matrix_Free(T);
 575     }
 576
 577     if (sample)
 578         assert(in_domain(Porig, sample->p));
 579
 580     return sample != NULL;
 581 }