util.c

   1 #include <polylib/polylibgmp.h>
   2 #include <stdlib.h>
   3 #include <assert.h>
   4 #include "config.h"
   5
   6 #ifndef HAVE_ENUMERATE4
   7 #define Polyhedron_Enumerate(a,b,c,d) Polyhedron_Enumerate(a,b,c)
   8 #endif
   9
  10 void manual_count(Polyhedron *P, Value* result)
  11 {
  12     Polyhedron *U = Universe_Polyhedron(0);
  13     Enumeration *en = Polyhedron_Enumerate(P,U,1024,NULL);
  14     Value *v = compute_poly(en,NULL);
  15     value_assign(*result, *v);
  16     value_clear(*v);
  17     free(v);
  18     Enumeration_Free(en);
  19     Polyhedron_Free(U);
  20 }
  21
  22 #include "ev_operations.h"
  23 #include <util.h>
  24 #include <barvinok.h>
  25
  26 /* Return random value between 0 and max-1 inclusive
  27  */
  28 int random_int(int max) {
  29     return (int) (((double)(max))*rand()/(RAND_MAX+1.0));
  30 }
  31
  32 /* Inplace polarization
  33  */
  34 void Polyhedron_Polarize(Polyhedron *P)
  35 {
  36     unsigned NbRows = P->NbConstraints + P->NbRays;
  37     int i;
  38     Value **q;
  39
  40     q = (Value **)malloc(NbRows * sizeof(Value *));
  41     assert(q);
  42     for (i = 0; i < P->NbRays; ++i)
  43         q[i] = P->Ray[i];
  44     for (; i < NbRows; ++i)
  45         q[i] = P->Constraint[i-P->NbRays];
  46     P->NbConstraints = NbRows - P->NbConstraints;
  47     P->NbRays = NbRows - P->NbRays;
  48     free(P->Constraint);
  49     P->Constraint = q;
  50     P->Ray = q + P->NbConstraints;
  51 }
  52
  53 /*
  54  * Rather general polar
  55  * We can optimize it significantly if we assume that
  56  * P includes zero
  57  *
  58  * Also, we calculate the polar as defined in Schrijver
  59  * The opposite should probably work as well and would
  60  * eliminate the need for multiplying by -1
  61  */
  62 Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays)
  63 {
  64     int i;
  65     Value mone;
  66     unsigned dim = P->Dimension + 2;
  67     Matrix *M = Matrix_Alloc(P->NbRays, dim);
  68
  69     assert(M);
  70     value_init(mone);
  71     value_set_si(mone, -1);
  72     for (i = 0; i < P->NbRays; ++i) {
  73         Vector_Scale(P->Ray[i], M->p[i], mone, dim);
  74         value_multiply(M->p[i][0], M->p[i][0], mone);
  75         value_multiply(M->p[i][dim-1], M->p[i][dim-1], mone);
  76     }
  77     P = Constraints2Polyhedron(M, NbMaxRays);
  78     assert(P);
  79     Matrix_Free(M);
  80     value_clear(mone);
  81     return P;
  82 }
  83
  84 /*
  85  * Returns the supporting cone of P at the vertex with index v
  86  */
  87 Polyhedron* supporting_cone(Polyhedron *P, int v)
  88 {
  89     Matrix *M;
  90     Value tmp;
  91     int i, n, j;
  92     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
  93     unsigned dim = P->Dimension + 2;
  94
  95     assert(v >=0 && v < P->NbRays);
  96     assert(value_pos_p(P->Ray[v][dim-1]));
  97     assert(supporting);
  98
  99     value_init(tmp);
 100     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
 101         Inner_Product(P->Constraint[i] + 1, P->Ray[v] + 1, dim - 1, &tmp);
 102         if ((supporting[i] = value_zero_p(tmp)))
 103             ++n;
 104     }
 105     assert(n >= dim - 2);
 106     value_clear(tmp);
 107     M = Matrix_Alloc(n, dim);
 108     assert(M);
 109     for (i = 0, j = 0; i < P->NbConstraints; ++i)
 110         if (supporting[i]) {
 111             value_set_si(M->p[j][dim-1], 0);
 112             Vector_Copy(P->Constraint[i], M->p[j++], dim-1);
 113         }
 114     free(supporting);
 115     P = Constraints2Polyhedron(M, P->NbRays+1);
 116     assert(P);
 117     Matrix_Free(M);
 118     return P;
 119 }
 120
 121 void value_lcm(Value i, Value j, Value* lcm)
 122 {
 123     Value aux;
 124     value_init(aux);
 125     value_multiply(aux,i,j);
 126     Gcd(i,j,lcm);
 127     value_division(*lcm,aux,*lcm);
 128     value_clear(aux);
 129 }
 130
 131 Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v)
 132 {
 133     Matrix *M;
 134     Value lcm, tmp, tmp2;
 135     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
 136     unsigned dim = P->Dimension + 2;
 137     unsigned nparam = v->Vertex->NbColumns - 2;
 138     unsigned nvar = dim - nparam - 2;
 139     int i, n, j;
 140     Vector *row;
 141
 142     assert(supporting);
 143     row = Vector_Alloc(nparam+1);
 144     assert(row);
 145     value_init(lcm);
 146     value_init(tmp);
 147     value_init(tmp2);
 148     value_set_si(lcm, 1);
 149     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
 150         Vector_Set(row->p, 0, nparam+1);
 151         for (j = 0 ; j < nvar; ++j) {
 152             value_set_si(tmp, 1);
 153             value_assign(tmp2,  P->Constraint[i][j+1]);
 154             if (value_ne(lcm, v->Vertex->p[j][nparam+1])) {
 155                 value_assign(tmp, lcm);
 156                 value_lcm(lcm, v->Vertex->p[j][nparam+1], &lcm);
 157                 value_division(tmp, lcm, tmp);
 158                 value_multiply(tmp2, tmp2, lcm);
 159                 value_division(tmp2, tmp2, v->Vertex->p[j][nparam+1]);
 160             }
 161             Vector_Combine(row->p, v->Vertex->p[j], row->p,
 162                            tmp, tmp2, nparam+1);
 163         }
 164         value_set_si(tmp, 1);
 165         Vector_Combine(row->p, P->Constraint[i]+1+nvar, row->p, tmp, lcm, nparam+1);
 166         for (j = 0; j < nparam+1; ++j)
 167             if (value_notzero_p(row->p[j]))
 168                 break;
 169         if ((supporting[i] = (j == nparam + 1)))
 170             ++n;
 171     }
 172     assert(n >= nvar);
 173     value_clear(tmp);
 174     value_clear(tmp2);
 175     value_clear(lcm);
 176     Vector_Free(row);
 177     M = Matrix_Alloc(n, nvar+2);
 178     assert(M);
 179     for (i = 0, j = 0; i < P->NbConstraints; ++i)
 180         if (supporting[i]) {
 181             value_set_si(M->p[j][nvar+1], 0);
 182             Vector_Copy(P->Constraint[i], M->p[j++], nvar+1);
 183         }
 184     free(supporting);
 185     P = Constraints2Polyhedron(M, P->NbRays+1);
 186     assert(P);
 187     Matrix_Free(M);
 188     return P;
 189 }
 190
 191 Polyhedron* triangularize_cone(Polyhedron *P, unsigned NbMaxCons)
 192 {
 193     const static int MAX_TRY=10;
 194     int i, j, r, n, t;
 195     Value tmp;
 196     unsigned dim = P->Dimension;
 197     Matrix *M = Matrix_Alloc(P->NbRays+1, dim+3);
 198     Matrix *M2;
 199     Polyhedron *L, *R, *T;
 200     assert(P->NbEq == 0);
 201
 202     R = NULL;
 203     value_init(tmp);
 204
 205     Vector_Set(M->p[0]+1, 0, dim+1);
 206     value_set_si(M->p[0][0], 1);
 207     value_set_si(M->p[0][dim+2], 1);
 208     Vector_Set(M->p[P->NbRays]+1, 0, dim+2);
 209     value_set_si(M->p[P->NbRays][0], 1);
 210     value_set_si(M->p[P->NbRays][dim+1], 1);
 211
 212     for (i = 0, r = 1; i < P->NbRays; ++i) {
 213         if (value_notzero_p(P->Ray[i][dim+1]))
 214             continue;
 215         Vector_Copy(P->Ray[i], M->p[r], dim+1);
 216         Inner_Product(M->p[r]+1, M->p[r]+1, dim, &tmp);
 217         value_assign(M->p[r][dim+1], tmp);
 218         value_set_si(M->p[r][dim+2], 0);
 219         ++r;
 220     }
 221
 222     L = Rays2Polyhedron(M, NbMaxCons);
 223
 224     M2 = Matrix_Alloc(dim+1, dim+2);
 225     Vector_Set(M2->p[0]+1, 0, dim);
 226     value_set_si(M2->p[0][0], 1);
 227     value_set_si(M2->p[0][dim+1], 1);
 228
 229     t = 1;
 230     if (0) {
 231 try_again:
 232         /* Usually R should still be 0 */
 233         Domain_Free(R);
 234         Polyhedron_Free(L);
 235         for (r = 1; r < P->NbRays; ++r) {
 236             value_set_si(M->p[r][dim+1], random_int((t+1)*dim)+1);
 237         }
 238         L = Rays2Polyhedron(M, NbMaxCons);
 239         ++t;
 240     }
 241     assert(t <= MAX_TRY);
 242
 243     R = NULL;
 244     n = 0;
 245
 246     for (i = 0; i < L->NbConstraints; ++i) {
 247         if (value_negz_p(L->Constraint[i][dim+1]))
 248             continue;
 249         if (value_notzero_p(L->Constraint[i][dim+2]))
 250             continue;
 251         for (j = 1, r = 1; j < M->NbRows; ++j) {
 252             Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
 253             if (value_notzero_p(tmp))
 254                 continue;
 255             if (r > dim)
 256                 goto try_again;
 257             Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
 258             value_set_si(M2->p[r][0], 1);
 259             value_set_si(M2->p[r][dim+1], 0);
 260             ++r;
 261         }
 262         assert(r == dim+1);
 263         T = Rays2Polyhedron(M2, P->NbConstraints+1);
 264         T->next = R;
 265         R = T;
 266         ++n;
 267     }
 268     Matrix_Free(M2);
 269
 270     Polyhedron_Free(L);
 271     value_clear(tmp);
 272     Matrix_Free(M);
 273
 274     return R;
 275 }
 276
 277 void check_triangulization(Polyhedron *P, Polyhedron *T)
 278 {
 279     Polyhedron *C, *D, *E, *F, *G, *U;
 280     for (C = T; C; C = C->next) {
 281         if (C == T)
 282             U = C;
 283         else
 284             U = DomainConvex(DomainUnion(U, C, 100), 100);
 285         for (D = C->next; D; D = D->next) {
 286             F = C->next;
 287             G = D->next;
 288             C->next = NULL;
 289             D->next = NULL;
 290             E = DomainIntersection(C, D, 600);
 291             assert(E->NbRays == 0 || E->NbEq >= 1);
 292             Polyhedron_Free(E);
 293             C->next = F;
 294             D->next = G;
 295         }
 296     }
 297     assert(PolyhedronIncludes(U, P));
 298     assert(PolyhedronIncludes(P, U));
 299 }
 300
 301 void Euclid(Value a, Value b, Value *x, Value *y, Value *g)
 302 {
 303     Value c, d, e, f, tmp;
 304
 305     value_init(c);
 306     value_init(d);
 307     value_init(e);
 308     value_init(f);
 309     value_init(tmp);
 310     value_absolute(c, a);
 311     value_absolute(d, b);
 312     value_set_si(e, 1);
 313     value_set_si(f, 0);
 314     while(value_pos_p(d)) {
 315         value_division(tmp, c, d);
 316         value_multiply(tmp, tmp, f);
 317         value_substract(e, e, tmp);
 318         value_division(tmp, c, d);
 319         value_multiply(tmp, tmp, d);
 320         value_substract(c, c, tmp);
 321         value_swap(c, d);
 322         value_swap(e, f);
 323     }
 324     value_assign(*g, c);
 325     if (value_zero_p(a))
 326         value_set_si(*x, 0);
 327     else if (value_pos_p(a))
 328         value_assign(*x, e);
 329     else value_oppose(*x, e);
 330     if (value_zero_p(b))
 331         value_set_si(*y, 0);
 332     else {
 333         value_multiply(tmp, a, *x);
 334         value_substract(tmp, c, tmp);
 335         value_division(*y, tmp, b);
 336     }
 337     value_clear(c);
 338     value_clear(d);
 339     value_clear(e);
 340     value_clear(f);
 341     value_clear(tmp);
 342 }
 343
 344 Matrix * unimodular_complete(Vector *row)
 345 {
 346     Value g, b, c, old, tmp;
 347     Matrix *m;
 348     unsigned i, j;
 349
 350     value_init(b);
 351     value_init(c);
 352     value_init(g);
 353     value_init(old);
 354     value_init(tmp);
 355     m = Matrix_Alloc(row->Size, row->Size);
 356     for (j = 0; j < row->Size; ++j) {
 357         value_assign(m->p[0][j], row->p[j]);
 358     }
 359     value_assign(g, row->p[0]);
 360     for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
 361         for (j = 0; j < row->Size; ++j) {
 362             if (j == i-1)
 363                 value_set_si(m->p[i][j], 1);
 364             else
 365                 value_set_si(m->p[i][j], 0);
 366         }
 367         value_assign(g, row->p[i]);
 368     }
 369     for (; i < row->Size; ++i) {
 370         value_assign(old, g);
 371         Euclid(old, row->p[i], &c, &b, &g);
 372         value_oppose(b, b);
 373         for (j = 0; j < row->Size; ++j) {
 374             if (j < i) {
 375                 value_multiply(tmp, row->p[j], b);
 376                 value_division(m->p[i][j], tmp, old);
 377             } else if (j == i)
 378                 value_assign(m->p[i][j], c);
 379             else
 380                 value_set_si(m->p[i][j], 0);
 381         }
 382     }
 383     value_clear(b);
 384     value_clear(c);
 385     value_clear(g);
 386     value_clear(old);
 387     value_clear(tmp);
 388     return m;
 389 }
 390
 391 /*
 392  * Returns a full-dimensional polyhedron with the same number
 393  * of integer points as P
 394  */
 395 Polyhedron *remove_equalities(Polyhedron *P)
 396 {
 397     Value g;
 398     Vector *v;
 399     Polyhedron *p = Polyhedron_Copy(P), *q;
 400     unsigned dim = p->Dimension;
 401     Matrix *m1, *m2;
 402     int i;
 403
 404     value_init(g);
 405     while (p->NbEq > 0) {
 406         assert(dim > 0);
 407         Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
 408         Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
 409         Vector_Gcd(p->Constraint[0]+1, dim, &g);
 410         if (value_notone_p(g) && value_notmone_p(g)) {
 411             Polyhedron_Free(p);
 412             p = Empty_Polyhedron(0);
 413             break;
 414         }
 415         v = Vector_Alloc(dim);
 416         Vector_Copy(p->Constraint[0]+1, v->p, dim);
 417         m1 = unimodular_complete(v);
 418         m2 = Matrix_Alloc(dim, dim+1);
 419         for (i = 0; i < dim-1 ; ++i) {
 420             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 421             value_set_si(m2->p[i][dim], 0);
 422         }
 423         Vector_Set(m2->p[dim-1], 0, dim);
 424         value_set_si(m2->p[dim-1][dim], 1);
 425         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 426         Vector_Free(v);
 427         Matrix_Free(m1);
 428         Matrix_Free(m2);
 429         Polyhedron_Free(p);
 430         p = q;
 431         --dim;
 432     }
 433     value_clear(g);
 434     return p;
 435 }
 436
 437 /*
 438  * Returns a full-dimensional polyhedron with the same number
 439  * of integer points as P
 440  * nvar specifies the number of variables
 441  * The remaining dimensions are assumed to be parameters
 442  * Destroys P
 443  * factor is NbEq x (nparam+2) matrix, containing stride constraints
 444  * on the parameters; column nparam is the constant;
 445  * column nparam+1 is the stride
 446  */
 447 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
 448 {
 449     Value g;
 450     Vector *v;
 451     Polyhedron *p = P, *q;
 452     unsigned dim = p->Dimension;
 453     Matrix *m1, *m2, *f;
 454     int i, j, skip;
 455
 456     value_init(g);
 457     f = Matrix_Alloc(p->NbEq, dim-nvar+2);
 458     j = 0;
 459     *factor = f;
 460     skip = 0;
 461     while (nvar > 0 && p->NbEq - skip > 0) {
 462         assert(dim > 0);
 463
 464         while (value_zero_p(p->Constraint[skip][0]) &&
 465                First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
 466             ++skip;
 467         if (p->NbEq == skip)
 468             break;
 469
 470         Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
 471         Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
 472         Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
 473         Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
 474         value_assign(f->p[j][dim-nvar+1], g);
 475         v = Vector_Alloc(dim);
 476         Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
 477         Vector_Set(v->p+nvar, 0, dim-nvar);
 478         m1 = unimodular_complete(v);
 479         m2 = Matrix_Alloc(dim, dim+1);
 480         for (i = 0; i < dim-1 ; ++i) {
 481             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 482             value_set_si(m2->p[i][dim], 0);
 483         }
 484         Vector_Set(m2->p[dim-1], 0, dim);
 485         value_set_si(m2->p[dim-1][dim], 1);
 486         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 487         Vector_Free(v);
 488         Matrix_Free(m1);
 489         Matrix_Free(m2);
 490         Polyhedron_Free(p);
 491         p = q;
 492         --dim;
 493         --nvar;
 494         ++j;
 495     }
 496     value_clear(g);
 497     return p;
 498 }
 499
 500 struct single {
 501     int nr;
 502     int pos[2];
 503 };
 504
 505 static void free_singles(int **singles, int dim)
 506 {
 507     int i;
 508     for (i = 0; i < dim; ++i)
 509         free(singles[i]);
 510     free(singles);
 511 }
 512
 513 static int **find_singles(Polyhedron *P, int dim, int max, int *nsingle)
 514 {
 515     int i, j, prev;
 516     int **singles = (int **) malloc(dim * sizeof(int *));
 517     assert(singles);
 518
 519     for (i = 0; i < dim; ++i) {
 520         singles[i] = (int *) malloc((max + 1) *sizeof(int));
 521         singles[i][0] = 0;
 522     }
 523
 524     for (i = 0; i < P->NbConstraints; ++i) {
 525         for (j = 0, prev = -1; j < dim; ++j) {
 526             if (value_notzero_p(P->Constraint[i][j+1])) {
 527                 if (prev == -1)
 528                     prev = j;
 529                 else {
 530                     if (prev != -2)
 531                         singles[prev][0] = -1;
 532                     singles[j][0] = -1;
 533                     prev = -2;
 534                 }
 535             }
 536         }
 537         if (prev >= 0 && singles[prev][0] >= 0)
 538             singles[prev][++singles[prev][0]] = i;
 539     }
 540     *nsingle = 0;
 541     for (j = 0; j < dim; ++j)
 542         if (singles[j][0] > 0)
 543             ++*nsingle;
 544     if (!*nsingle) {
 545         free_singles(singles, dim);
 546         singles = 0;
 547     }
 548     return singles;
 549 }
 550
 551 /*
 552  * The number of points in P is equal to factor time
 553  * the number of points in the polyhedron returned.
 554  * The return value is zero if no reduction can be found.
 555  */
 556 Polyhedron* Polyhedron_Reduce(Polyhedron *P, Value* factor)
 557 {
 558     int i, j, nsingle, k, p;
 559     unsigned dim = P->Dimension;
 560     int **singles;
 561
 562     value_set_si(*factor, 1);
 563
 564     assert (P->NbEq == 0);
 565
 566     singles = find_singles(P, dim, 2, &nsingle);
 567
 568     if (nsingle == 0)
 569         return 0;
 570
 571     {
 572         Value tmp, pos, neg;
 573         Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
 574
 575         value_init(tmp);
 576         value_init(pos);
 577         value_init(neg);
 578
 579         for (i = 0, j = 0; i < dim; ++i) {
 580             if (singles[i][0] != 2)
 581                 value_set_si(m->p[j++][i], 1);
 582             else {
 583                 for (k = 0; k <= 1; ++k) {
 584                     p = singles[i][1+k];
 585                     value_oppose(tmp, P->Constraint[p][dim+1]);
 586                     if (value_pos_p(P->Constraint[p][i+1]))
 587                         mpz_cdiv_q(pos, tmp, P->Constraint[p][i+1]);
 588                     else
 589                         mpz_fdiv_q(neg, tmp, P->Constraint[p][i+1]);
 590                 }
 591                 value_substract(tmp, neg, pos);
 592                 value_increment(tmp, tmp);
 593                 value_multiply(*factor, *factor, tmp);
 594             }
 595         }
 596         value_set_si(m->p[dim-nsingle][dim], 1);
 597         P = Polyhedron_Image(P, m, P->NbConstraints);
 598         Matrix_Free(m);
 599         free_singles(singles, dim);
 600
 601         value_clear(tmp);
 602         value_clear(pos);
 603         value_clear(neg);
 604     }
 605
 606     return P;
 607 }
 608
 609 struct section { Polyhedron * D; evalue E; };
 610
 611 static Polyhedron* ParamPolyhedron_Reduce_mod(Polyhedron *P, unsigned nvar,
 612                                               evalue* factor)
 613 {
 614     int nsingle;
 615     int **singles;
 616     unsigned dim = P->Dimension;
 617
 618     singles = find_singles(P, nvar, P->NbConstraints, &nsingle);
 619
 620     if (nsingle == 0)
 621         return 0;
 622
 623     {
 624         int i, j, p, n;
 625         Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
 626         Value tmp;
 627         evalue mone;
 628         value_init(mone.d);
 629         evalue_set_si(&mone, -1, 1);
 630
 631         value_init(tmp);
 632
 633         for (i = 0, j = 0; i < dim; ++i) {
 634             if (i >= nvar || singles[i][0] < 2)
 635                 value_set_si(m->p[j++][i], 1);
 636             else {
 637                 struct section *s;
 638                 Matrix *M;
 639                 int nd = 0;
 640                 int k, l, k2, l2, q;
 641                 evalue *L, *U;
 642                 evalue F;
 643                 /* put those with positive coefficients first; number: p */
 644                 for (p = 0, n = singles[i][0]-1; p <= n; ) {
 645                     while (value_pos_p(P->Constraint[singles[i][1+p]][i+1]))
 646                         ++p;
 647                     while (value_neg_p(P->Constraint[singles[i][1+n]][i+1]))
 648                         --n;
 649                     if (p < n) {
 650                         int t = singles[i][1+p];
 651                         singles[i][1+p] = singles[i][1+n];
 652                         singles[i][1+n] = t;
 653                         ++p;
 654                         --n;
 655                     }
 656                 }
 657                 n = singles[i][0]-p;
 658                 assert (p >= 1 && n >= 1);
 659                 s = (struct section *) malloc(p * n * sizeof(struct section));
 660                 M = Matrix_Alloc((p-1) + (n-1), dim-nvar+2);
 661                 for (k = 0; k < p; ++k) {
 662                     for (k2 = 0; k2 < p; ++k2) {
 663                         if (k2 == k)
 664                             continue;
 665                         q = k2 - (k2 > k);
 666                         value_oppose(tmp, P->Constraint[singles[i][1+k2]][i+1]);
 667                         value_set_si(M->p[q][0], 1);
 668                         Vector_Combine(P->Constraint[singles[i][1+k]]+1+nvar,
 669                                        P->Constraint[singles[i][1+k2]]+1+nvar,
 670                                        M->p[q]+1,
 671                                        tmp,
 672                                        P->Constraint[singles[i][1+k]][i+1],
 673                                        dim-nvar+1);
 674                     }
 675                     for (l = p; l < p+n; ++l) {
 676                         value_oppose(tmp, P->Constraint[singles[i][1+l]][i+1]);
 677                         for (l2 = p; l2 < p+n; ++l2) {
 678                             if (l2 == l)
 679                                 continue;
 680                             q = l2-1 - (l2 > l);
 681                             value_set_si(M->p[q][0], 1);
 682                             Vector_Combine(P->Constraint[singles[i][1+l2]]+1+nvar,
 683                                            P->Constraint[singles[i][1+l]]+1+nvar,
 684                                            M->p[q]+1,
 685                                            tmp,
 686                                            P->Constraint[singles[i][1+l2]][i+1],
 687                                            dim-nvar+1);
 688                         }
 689                         s[nd].D = Constraints2Polyhedron(M, P->NbRays);
 690                         if (emptyQ(s[nd].D)) {
 691                             Polyhedron_Free(s[nd].D);
 692                             continue;
 693                         }
 694                         L = bv_ceil3(P->Constraint[singles[i][1+k]]+1+nvar,
 695                                   dim-nvar+1,
 696                                   P->Constraint[singles[i][1+k]][i+1]);
 697                         U = bv_ceil3(P->Constraint[singles[i][1+l]]+1+nvar,
 698                                   dim-nvar+1,
 699                                   P->Constraint[singles[i][1+l]][i+1]);
 700                         eadd(L, U);
 701                         eadd(&mone, U);
 702                         emul(&mone, U);
 703                         s[nd].E = *U;
 704                         free_evalue_refs(L);
 705                         free(L);
 706                         free(U);
 707                         ++nd;
 708                     }
 709                 }
 710
 711                 Matrix_Free(M);
 712
 713                 value_init(F.d);
 714                 value_set_si(F.d, 0);
 715                 F.x.p = new_enode(partition, 2*nd, -1);
 716                 for (k = 0; k < nd; ++k) {
 717                     EVALUE_SET_DOMAIN(F.x.p->arr[2*k], s[k].D);
 718                     value_clear(F.x.p->arr[2*k+1].d);
 719                     F.x.p->arr[2*k+1] = s[k].E;
 720                 }
 721                 free(s);
 722
 723                 emul(&F, factor);
 724                 free_evalue_refs(&F);
 725             }
 726         }
 727         value_set_si(m->p[dim-nsingle][dim], 1);
 728         P = Polyhedron_Image(P, m, P->NbConstraints);
 729         Matrix_Free(m);
 730         free_singles(singles, nvar);
 731
 732         value_clear(tmp);
 733
 734         free_evalue_refs(&mone);
 735     }
 736
 737     reduce_evalue(factor);
 738
 739     return P;
 740 }
 741
 742 #ifdef USE_MODULO
 743 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
 744                                    evalue* factor)
 745 {
 746     return ParamPolyhedron_Reduce_mod(P, nvar, factor);
 747 }
 748 #else
 749 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
 750                                    evalue* factor)
 751 {
 752     Polyhedron *R = ParamPolyhedron_Reduce_mod(P, nvar, factor);
 753     evalue_mod2table(factor, P->Dimension - nvar);
 754     reduce_evalue(factor);
 755     return R;
 756 }
 757 #endif
 758
 759 Bool isIdentity(Matrix *M)
 760 {
 761     unsigned i, j;
 762     if (M->NbRows != M->NbColumns)
 763         return False;
 764
 765     for (i = 0;i < M->NbRows; i ++)
 766         for (j = 0; j < M->NbColumns; j ++)
 767             if (i == j) {
 768                 if(value_notone_p(M->p[i][j]))
 769                     return False;
 770             } else {
 771                 if(value_notzero_p(M->p[i][j]))
 772                     return False;
 773             }
 774     return True;
 775 }
 776
 777 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
 778 {
 779   Param_Domain *P;
 780   Param_Vertices *V;
 781
 782   for(P=PP->D;P;P=P->next) {
 783
 784     /* prints current val. dom. */
 785     printf( "---------------------------------------\n" );
 786     printf( "Domain :\n");
 787     Print_Domain( stdout, P->Domain, param_names );
 788
 789     /* scan the vertices */
 790     printf( "Vertices :\n");
 791     FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
 792
 793       /* prints each vertex */
 794       Print_Vertex( stdout, V->Vertex, param_names );
 795       printf( "\n" );
 796     }
 797     END_FORALL_PVertex_in_ParamPolyhedron;
 798   }
 799 }
 800
 801 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
 802 {
 803     for (; en; en = en->next) {
 804         Print_Domain(Dst, en->ValidityDomain, params);
 805         print_evalue(Dst, &en->EP, params);
 806     }
 807 }
 808
 809 void Enumeration_mod2table(Enumeration *en, unsigned nparam)
 810 {
 811     for (; en; en = en->next)
 812         evalue_mod2table(&en->EP, nparam);
 813 }
 814
 815 void Free_ParamNames(char **params, int m)
 816 {
 817     while (--m >= 0)
 818         free(params[m]);
 819     free(params);
 820 }
 821
 822 int DomainIncludes(Polyhedron *Pol1, Polyhedron *Pol2)
 823 {
 824     Polyhedron *P2;
 825     for ( ; Pol1; Pol1 = Pol1->next) {
 826         for (P2 = Pol2; P2; P2 = P2->next)
 827             if (!PolyhedronIncludes(Pol1, P2))
 828                 break;
 829         if (!P2)
 830             return 1;
 831     }
 832     return 0;
 833 }