edomain.cc

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