util.c

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