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         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 M;
 457     Polyhedron_Matrix_View(*P, &M, (*P)->NbEq);
 458
 459     Matrix *T = compress_variables(&M, nparam);
 460
 461     if (!T) {
 462         *P = NULL;
 463         return NULL;
 464     }
 465     if (isIdentity(T)) {
 466         Matrix_Free(T);
 467         T = NULL;
 468     } else
 469         *P = Polyhedron_Preimage(*P, T, MaxRays);
 470
 471     return T;
 472 }
 473
 474 bool EDomain::not_empty(lexmin_options *options)
 475 {
 476     Polyhedron *P = D;
 477     Polyhedron *Porig = P;
 478
 479     POL_ENSURE_VERTICES(P);
 480     if (emptyQ2(P))
 481         return false;
 482
 483     for (int i = 0; i < P->NbRays; ++i)
 484         if (value_one_p(P->Ray[i][1+P->Dimension])) {
 485             sample = Vector_Alloc(P->Dimension + 1);
 486             Vector_Copy(P->Ray[i]+1, sample->p, P->Dimension+1);
 487             return true;
 488         }
 489
 490     if (options->emptiness_check == BV_LEXMIN_EMPTINESS_CHECK_COUNT) {
 491         bool notzero;
 492         Value cb;
 493         value_init(cb);
 494         barvinok_count_with_options(P, &cb, options->verify->barvinok);
 495         notzero = value_notzero_p(cb);
 496         value_clear(cb);
 497         return notzero;
 498     }
 499
 500     Matrix *T = NULL;
 501     while (P && !emptyQ2(P) && P->NbEq > 0) {
 502         Polyhedron *Q = P;
 503         Matrix *T2 = remove_equalities(&P, 0, options->verify->barvinok->MaxRays);
 504         if (!T)
 505             T = T2;
 506         else {
 507             if (T2) {
 508                 Matrix *T3 = Matrix_Alloc(T->NbRows, T2->NbColumns);
 509                 Matrix_Product(T, T2, T3);
 510                 Matrix_Free(T);
 511                 Matrix_Free(T2);
 512                 T = T3;
 513             }
 514             if (Q != Porig)
 515                 Polyhedron_Free(Q);
 516         }
 517     }
 518     if (P)
 519         sample = Polyhedron_Sample(P, options->verify->barvinok);
 520     if (sample) {
 521         if (T) {
 522             Vector *P_sample = Vector_Alloc(Porig->Dimension + 1);
 523             Matrix_Vector_Product(T, sample->p, P_sample->p);
 524             Vector_Free(sample);
 525             sample = P_sample;
 526         }
 527     }
 528     if (T) {
 529         Polyhedron_Free(P);
 530         Matrix_Free(T);
 531     }
 532
 533     if (sample)
 534         assert(in_domain(Porig, sample->p));
 535
 536     return sample != NULL;
 537 }