util.c

   1 #include <polylib/polylibgmp.h>
   2 #include "ev_operations.h"
   3 #include <stdlib.h>
   4 #include <assert.h>
   5 #include <util.h>
   6 #include "config.h"
   7
   8 #ifndef HAVE_ENUMERATE4
   9 #define Polyhedron_Enumerate(a,b,c,d) Polyhedron_Enumerate(a,b,c)
  10 #endif
  11
  12 /* Return random value between 0 and max-1 inclusive
  13  */
  14 int random_int(int max) {
  15     return (int) (((double)(max))*rand()/(RAND_MAX+1.0));
  16 }
  17
  18 /* Inplace polarization
  19  */
  20 void Polyhedron_Polarize(Polyhedron *P)
  21 {
  22     unsigned NbRows = P->NbConstraints + P->NbRays;
  23     int i;
  24     Value **q;
  25
  26     q = (Value **)malloc(NbRows * sizeof(Value *));
  27     assert(q);
  28     for (i = 0; i < P->NbRays; ++i)
  29         q[i] = P->Ray[i];
  30     for (; i < NbRows; ++i)
  31         q[i] = P->Constraint[i-P->NbRays];
  32     P->NbConstraints = NbRows - P->NbConstraints;
  33     P->NbRays = NbRows - P->NbRays;
  34     free(P->Constraint);
  35     P->Constraint = q;
  36     P->Ray = q + P->NbConstraints;
  37 }
  38
  39 /*
  40  * Rather general polar
  41  * We can optimize it significantly if we assume that
  42  * P includes zero
  43  *
  44  * Also, we calculate the polar as defined in Schrijver
  45  * The opposite should probably work as well and would
  46  * eliminate the need for multiplying by -1
  47  */
  48 Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays)
  49 {
  50     int i;
  51     Value mone;
  52     unsigned dim = P->Dimension + 2;
  53     Matrix *M = Matrix_Alloc(P->NbRays, dim);
  54
  55     assert(M);
  56     value_init(mone);
  57     value_set_si(mone, -1);
  58     for (i = 0; i < P->NbRays; ++i) {
  59         Vector_Scale(P->Ray[i], M->p[i], mone, dim);
  60         value_multiply(M->p[i][0], M->p[i][0], mone);
  61         value_multiply(M->p[i][dim-1], M->p[i][dim-1], mone);
  62     }
  63     P = Constraints2Polyhedron(M, NbMaxRays);
  64     assert(P);
  65     Matrix_Free(M);
  66     value_clear(mone);
  67     return P;
  68 }
  69
  70 /*
  71  * Returns the supporting cone of P at the vertex with index v
  72  */
  73 Polyhedron* supporting_cone(Polyhedron *P, int v)
  74 {
  75     Matrix *M;
  76     Value tmp;
  77     int i, n, j;
  78     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
  79     unsigned dim = P->Dimension + 2;
  80
  81     assert(v >=0 && v < P->NbRays);
  82     assert(value_pos_p(P->Ray[v][dim-1]));
  83     assert(supporting);
  84
  85     value_init(tmp);
  86     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
  87         Inner_Product(P->Constraint[i] + 1, P->Ray[v] + 1, dim - 1, &tmp);
  88         if ((supporting[i] = value_zero_p(tmp)))
  89             ++n;
  90     }
  91     assert(n >= dim - 2);
  92     value_clear(tmp);
  93     M = Matrix_Alloc(n, dim);
  94     assert(M);
  95     for (i = 0, j = 0; i < P->NbConstraints; ++i)
  96         if (supporting[i]) {
  97             value_set_si(M->p[j][dim-1], 0);
  98             Vector_Copy(P->Constraint[i], M->p[j++], dim-1);
  99         }
 100     free(supporting);
 101     P = Constraints2Polyhedron(M, P->NbRays+1);
 102     assert(P);
 103     Matrix_Free(M);
 104     return P;
 105 }
 106
 107 void value_lcm(Value i, Value j, Value* lcm)
 108 {
 109     Value aux;
 110     value_init(aux);
 111     value_multiply(aux,i,j);
 112     Gcd(i,j,lcm);
 113     value_division(*lcm,aux,*lcm);
 114     value_clear(aux);
 115 }
 116
 117 Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v)
 118 {
 119     Matrix *M;
 120     Value lcm, tmp, tmp2;
 121     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
 122     unsigned dim = P->Dimension + 2;
 123     unsigned nparam = v->Vertex->NbColumns - 2;
 124     unsigned nvar = dim - nparam - 2;
 125     int i, n, j;
 126     Vector *row;
 127
 128     assert(supporting);
 129     row = Vector_Alloc(nparam+1);
 130     assert(row);
 131     value_init(lcm);
 132     value_init(tmp);
 133     value_init(tmp2);
 134     value_set_si(lcm, 1);
 135     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
 136         Vector_Set(row->p, 0, nparam+1);
 137         for (j = 0 ; j < nvar; ++j) {
 138             value_set_si(tmp, 1);
 139             value_assign(tmp2,  P->Constraint[i][j+1]);
 140             if (value_ne(lcm, v->Vertex->p[j][nparam+1])) {
 141                 value_assign(tmp, lcm);
 142                 value_lcm(lcm, v->Vertex->p[j][nparam+1], &lcm);
 143                 value_division(tmp, lcm, tmp);
 144                 value_multiply(tmp2, tmp2, lcm);
 145                 value_division(tmp2, tmp2, v->Vertex->p[j][nparam+1]);
 146             }
 147             Vector_Combine(row->p, v->Vertex->p[j], row->p,
 148                            tmp, tmp2, nparam+1);
 149         }
 150         value_set_si(tmp, 1);
 151         Vector_Combine(row->p, P->Constraint[i]+1+nvar, row->p, tmp, lcm, nparam+1);
 152         for (j = 0; j < nparam+1; ++j)
 153             if (value_notzero_p(row->p[j]))
 154                 break;
 155         if ((supporting[i] = (j == nparam + 1)))
 156             ++n;
 157     }
 158     assert(n >= nvar);
 159     value_clear(tmp);
 160     value_clear(tmp2);
 161     value_clear(lcm);
 162     Vector_Free(row);
 163     M = Matrix_Alloc(n, nvar+2);
 164     assert(M);
 165     for (i = 0, j = 0; i < P->NbConstraints; ++i)
 166         if (supporting[i]) {
 167             value_set_si(M->p[j][nvar+1], 0);
 168             Vector_Copy(P->Constraint[i], M->p[j++], nvar+1);
 169         }
 170     free(supporting);
 171     P = Constraints2Polyhedron(M, P->NbRays+1);
 172     assert(P);
 173     Matrix_Free(M);
 174     return P;
 175 }
 176
 177 Polyhedron* triangularize_cone(Polyhedron *P, unsigned NbMaxCons)
 178 {
 179     const static int MAX_TRY=10;
 180     int i, j, r, n, t;
 181     Value tmp;
 182     unsigned dim = P->Dimension;
 183     Matrix *M = Matrix_Alloc(P->NbRays+1, dim+3);
 184     Matrix *M2;
 185     Polyhedron *L, *R, *T;
 186     assert(P->NbEq == 0);
 187
 188     R = NULL;
 189     value_init(tmp);
 190
 191     Vector_Set(M->p[0]+1, 0, dim+1);
 192     value_set_si(M->p[0][0], 1);
 193     value_set_si(M->p[0][dim+2], 1);
 194     Vector_Set(M->p[P->NbRays]+1, 0, dim+2);
 195     value_set_si(M->p[P->NbRays][0], 1);
 196     value_set_si(M->p[P->NbRays][dim+1], 1);
 197
 198     for (i = 0, r = 1; i < P->NbRays; ++i) {
 199         if (value_notzero_p(P->Ray[i][dim+1]))
 200             continue;
 201         Vector_Copy(P->Ray[i], M->p[r], dim+1);
 202         Inner_Product(M->p[r]+1, M->p[r]+1, dim, &tmp);
 203         value_assign(M->p[r][dim+1], tmp);
 204         value_set_si(M->p[r][dim+2], 0);
 205         ++r;
 206     }
 207
 208     L = Rays2Polyhedron(M, NbMaxCons);
 209
 210     M2 = Matrix_Alloc(dim+1, dim+2);
 211     Vector_Set(M2->p[0]+1, 0, dim);
 212     value_set_si(M2->p[0][0], 1);
 213     value_set_si(M2->p[0][dim+1], 1);
 214
 215     t = 1;
 216     if (0) {
 217 try_again:
 218         /* Usually R should still be 0 */
 219         Domain_Free(R);
 220         Polyhedron_Free(L);
 221         for (r = 1; r < P->NbRays; ++r) {
 222             value_set_si(M->p[r][dim+1], random_int((t+1)*dim)+1);
 223         }
 224         L = Rays2Polyhedron(M, NbMaxCons);
 225         ++t;
 226     }
 227     assert(t <= MAX_TRY);
 228
 229     R = NULL;
 230     n = 0;
 231
 232     for (i = 0; i < L->NbConstraints; ++i) {
 233         if (value_negz_p(L->Constraint[i][dim+1]))
 234             continue;
 235         if (value_notzero_p(L->Constraint[i][dim+2]))
 236             continue;
 237         for (j = 1, r = 1; j < M->NbRows; ++j) {
 238             Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
 239             if (value_notzero_p(tmp))
 240                 continue;
 241             if (r > dim)
 242                 goto try_again;
 243             Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
 244             value_set_si(M2->p[r][0], 1);
 245             value_set_si(M2->p[r][dim+1], 0);
 246             ++r;
 247         }
 248         assert(r == dim+1);
 249         T = Rays2Polyhedron(M2, P->NbConstraints+1);
 250         T->next = R;
 251         R = T;
 252         ++n;
 253     }
 254     Matrix_Free(M2);
 255
 256     Polyhedron_Free(L);
 257     value_clear(tmp);
 258     Matrix_Free(M);
 259
 260     return R;
 261 }
 262
 263 void check_triangulization(Polyhedron *P, Polyhedron *T)
 264 {
 265     Polyhedron *C, *D, *E, *F, *G, *U;
 266     for (C = T; C; C = C->next) {
 267         if (C == T)
 268             U = C;
 269         else
 270             U = DomainConvex(DomainUnion(U, C, 100), 100);
 271         for (D = C->next; D; D = D->next) {
 272             F = C->next;
 273             G = D->next;
 274             C->next = NULL;
 275             D->next = NULL;
 276             E = DomainIntersection(C, D, 600);
 277             assert(E->NbRays == 0 || E->NbEq >= 1);
 278             Polyhedron_Free(E);
 279             C->next = F;
 280             D->next = G;
 281         }
 282     }
 283     assert(PolyhedronIncludes(U, P));
 284     assert(PolyhedronIncludes(P, U));
 285 }
 286
 287 void Euclid(Value a, Value b, Value *x, Value *y, Value *g)
 288 {
 289     Value c, d, e, f, tmp;
 290
 291     value_init(c);
 292     value_init(d);
 293     value_init(e);
 294     value_init(f);
 295     value_init(tmp);
 296     value_absolute(c, a);
 297     value_absolute(d, b);
 298     value_set_si(e, 1);
 299     value_set_si(f, 0);
 300     while(value_pos_p(d)) {
 301         value_division(tmp, c, d);
 302         value_multiply(tmp, tmp, f);
 303         value_substract(e, e, tmp);
 304         value_division(tmp, c, d);
 305         value_multiply(tmp, tmp, d);
 306         value_substract(c, c, tmp);
 307         value_swap(c, d);
 308         value_swap(e, f);
 309     }
 310     value_assign(*g, c);
 311     if (value_zero_p(a))
 312         value_set_si(*x, 0);
 313     else if (value_pos_p(a))
 314         value_assign(*x, e);
 315     else value_oppose(*x, e);
 316     if (value_zero_p(b))
 317         value_set_si(*y, 0);
 318     else {
 319         value_multiply(tmp, a, *x);
 320         value_substract(tmp, c, tmp);
 321         value_division(*y, tmp, b);
 322     }
 323     value_clear(c);
 324     value_clear(d);
 325     value_clear(e);
 326     value_clear(f);
 327     value_clear(tmp);
 328 }
 329
 330 Matrix * unimodular_complete(Vector *row)
 331 {
 332     Value g, b, c, old, tmp;
 333     Matrix *m;
 334     unsigned i, j;
 335
 336     value_init(b);
 337     value_init(c);
 338     value_init(g);
 339     value_init(old);
 340     value_init(tmp);
 341     m = Matrix_Alloc(row->Size, row->Size);
 342     for (j = 0; j < row->Size; ++j) {
 343         value_assign(m->p[0][j], row->p[j]);
 344     }
 345     value_assign(g, row->p[0]);
 346     for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
 347         for (j = 0; j < row->Size; ++j) {
 348             if (j == i-1)
 349                 value_set_si(m->p[i][j], 1);
 350             else
 351                 value_set_si(m->p[i][j], 0);
 352         }
 353         value_assign(g, row->p[i]);
 354     }
 355     for (; i < row->Size; ++i) {
 356         value_assign(old, g);
 357         Euclid(old, row->p[i], &c, &b, &g);
 358         value_oppose(b, b);
 359         for (j = 0; j < row->Size; ++j) {
 360             if (j < i) {
 361                 value_multiply(tmp, row->p[j], b);
 362                 value_division(m->p[i][j], tmp, old);
 363             } else if (j == i)
 364                 value_assign(m->p[i][j], c);
 365             else
 366                 value_set_si(m->p[i][j], 0);
 367         }
 368     }
 369     value_clear(b);
 370     value_clear(c);
 371     value_clear(g);
 372     value_clear(old);
 373     value_clear(tmp);
 374     return m;
 375 }
 376
 377 /*
 378  * Returns a full-dimensional polyhedron with the same number
 379  * of integer points as P
 380  */
 381 Polyhedron *remove_equalities(Polyhedron *P)
 382 {
 383     Value g;
 384     Vector *v;
 385     Polyhedron *p = Polyhedron_Copy(P), *q;
 386     unsigned dim = p->Dimension;
 387     Matrix *m1, *m2;
 388     int i;
 389
 390     value_init(g);
 391     while (p->NbEq > 0) {
 392         assert(dim > 0);
 393         Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
 394         Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
 395         Vector_Gcd(p->Constraint[0]+1, dim, &g);
 396         if (value_notone_p(g) && value_notmone_p(g)) {
 397             Polyhedron_Free(p);
 398             p = Empty_Polyhedron(0);
 399             break;
 400         }
 401         v = Vector_Alloc(dim);
 402         Vector_Copy(p->Constraint[0]+1, v->p, dim);
 403         m1 = unimodular_complete(v);
 404         m2 = Matrix_Alloc(dim, dim+1);
 405         for (i = 0; i < dim-1 ; ++i) {
 406             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 407             value_set_si(m2->p[i][dim], 0);
 408         }
 409         Vector_Set(m2->p[dim-1], 0, dim);
 410         value_set_si(m2->p[dim-1][dim], 1);
 411         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 412         Vector_Free(v);
 413         Matrix_Free(m1);
 414         Matrix_Free(m2);
 415         Polyhedron_Free(p);
 416         p = q;
 417         --dim;
 418     }
 419     value_clear(g);
 420     return p;
 421 }
 422
 423 /*
 424  * Returns a full-dimensional polyhedron with the same number
 425  * of integer points as P
 426  * nvar specifies the number of variables
 427  * The remaining dimensions are assumed to be parameters
 428  * Destroys P
 429  */
 430 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
 431 {
 432     Value g;
 433     Vector *v;
 434     Polyhedron *p = P, *q;
 435     unsigned dim = p->Dimension;
 436     Matrix *m1, *m2, *f;
 437     int i, j, skip;
 438
 439     value_init(g);
 440     f = Matrix_Alloc(p->NbEq, dim-nvar+2);
 441     j = 0;
 442     *factor = f;
 443     skip = 0;
 444     while (nvar > 0 && p->NbEq - skip > 0) {
 445         assert(dim > 0);
 446
 447         while (value_zero_p(p->Constraint[skip][0]) &&
 448                First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
 449             ++skip;
 450         if (p->NbEq == skip)
 451             break;
 452
 453         Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
 454         Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
 455         Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
 456         Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
 457         value_assign(f->p[j][dim-nvar+1], g);
 458         v = Vector_Alloc(dim);
 459         Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
 460         Vector_Set(v->p+nvar, 0, dim-nvar);
 461         m1 = unimodular_complete(v);
 462         m2 = Matrix_Alloc(dim, dim+1);
 463         for (i = 0; i < dim-1 ; ++i) {
 464             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 465             value_set_si(m2->p[i][dim], 0);
 466         }
 467         Vector_Set(m2->p[dim-1], 0, dim);
 468         value_set_si(m2->p[dim-1][dim], 1);
 469         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 470         Vector_Free(v);
 471         Matrix_Free(m1);
 472         Matrix_Free(m2);
 473         Polyhedron_Free(p);
 474         p = q;
 475         --dim;
 476         --nvar;
 477         ++j;
 478     }
 479     value_clear(g);
 480     return p;
 481 }
 482
 483 struct single {
 484     int nr;
 485     int pos[2];
 486 };
 487
 488 /*
 489  * The number of points in P is equal to factor time
 490  * the number of points in the polyhedron returned.
 491  * The return value is zero if no reduction can be found.
 492  */
 493 Polyhedron* Polyhedron_Reduce(Polyhedron *P, Value* factor)
 494 {
 495     int i, j, prev, nsingle, k, p;
 496     unsigned dim = P->Dimension;
 497     struct single *singles;
 498     int *bad;
 499     Value tmp, pos, neg;
 500
 501     value_init(tmp);
 502     value_init(pos);
 503     value_init(neg);
 504
 505     value_set_si(*factor, 1);
 506
 507     singles = (struct single *)malloc(dim * sizeof(struct single));
 508     assert(singles);
 509     for (i = 0; i < dim; ++i)
 510         singles[i].nr = 0;
 511     bad = (int *)calloc(dim, sizeof(int));
 512     assert(bad);
 513
 514     assert (P->NbEq == 0);
 515
 516     nsingle = 0;
 517     for (i = 0; i < P->NbConstraints; ++i) {
 518         for (j = 0, prev = -1; j < dim; ++j) {
 519             if (value_notzero_p(P->Constraint[i][j+1])) {
 520                 if (prev == -1)
 521                     prev = j;
 522                 else {
 523                     if (prev != -2)
 524                         bad[prev] = 1;
 525                     bad[j] = 1;
 526                     prev = -2;
 527                 }
 528             }
 529         }
 530         if (prev >= 0)
 531             singles[prev].pos[singles[prev].nr++] = i;
 532     }
 533     for (j = 0; j < dim; ++j) {
 534         if (bad[j])
 535             singles[j].nr = 0;
 536         else if (singles[j].nr == 2)
 537             ++nsingle;
 538     }
 539     if (nsingle) {
 540         Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
 541         for (i = 0, j = 0; i < dim; ++i) {
 542             if (singles[i].nr != 2)
 543                 value_set_si(m->p[j++][i], 1);
 544             else {
 545                 for (k = 0; k <= 1; ++k) {
 546                     p = singles[i].pos[k];
 547                     value_oppose(tmp, P->Constraint[p][dim+1]);
 548                     if (value_pos_p(P->Constraint[p][i+1]))
 549                         mpz_cdiv_q(pos, tmp, P->Constraint[p][i+1]);
 550                     else
 551                         mpz_fdiv_q(neg, tmp, P->Constraint[p][i+1]);
 552                 }
 553                 value_substract(tmp, neg, pos);
 554                 value_increment(tmp, tmp);
 555                 value_multiply(*factor, *factor, tmp);
 556             }
 557         }
 558         value_set_si(m->p[dim-nsingle][dim], 1);
 559         P = Polyhedron_Image(P, m, P->NbConstraints);
 560         Matrix_Free(m);
 561     } else
 562         P = NULL;
 563     free(singles);
 564     free(bad);
 565
 566     value_clear(tmp);
 567     value_clear(pos);
 568     value_clear(neg);
 569
 570     return P;
 571 }
 572
 573 void manual_count(Polyhedron *P, Value* result)
 574 {
 575     Polyhedron *U = Universe_Polyhedron(0);
 576     Enumeration *en = Polyhedron_Enumerate(P,U,1024,NULL);
 577     Value *v = compute_poly(en,NULL);
 578     value_assign(*result, *v);
 579     value_clear(*v);
 580     free(v);
 581     Enumeration_Free(en);
 582     Polyhedron_Free(U);
 583 }
 584
 585 Bool isIdentity(Matrix *M)
 586 {
 587     unsigned i, j;
 588     if (M->NbRows != M->NbColumns)
 589         return False;
 590
 591     for (i = 0;i < M->NbRows; i ++)
 592         for (j = 0; j < M->NbColumns; j ++)
 593             if (i == j) {
 594                 if(value_notone_p(M->p[i][j]))
 595                     return False;
 596             } else {
 597                 if(value_notzero_p(M->p[i][j]))
 598                     return False;
 599             }
 600     return True;
 601 }
 602
 603 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
 604 {
 605   Param_Domain *P;
 606   Param_Vertices *V;
 607
 608   for(P=PP->D;P;P=P->next) {
 609
 610     /* prints current val. dom. */
 611     printf( "---------------------------------------\n" );
 612     printf( "Domain :\n");
 613     Print_Domain( stdout, P->Domain, param_names );
 614
 615     /* scan the vertices */
 616     printf( "Vertices :\n");
 617     FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
 618
 619       /* prints each vertex */
 620       Print_Vertex( stdout, V->Vertex, param_names );
 621       printf( "\n" );
 622     }
 623     END_FORALL_PVertex_in_ParamPolyhedron;
 624   }
 625 }
 626
 627 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
 628 {
 629     for (; en; en = en->next) {
 630         Print_Domain(Dst, en->ValidityDomain, params);
 631         print_evalue(Dst, &en->EP, params);
 632     }
 633 }
 634
 635 void Free_ParamNames(char **params, int m)
 636 {
 637     while (--m >= 0)
 638         free(params[m]);
 639     free(params);
 640 }