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 /*
 164 void EDomain::print(FILE *out, char **p)
 165 {
 166     fdostream os(dup(fileno(out)));
 167     for (int i = 0; i < floors.size(); ++i) {
 168         os << "floor " << i << ": [";
 169         evalue_print(os, floors[i]->e, p);
 170         os << "]" << endl;
 171     }
 172     Polyhedron_Print(out, P_VALUE_FMT, D);
 173 }
 174 */
 175
 176 static int type_offset(enode *p)
 177 {
 178    return p->type == fractional ? 1 :
 179           p->type == flooring ? 1 : 0;
 180 }
 181
 182 static void add_coeff(Value *cons, int len, evalue *coeff, int pos)
 183 {
 184     Value tmp;
 185
 186     assert(value_notzero_p(coeff->d));
 187
 188     value_init(tmp);
 189
 190     value_lcm(cons[0], coeff->d, &tmp);
 191     value_division(tmp, tmp, cons[0]);
 192     Vector_Scale(cons, cons, tmp, len);
 193     value_division(tmp, cons[0], coeff->d);
 194     value_addmul(cons[pos], tmp, coeff->x.n);
 195
 196     value_clear(tmp);
 197 }
 198
 199 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len);
 200
 201 static void add_fract(evalue *e, Value *cons, int len, int pos)
 202 {
 203     evalue mone;
 204     value_init(mone.d);
 205     evalue_set_si(&mone, -1, 1);
 206
 207     /* contribution of alpha * fract(X) is
 208      *      alpha * X ...
 209      */
 210     assert(e->x.p->size == 3);
 211     evalue argument;
 212     value_init(argument.d);
 213     evalue_copy(&argument, &e->x.p->arr[0]);
 214     emul(&e->x.p->arr[2], &argument);
 215     evalue2constraint_r(NULL, &argument, cons, len);
 216     free_evalue_refs(&argument);
 217
 218     /*      - alpha * floor(X) */
 219     emul(&mone, &e->x.p->arr[2]);
 220     add_coeff(cons, len, &e->x.p->arr[2], pos);
 221     emul(&mone, &e->x.p->arr[2]);
 222
 223     free_evalue_refs(&mone);
 224 }
 225
 226 static int evalue2constraint_r(EDomain *D, evalue *E, Value *cons, int len)
 227 {
 228     int r = 0;
 229     if (value_zero_p(E->d) && E->x.p->type == fractional) {
 230         int i;
 231         assert(E->x.p->size == 3);
 232         r = evalue2constraint_r(D, &E->x.p->arr[1], cons, len);
 233         assert(value_notzero_p(E->x.p->arr[2].d));
 234         if (D && (i = D->find_floor(&E->x.p->arr[0])) >= 0) {
 235             add_fract(E, cons, len, 1+D->D->Dimension-D->floors.size()+i);
 236         } else {
 237             if (value_pos_p(E->x.p->arr[2].x.n)) {
 238                 evalue coeff;
 239                 value_init(coeff.d);
 240                 value_init(coeff.x.n);
 241                 value_set_si(coeff.d, 1);
 242                 evalue_denom(&E->x.p->arr[0], &coeff.d);
 243                 value_decrement(coeff.x.n, coeff.d);
 244                 emul(&E->x.p->arr[2], &coeff);
 245                 add_coeff(cons, len, &coeff, len-1);
 246                 free_evalue_refs(&coeff);
 247             }
 248             r = 1;
 249         }
 250     } else if (value_zero_p(E->d)) {
 251         assert(E->x.p->type == polynomial);
 252         assert(E->x.p->size == 2);
 253         r = evalue2constraint_r(D, &E->x.p->arr[0], cons, len) || r;
 254         add_coeff(cons, len, &E->x.p->arr[1], E->x.p->pos);
 255     } else {
 256         add_coeff(cons, len, E, len-1);
 257     }
 258     return r;
 259 }
 260
 261 EDomain_floor::EDomain_floor(const evalue *f, int dim)
 262 {
 263     e = new evalue;
 264     value_init(e->d);
 265     evalue_copy(e, f);
 266     v = Vector_Alloc(1+dim+1);
 267     value_set_si(v->p[0], 1);
 268     evalue2constraint_r(NULL, e, v->p, v->Size);
 269     refcount = 1;
 270     substituted = false;
 271 }
 272
 273 void EDomain_floor::eval(Value *values, Value *res) const
 274 {
 275     Inner_Product(v->p+1, values, v->Size-2, res);
 276     value_addto(*res, *res, v->p[v->Size-1]);
 277     value_pdivision(*res, *res, v->p[0]);
 278 }
 279
 280 int evalue2constraint(EDomain *D, evalue *E, Value *cons, int len)
 281 {
 282     Vector_Set(cons, 0, len);
 283     value_set_si(cons[0], 1);
 284     return evalue2constraint_r(D, E, cons, len);
 285 }
 286
 287 bool EDomain::contains(Value *point, int len) const
 288 {
 289     assert(len <= D->Dimension);
 290     if (len == D->Dimension)
 291         return in_domain(D, point);
 292
 293     Vector *all_val = Vector_Alloc(D->Dimension);
 294     Vector_Copy(point, all_val->p, len);
 295     for (int i = len-dimension(); i < floors.size(); ++i)
 296         floors[i]->eval(all_val->p, &all_val->p[dimension()+i]);
 297     bool in = in_domain(D, all_val->p);
 298     Vector_Free(all_val);
 299     return in;
 300 }
 301
 302 ge_constraint *EDomain::compute_ge_constraint(evalue *constraint) const
 303 {
 304     ge_constraint *ge = new ge_constraint(this);
 305     evalue mone;
 306     value_init(mone.d);
 307     evalue_set_si(&mone, -1, 1);
 308     int fract = 0;
 309     for (evalue *e = constraint; value_zero_p(e->d);
 310          e = &e->x.p->arr[type_offset(e->x.p)]) {
 311         int i;
 312         if (e->x.p->type != fractional)
 313             continue;
 314         if (find_floor(&e->x.p->arr[0]) == -1)
 315             ++fract;
 316     }
 317
 318     int rows = D->NbConstraints+2*fract+1;
 319     int cols = 2+D->Dimension+fract;
 320     ge->M = Matrix_Alloc(rows, cols);
 321     for (int i = 0; i < D->NbConstraints; ++i) {
 322         Vector_Copy(D->Constraint[i], ge->M->p[i], 1+D->Dimension);
 323         value_assign(ge->M->p[i][1+D->Dimension+fract],
 324                      D->Constraint[i][1+D->Dimension]);
 325     }
 326     value_set_si(ge->M->p[rows-1][0], 1);
 327     fract = 0;
 328     evalue *e;
 329     for (e = constraint; value_zero_p(e->d); e = &e->x.p->arr[type_offset(e->x.p)]) {
 330         if (e->x.p->type == fractional) {
 331             int i, pos;
 332
 333             i = find_floor(&e->x.p->arr[0]);
 334             if (i >= 0)
 335                 pos = D->Dimension-floors.size()+i;
 336             else
 337                 pos = D->Dimension+fract;
 338
 339             add_fract(e, ge->M->p[rows-1], cols, 1+pos);
 340
 341             if (pos < D->Dimension)
 342                 continue;
 343
 344             EDomain_floor *new_floor;
 345             new_floor = new EDomain_floor(&e->x.p->arr[0], dimension());
 346
 347             /* constraints for the new floor */
 348             int row = D->NbConstraints+2*fract;
 349             Vector_Copy(new_floor->v->p+1, ge->M->p[row]+1, dimension());
 350             value_assign(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]);
 351             value_oppose(ge->M->p[row][1+D->Dimension+fract], new_floor->v->p[0]);
 352             value_set_si(ge->M->p[row][0], 1);
 353             assert(value_eq(ge->M->p[row][cols-1], new_floor->v->p[1+dimension()]));
 354             assert(Vector_Equal(new_floor->v->p+1, ge->M->p[row]+1, dimension()));
 355
 356             Vector_Scale(ge->M->p[row]+1, ge->M->p[row+1]+1, mone.x.n, cols-1);
 357             value_set_si(ge->M->p[row+1][0], 1);
 358             value_addto(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1],
 359                         ge->M->p[row+1][1+D->Dimension+fract]);
 360             value_decrement(ge->M->p[row+1][cols-1], ge->M->p[row+1][cols-1]);
 361
 362             ge->new_floors.push_back(new_floor);
 363
 364             ++fract;
 365         } else {
 366             assert(e->x.p->type == polynomial);
 367             assert(e->x.p->size == 2);
 368             add_coeff(ge->M->p[rows-1], cols, &e->x.p->arr[1], e->x.p->pos);
 369         }
 370     }
 371     add_coeff(ge->M->p[rows-1], cols, e, cols-1);
 372     ge->simplified = ConstraintSimplify(ge->M->p[rows-1], ge->M->p[rows-1],
 373                                          cols, &ge->M->p[rows-1][0]);
 374     value_set_si(ge->M->p[rows-1][0], 1);
 375     free_evalue_refs(&mone);
 376     return ge;
 377 }
 378
 379 /* "align" matrix to have nrows by inserting
 380  * the necessary number of rows and an equal number of columns at the end
 381  * right before the constant row/column
 382  */
 383 static Matrix *align_matrix_initial(Matrix *M, int nrows)
 384 {
 385     int i;
 386     int newrows = nrows - M->NbRows;
 387     Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
 388     for (i = 0; i < M->NbRows-1; ++i) {
 389         Vector_Copy(M->p[i], M2->p[i], M->NbColumns-1);
 390         value_assign(M2->p[i][M2->NbColumns-1], M->p[i][M->NbColumns-1]);
 391     }
 392     for (i = 0; i <= newrows; ++i)
 393         value_assign(M2->p[M->NbRows-1+i][M->NbColumns-1+i],
 394                      M->p[M->NbRows-1][M->NbColumns-1]);
 395     return M2;
 396 }
 397
 398 static Matrix *InsertColumns(Matrix *M, int pos, int n)
 399 {
 400     Matrix *R;
 401     int i;
 402
 403     R = Matrix_Alloc(M->NbRows, M->NbColumns+n);
 404     for (i = 0; i < M->NbRows; ++i) {
 405         Vector_Copy(M->p[i], R->p[i], pos);
 406         Vector_Copy(M->p[i]+pos, R->p[i]+pos+n, M->NbColumns-pos);
 407     }
 408     return R;
 409 }
 410
 411 void EDomain_floor::substitute(evalue **subs, Matrix *T)
 412 {
 413     /* This is a hack.  The EDomain_floor elements are possibly shared
 414      * by many EDomain structures and we want to perform the substitution
 415      * only once.
 416      */
 417     if (substituted)
 418         return;
 419     substituted = true;
 420
 421     evalue_substitute(e, subs);
 422
 423     assert(T->NbRows == v->Size-1);
 424     Vector *tmp = Vector_Alloc(1+T->NbColumns);
 425     Vector_Matrix_Product(v->p+1, T, tmp->p+1);
 426     value_multiply(tmp->p[0], v->p[0], T->p[T->NbRows-1][T->NbColumns-1]);
 427     Vector_Free(v);
 428     v = tmp;
 429 }
 430
 431 /* T is a homogeneous matrix that maps the original variables to the new variables.
 432  * this has constraints in the new variables and this method
 433  * transforms this to constraints in the original variables.
 434  */
 435 void EDomain::substitute(evalue **subs, Matrix *T, Matrix *Eq, unsigned MaxRays)
 436 {
 437     int nexist = floors.size();
 438     Matrix *M = align_matrix_initial(T, T->NbRows+nexist);
 439     Polyhedron *new_D = DomainPreimage(D, M, MaxRays);
 440     Polyhedron_Free(D);
 441     D = new_D;
 442     Matrix_Free(M);
 443
 444     M = nexist ? InsertColumns(Eq, 1+dimension(), nexist) : Eq;
 445     new_D = DomainAddConstraints(D, M, MaxRays);
 446     Polyhedron_Free(D);
 447     D = new_D;
 448     if (nexist)
 449         Matrix_Free(M);
 450     for (int i = 0; i < floors.size(); ++i)
 451         floors[i]->substitute(subs, T);
 452 }
 453
 454 static Matrix *remove_equalities(Polyhedron **P, unsigned nparam, unsigned MaxRays)
 455 {
 456     /* Matrix "view" of equalities */
 457     Matrix M;
 458     M.NbRows = (*P)->NbEq;
 459     M.NbColumns = (*P)->Dimension+2;
 460     M.p_Init = (*P)->p_Init;
 461     M.p = (*P)->Constraint;
 462
 463     Matrix *T = compress_variables(&M, nparam);
 464
 465     if (!T) {
 466         *P = NULL;
 467         return NULL;
 468     }
 469     if (isIdentity(T)) {
 470         Matrix_Free(T);
 471         T = NULL;
 472     } else
 473         *P = Polyhedron_Preimage(*P, T, MaxRays);
 474
 475     return T;
 476 }
 477
 478 bool EDomain::not_empty(lexmin_options *options)
 479 {
 480     Polyhedron *P = D;
 481     Polyhedron *Porig = P;
 482
 483     POL_ENSURE_VERTICES(P);
 484     if (emptyQ2(P))
 485         return false;
 486
 487     for (int i = 0; i < P->NbRays; ++i)
 488         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 489             sample = Vector_Alloc(P->Dimension + 1);
 490             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 491             return true;
 492         }
 493
 494     if (options->emptiness_check == BV_LEXMIN_EMPTINESS_CHECK_COUNT) {
 495         bool notzero;
 496         Value cb;
 497         value_init(cb);
 498         barvinok_count_with_options(P, &cb, options->verify.barvinok);
 499         notzero = value_notzero_p(cb);
 500         value_clear(cb);
 501         return notzero;
 502     }
 503
 504     Matrix *T = NULL;
 505     while (P && !emptyQ2(P) && P->NbEq > 0) {
 506         Polyhedron *Q = P;
 507         Matrix *T2 = remove_equalities(&P, 0, options->verify.barvinok->MaxRays);
 508         if (!T)
 509             T = T2;
 510         else {
 511             if (T2) {
 512                 Matrix *T3 = Matrix_Alloc(T->NbRows, T2->NbColumns);
 513                 Matrix_Product(T, T2, T3);
 514                 Matrix_Free(T);
 515                 Matrix_Free(T2);
 516                 T = T3;
 517             }
 518             if (Q != Porig)
 519                 Polyhedron_Free(Q);
 520         }
 521     }
 522     if (P)
 523         sample = Polyhedron_Sample(P, options->verify.barvinok);
 524     if (sample) {
 525         if (T) {
 526             Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
 527             Matrix_Vector_Product(T, sample->p, P_sample->p);
 528             Vector_Free(sample);
 529             sample = P_sample;
 530         }
 531     }
 532     if (T) {
 533         Polyhedron_Free(P);
 534         Matrix_Free(T);
 535     }
 536
 537     if (sample)
 538         assert(in_domain(Porig, sample->p));
 539
 540     return sample != NULL;
 541 }