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, *M3;
 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     M3 = Matrix_Copy(M);
 223     L = Rays2Polyhedron(M3, NbMaxCons);
 224     Matrix_Free(M3);
 225
 226     M2 = Matrix_Alloc(dim+1, dim+2);
 227
 228     t = 1;
 229     if (0) {
 230 try_again:
 231         /* Usually R should still be 0 */
 232         Domain_Free(R);
 233         Polyhedron_Free(L);
 234         for (r = 1; r < P->NbRays; ++r) {
 235             value_set_si(M->p[r][dim+1], random_int((t+1)*dim)+1);
 236         }
 237         M3 = Matrix_Copy(M);
 238         L = Rays2Polyhedron(M3, NbMaxCons);
 239         Matrix_Free(M3);
 240         ++t;
 241     }
 242     assert(t <= MAX_TRY);
 243
 244     R = NULL;
 245     n = 0;
 246
 247     for (i = 0; i < L->NbConstraints; ++i) {
 248         if (value_negz_p(L->Constraint[i][dim+1]))
 249             continue;
 250         if (value_notzero_p(L->Constraint[i][dim+2]))
 251             continue;
 252         for (j = 1, r = 1; j < M->NbRows; ++j) {
 253             Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
 254             if (value_notzero_p(tmp))
 255                 continue;
 256             if (r > dim)
 257                 goto try_again;
 258             Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
 259             value_set_si(M2->p[r][0], 1);
 260             value_set_si(M2->p[r][dim+1], 0);
 261             ++r;
 262         }
 263         assert(r == dim+1);
 264         Vector_Set(M2->p[0]+1, 0, dim);
 265         value_set_si(M2->p[0][0], 1);
 266         value_set_si(M2->p[0][dim+1], 1);
 267         T = Rays2Polyhedron(M2, P->NbConstraints+1);
 268         T->next = R;
 269         R = T;
 270         ++n;
 271     }
 272     Matrix_Free(M2);
 273
 274     Polyhedron_Free(L);
 275     value_clear(tmp);
 276     Matrix_Free(M);
 277
 278     return R;
 279 }
 280
 281 void check_triangulization(Polyhedron *P, Polyhedron *T)
 282 {
 283     Polyhedron *C, *D, *E, *F, *G, *U;
 284     for (C = T; C; C = C->next) {
 285         if (C == T)
 286             U = C;
 287         else
 288             U = DomainConvex(DomainUnion(U, C, 100), 100);
 289         for (D = C->next; D; D = D->next) {
 290             F = C->next;
 291             G = D->next;
 292             C->next = NULL;
 293             D->next = NULL;
 294             E = DomainIntersection(C, D, 600);
 295             assert(E->NbRays == 0 || E->NbEq >= 1);
 296             Polyhedron_Free(E);
 297             C->next = F;
 298             D->next = G;
 299         }
 300     }
 301     assert(PolyhedronIncludes(U, P));
 302     assert(PolyhedronIncludes(P, U));
 303 }
 304
 305 void Euclid(Value a, Value b, Value *x, Value *y, Value *g)
 306 {
 307     Value c, d, e, f, tmp;
 308
 309     value_init(c);
 310     value_init(d);
 311     value_init(e);
 312     value_init(f);
 313     value_init(tmp);
 314     value_absolute(c, a);
 315     value_absolute(d, b);
 316     value_set_si(e, 1);
 317     value_set_si(f, 0);
 318     while(value_pos_p(d)) {
 319         value_division(tmp, c, d);
 320         value_multiply(tmp, tmp, f);
 321         value_substract(e, e, tmp);
 322         value_division(tmp, c, d);
 323         value_multiply(tmp, tmp, d);
 324         value_substract(c, c, tmp);
 325         value_swap(c, d);
 326         value_swap(e, f);
 327     }
 328     value_assign(*g, c);
 329     if (value_zero_p(a))
 330         value_set_si(*x, 0);
 331     else if (value_pos_p(a))
 332         value_assign(*x, e);
 333     else value_oppose(*x, e);
 334     if (value_zero_p(b))
 335         value_set_si(*y, 0);
 336     else {
 337         value_multiply(tmp, a, *x);
 338         value_substract(tmp, c, tmp);
 339         value_division(*y, tmp, b);
 340     }
 341     value_clear(c);
 342     value_clear(d);
 343     value_clear(e);
 344     value_clear(f);
 345     value_clear(tmp);
 346 }
 347
 348 Matrix * unimodular_complete(Vector *row)
 349 {
 350     Value g, b, c, old, tmp;
 351     Matrix *m;
 352     unsigned i, j;
 353
 354     value_init(b);
 355     value_init(c);
 356     value_init(g);
 357     value_init(old);
 358     value_init(tmp);
 359     m = Matrix_Alloc(row->Size, row->Size);
 360     for (j = 0; j < row->Size; ++j) {
 361         value_assign(m->p[0][j], row->p[j]);
 362     }
 363     value_assign(g, row->p[0]);
 364     for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
 365         for (j = 0; j < row->Size; ++j) {
 366             if (j == i-1)
 367                 value_set_si(m->p[i][j], 1);
 368             else
 369                 value_set_si(m->p[i][j], 0);
 370         }
 371         value_assign(g, row->p[i]);
 372     }
 373     for (; i < row->Size; ++i) {
 374         value_assign(old, g);
 375         Euclid(old, row->p[i], &c, &b, &g);
 376         value_oppose(b, b);
 377         for (j = 0; j < row->Size; ++j) {
 378             if (j < i) {
 379                 value_multiply(tmp, row->p[j], b);
 380                 value_division(m->p[i][j], tmp, old);
 381             } else if (j == i)
 382                 value_assign(m->p[i][j], c);
 383             else
 384                 value_set_si(m->p[i][j], 0);
 385         }
 386     }
 387     value_clear(b);
 388     value_clear(c);
 389     value_clear(g);
 390     value_clear(old);
 391     value_clear(tmp);
 392     return m;
 393 }
 394
 395 /*
 396  * Returns a full-dimensional polyhedron with the same number
 397  * of integer points as P
 398  */
 399 Polyhedron *remove_equalities(Polyhedron *P)
 400 {
 401     Value g;
 402     Vector *v;
 403     Polyhedron *p = Polyhedron_Copy(P), *q;
 404     unsigned dim = p->Dimension;
 405     Matrix *m1, *m2;
 406     int i;
 407
 408     value_init(g);
 409     while (p->NbEq > 0) {
 410         assert(dim > 0);
 411         Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
 412         Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
 413         Vector_Gcd(p->Constraint[0]+1, dim, &g);
 414         if (value_notone_p(g) && value_notmone_p(g)) {
 415             Polyhedron_Free(p);
 416             p = Empty_Polyhedron(0);
 417             break;
 418         }
 419         v = Vector_Alloc(dim);
 420         Vector_Copy(p->Constraint[0]+1, v->p, dim);
 421         m1 = unimodular_complete(v);
 422         m2 = Matrix_Alloc(dim, dim+1);
 423         for (i = 0; i < dim-1 ; ++i) {
 424             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 425             value_set_si(m2->p[i][dim], 0);
 426         }
 427         Vector_Set(m2->p[dim-1], 0, dim);
 428         value_set_si(m2->p[dim-1][dim], 1);
 429         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 430         Vector_Free(v);
 431         Matrix_Free(m1);
 432         Matrix_Free(m2);
 433         Polyhedron_Free(p);
 434         p = q;
 435         --dim;
 436     }
 437     value_clear(g);
 438     return p;
 439 }
 440
 441 /*
 442  * Returns a full-dimensional polyhedron with the same number
 443  * of integer points as P
 444  * nvar specifies the number of variables
 445  * The remaining dimensions are assumed to be parameters
 446  * Destroys P
 447  * factor is NbEq x (nparam+2) matrix, containing stride constraints
 448  * on the parameters; column nparam is the constant;
 449  * column nparam+1 is the stride
 450  */
 451 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
 452 {
 453     Value g;
 454     Vector *v;
 455     Polyhedron *p = P, *q;
 456     unsigned dim = p->Dimension;
 457     Matrix *m1, *m2, *f;
 458     int i, j, skip;
 459
 460     value_init(g);
 461     f = Matrix_Alloc(p->NbEq, dim-nvar+2);
 462     j = 0;
 463     *factor = f;
 464     skip = 0;
 465     while (nvar > 0 && p->NbEq - skip > 0) {
 466         assert(dim > 0);
 467
 468         while (value_zero_p(p->Constraint[skip][0]) &&
 469                First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
 470             ++skip;
 471         if (p->NbEq == skip)
 472             break;
 473
 474         Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
 475         Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
 476         Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
 477         Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
 478         value_assign(f->p[j][dim-nvar+1], g);
 479         v = Vector_Alloc(dim);
 480         Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
 481         Vector_Set(v->p+nvar, 0, dim-nvar);
 482         m1 = unimodular_complete(v);
 483         m2 = Matrix_Alloc(dim, dim+1);
 484         for (i = 0; i < dim-1 ; ++i) {
 485             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 486             value_set_si(m2->p[i][dim], 0);
 487         }
 488         Vector_Set(m2->p[dim-1], 0, dim);
 489         value_set_si(m2->p[dim-1][dim], 1);
 490         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 491         Vector_Free(v);
 492         Matrix_Free(m1);
 493         Matrix_Free(m2);
 494         Polyhedron_Free(p);
 495         p = q;
 496         --dim;
 497         --nvar;
 498         ++j;
 499     }
 500     value_clear(g);
 501     return p;
 502 }
 503
 504 struct single {
 505     int nr;
 506     int pos[2];
 507 };
 508
 509 static void free_singles(int **singles, int dim)
 510 {
 511     int i;
 512     for (i = 0; i < dim; ++i)
 513         free(singles[i]);
 514     free(singles);
 515 }
 516
 517 static int **find_singles(Polyhedron *P, int dim, int max, int *nsingle)
 518 {
 519     int i, j, prev;
 520     int **singles = (int **) malloc(dim * sizeof(int *));
 521     assert(singles);
 522
 523     for (i = 0; i < dim; ++i) {
 524         singles[i] = (int *) malloc((max + 1) *sizeof(int));
 525         singles[i][0] = 0;
 526     }
 527
 528     for (i = 0; i < P->NbConstraints; ++i) {
 529         for (j = 0, prev = -1; j < dim; ++j) {
 530             if (value_notzero_p(P->Constraint[i][j+1])) {
 531                 if (prev == -1)
 532                     prev = j;
 533                 else {
 534                     if (prev != -2)
 535                         singles[prev][0] = -1;
 536                     singles[j][0] = -1;
 537                     prev = -2;
 538                 }
 539             }
 540         }
 541         if (prev >= 0 && singles[prev][0] >= 0)
 542             singles[prev][++singles[prev][0]] = i;
 543     }
 544     *nsingle = 0;
 545     for (j = 0; j < dim; ++j)
 546         if (singles[j][0] > 0)
 547             ++*nsingle;
 548     if (!*nsingle) {
 549         free_singles(singles, dim);
 550         singles = 0;
 551     }
 552     return singles;
 553 }
 554
 555 /*
 556  * The number of points in P is equal to factor time
 557  * the number of points in the polyhedron returned.
 558  * The return value is zero if no reduction can be found.
 559  */
 560 Polyhedron* Polyhedron_Reduce(Polyhedron *P, Value* factor)
 561 {
 562     int i, j, nsingle, k, p;
 563     unsigned dim = P->Dimension;
 564     int **singles;
 565
 566     value_set_si(*factor, 1);
 567
 568     assert (P->NbEq == 0);
 569
 570     singles = find_singles(P, dim, 2, &nsingle);
 571
 572     if (nsingle == 0)
 573         return 0;
 574
 575     {
 576         Value tmp, pos, neg;
 577         Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
 578
 579         value_init(tmp);
 580         value_init(pos);
 581         value_init(neg);
 582
 583         for (i = 0, j = 0; i < dim; ++i) {
 584             if (singles[i][0] != 2)
 585                 value_set_si(m->p[j++][i], 1);
 586             else {
 587                 for (k = 0; k <= 1; ++k) {
 588                     p = singles[i][1+k];
 589                     value_oppose(tmp, P->Constraint[p][dim+1]);
 590                     if (value_pos_p(P->Constraint[p][i+1]))
 591                         mpz_cdiv_q(pos, tmp, P->Constraint[p][i+1]);
 592                     else
 593                         mpz_fdiv_q(neg, tmp, P->Constraint[p][i+1]);
 594                 }
 595                 value_substract(tmp, neg, pos);
 596                 value_increment(tmp, tmp);
 597                 value_multiply(*factor, *factor, tmp);
 598             }
 599         }
 600         value_set_si(m->p[dim-nsingle][dim], 1);
 601         P = Polyhedron_Image(P, m, P->NbConstraints);
 602         Matrix_Free(m);
 603         free_singles(singles, dim);
 604
 605         value_clear(tmp);
 606         value_clear(pos);
 607         value_clear(neg);
 608     }
 609
 610     return P;
 611 }
 612
 613 /*
 614  * Replaces constraint a x >= c by x >= ceil(c/a)
 615  * where "a" is a common factor in the coefficients
 616  * old is the constraint; v points to an initialized
 617  * value that this procedure can use.
 618  * Return non-zero if something changed.
 619  * Result is placed in new.
 620  */
 621 int ConstraintSimplify(Value *old, Value *new, int len, Value* v)
 622 {
 623     Vector_Gcd(old+1, len - 2, v);
 624
 625     if (value_one_p(*v))
 626         return 0;
 627
 628     Vector_AntiScale(old+1, new+1, *v, len-2);
 629     mpz_fdiv_q(new[len-1], old[len-1], *v);
 630     return 1;
 631 }
 632
 633 /*
 634  * Project on final dim dimensions
 635  */
 636 static Polyhedron* Polyhedron_Project(Polyhedron *P, int dim)
 637 {
 638     int i;
 639     int remove = P->Dimension - dim;
 640
 641     if (P->Dimension == dim)
 642         return Polyhedron_Copy(P);
 643
 644     Matrix *T = Matrix_Alloc(dim+1, P->Dimension+1);
 645     for (i = 0; i < dim+1; ++i)
 646         value_set_si(T->p[i][i+remove], 1);
 647     Polyhedron *I = Polyhedron_Image(P, T, P->NbConstraints);
 648     Matrix_Free(T);
 649     return I;
 650 }
 651
 652 struct section { Polyhedron * D; evalue E; };
 653
 654 static Polyhedron* ParamPolyhedron_Reduce_mod(Polyhedron *P, unsigned nvar,
 655                                               evalue* factor)
 656 {
 657     int nsingle;
 658     int **singles;
 659     unsigned dim = P->Dimension;
 660
 661     singles = find_singles(P, nvar, P->NbConstraints, &nsingle);
 662
 663     if (nsingle == 0)
 664         return 0;
 665
 666     {
 667         Polyhedron *C, *T;
 668         int i, j, p, n;
 669         Matrix *m = Matrix_Alloc((dim-nsingle)+1, dim+1);
 670         Value tmp, g;
 671         evalue mone;
 672         value_init(mone.d);
 673         evalue_set_si(&mone, -1, 1);
 674         C = Polyhedron_Project(P, dim-nvar);
 675
 676         value_init(tmp);
 677         value_init(g);
 678
 679         for (i = 0, j = 0; i < dim; ++i) {
 680             if (i >= nvar || singles[i][0] < 2)
 681                 value_set_si(m->p[j++][i], 1);
 682             else {
 683                 struct section *s;
 684                 Matrix *M, *M2;
 685                 int nd = 0;
 686                 int k, l, k2, l2, q;
 687                 evalue *L, *U;
 688                 evalue F;
 689                 /* put those with positive coefficients first; number: p */
 690                 for (p = 0, n = singles[i][0]-1; p <= n; ) {
 691                     while (value_pos_p(P->Constraint[singles[i][1+p]][i+1]))
 692                         ++p;
 693                     while (value_neg_p(P->Constraint[singles[i][1+n]][i+1]))
 694                         --n;
 695                     if (p < n) {
 696                         int t = singles[i][1+p];
 697                         singles[i][1+p] = singles[i][1+n];
 698                         singles[i][1+n] = t;
 699                         ++p;
 700                         --n;
 701                     }
 702                 }
 703                 n = singles[i][0]-p;
 704                 assert (p >= 1 && n >= 1);
 705                 s = (struct section *) malloc(p * n * sizeof(struct section));
 706                 M = Matrix_Alloc((p-1) + (n-1), dim-nvar+2);
 707                 for (k = 0; k < p; ++k) {
 708                     for (k2 = 0; k2 < p; ++k2) {
 709                         if (k2 == k)
 710                             continue;
 711                         q = k2 - (k2 > k);
 712                         value_oppose(tmp, P->Constraint[singles[i][1+k2]][i+1]);
 713                         value_set_si(M->p[q][0], 1);
 714                         Vector_Combine(P->Constraint[singles[i][1+k]]+1+nvar,
 715                                        P->Constraint[singles[i][1+k2]]+1+nvar,
 716                                        M->p[q]+1,
 717                                        tmp,
 718                                        P->Constraint[singles[i][1+k]][i+1],
 719                                        dim-nvar+1);
 720                         if (k2 > k)
 721                             value_decrement(M->p[q][dim-nvar+1],
 722                                             M->p[q][dim-nvar+1]);
 723                         ConstraintSimplify(M->p[q], M->p[q],
 724                                            dim-nvar+2, &g);
 725                     }
 726                     for (l = p; l < p+n; ++l) {
 727                         value_oppose(tmp, P->Constraint[singles[i][1+l]][i+1]);
 728                         for (l2 = p; l2 < p+n; ++l2) {
 729                             if (l2 == l)
 730                                 continue;
 731                             q = l2-1 - (l2 > l);
 732                             value_set_si(M->p[q][0], 1);
 733                             Vector_Combine(P->Constraint[singles[i][1+l2]]+1+nvar,
 734                                            P->Constraint[singles[i][1+l]]+1+nvar,
 735                                            M->p[q]+1,
 736                                            tmp,
 737                                            P->Constraint[singles[i][1+l2]][i+1],
 738                                            dim-nvar+1);
 739                             if (l2 > l)
 740                                 value_decrement(M->p[q][dim-nvar+1],
 741                                                 M->p[q][dim-nvar+1]);
 742                             ConstraintSimplify(M->p[q], M->p[q],
 743                                                dim-nvar+2, &g);
 744                         }
 745                         M2 = Matrix_Copy(M);
 746                         s[nd].D = Constraints2Polyhedron(M2, P->NbRays);
 747                         Matrix_Free(M2);
 748                         if (emptyQ(s[nd].D)) {
 749                             Polyhedron_Free(s[nd].D);
 750                             continue;
 751                         }
 752                         T = DomainIntersection(s[nd].D, C, 0);
 753                         L = bv_ceil3(P->Constraint[singles[i][1+k]]+1+nvar,
 754                                   dim-nvar+1,
 755                                   P->Constraint[singles[i][1+k]][i+1], T);
 756                         U = bv_ceil3(P->Constraint[singles[i][1+l]]+1+nvar,
 757                                   dim-nvar+1,
 758                                   P->Constraint[singles[i][1+l]][i+1], T);
 759                         Domain_Free(T);
 760                         eadd(L, U);
 761                         eadd(&mone, U);
 762                         emul(&mone, U);
 763                         s[nd].E = *U;
 764                         free_evalue_refs(L);
 765                         free(L);
 766                         free(U);
 767                         ++nd;
 768                     }
 769                 }
 770
 771                 Matrix_Free(M);
 772
 773                 value_init(F.d);
 774                 value_set_si(F.d, 0);
 775                 F.x.p = new_enode(partition, 2*nd, -1);
 776                 for (k = 0; k < nd; ++k) {
 777                     EVALUE_SET_DOMAIN(F.x.p->arr[2*k], s[k].D);
 778                     value_clear(F.x.p->arr[2*k+1].d);
 779                     F.x.p->arr[2*k+1] = s[k].E;
 780                 }
 781                 free(s);
 782
 783                 emul(&F, factor);
 784                 free_evalue_refs(&F);
 785             }
 786         }
 787         value_set_si(m->p[dim-nsingle][dim], 1);
 788         P = Polyhedron_Image(P, m, P->NbConstraints);
 789         Matrix_Free(m);
 790         free_singles(singles, nvar);
 791
 792         value_clear(g);
 793         value_clear(tmp);
 794
 795         free_evalue_refs(&mone);
 796         Polyhedron_Free(C);
 797     }
 798
 799     reduce_evalue(factor);
 800
 801     return P;
 802 }
 803
 804 #ifdef USE_MODULO
 805 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
 806                                    evalue* factor)
 807 {
 808     return ParamPolyhedron_Reduce_mod(P, nvar, factor);
 809 }
 810 #else
 811 Polyhedron* ParamPolyhedron_Reduce(Polyhedron *P, unsigned nvar,
 812                                    evalue* factor)
 813 {
 814     Polyhedron *R;
 815     evalue tmp;
 816     value_init(tmp.d);
 817     evalue_set_si(&tmp, 1, 1);
 818     R = ParamPolyhedron_Reduce_mod(P, nvar, &tmp);
 819     evalue_mod2table(&tmp, P->Dimension - nvar);
 820     reduce_evalue(&tmp);
 821     emul(&tmp, factor);
 822     free_evalue_refs(&tmp);
 823     return R;
 824 }
 825 #endif
 826
 827 Bool isIdentity(Matrix *M)
 828 {
 829     unsigned i, j;
 830     if (M->NbRows != M->NbColumns)
 831         return False;
 832
 833     for (i = 0;i < M->NbRows; i ++)
 834         for (j = 0; j < M->NbColumns; j ++)
 835             if (i == j) {
 836                 if(value_notone_p(M->p[i][j]))
 837                     return False;
 838             } else {
 839                 if(value_notzero_p(M->p[i][j]))
 840                     return False;
 841             }
 842     return True;
 843 }
 844
 845 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
 846 {
 847   Param_Domain *P;
 848   Param_Vertices *V;
 849
 850   for(P=PP->D;P;P=P->next) {
 851
 852     /* prints current val. dom. */
 853     printf( "---------------------------------------\n" );
 854     printf( "Domain :\n");
 855     Print_Domain( stdout, P->Domain, param_names );
 856
 857     /* scan the vertices */
 858     printf( "Vertices :\n");
 859     FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
 860
 861       /* prints each vertex */
 862       Print_Vertex( stdout, V->Vertex, param_names );
 863       printf( "\n" );
 864     }
 865     END_FORALL_PVertex_in_ParamPolyhedron;
 866   }
 867 }
 868
 869 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
 870 {
 871     for (; en; en = en->next) {
 872         Print_Domain(Dst, en->ValidityDomain, params);
 873         print_evalue(Dst, &en->EP, params);
 874     }
 875 }
 876
 877 void Enumeration_mod2table(Enumeration *en, unsigned nparam)
 878 {
 879     for (; en; en = en->next) {
 880         evalue_mod2table(&en->EP, nparam);
 881         reduce_evalue(&en->EP);
 882     }
 883 }
 884
 885 size_t Enumeration_size(Enumeration *en)
 886 {
 887     size_t s = 0;
 888
 889     for (; en; en = en->next) {
 890         s += domain_size(en->ValidityDomain);
 891         s += evalue_size(&en->EP);
 892     }
 893     return s;
 894 }
 895
 896 void Free_ParamNames(char **params, int m)
 897 {
 898     while (--m >= 0)
 899         free(params[m]);
 900     free(params);
 901 }
 902
 903 int DomainIncludes(Polyhedron *Pol1, Polyhedron *Pol2)
 904 {
 905     Polyhedron *P2;
 906     for ( ; Pol1; Pol1 = Pol1->next) {
 907         for (P2 = Pol2; P2; P2 = P2->next)
 908             if (!PolyhedronIncludes(Pol1, P2))
 909                 break;
 910         if (!P2)
 911             return 1;
 912     }
 913     return 0;
 914 }
 915
 916 static Polyhedron *p_simplify_constraints(Polyhedron *P, Vector *row,
 917                                           Value *g, unsigned MaxRays)
 918 {
 919     Polyhedron *T, *R = P;
 920     int len = P->Dimension+2;
 921     int r;
 922
 923     /* Also look at equalities.
 924      * If an equality can be "simplified" then there
 925      * are no integer solutions anyway and the following loop
 926      * will add a conflicting constraint
 927      */
 928     for (r = 0; r < R->NbConstraints; ++r) {
 929         if (ConstraintSimplify(R->Constraint[r], row->p, len, g)) {
 930             T = R;
 931             R = AddConstraints(row->p, 1, R, MaxRays);
 932             if (T != P)
 933                 Polyhedron_Free(T);
 934             r = -1;
 935         }
 936     }
 937     if (R != P)
 938         Polyhedron_Free(P);
 939     return R;
 940 }
 941
 942 /*
 943  * Replaces constraint a x >= c by x >= ceil(c/a)
 944  * where "a" is a common factor in the coefficients
 945  * Destroys P and returns a newly allocated Polyhedron
 946  * or just returns P in case no changes were made
 947  */
 948 Polyhedron *DomainConstraintSimplify(Polyhedron *P, unsigned MaxRays)
 949 {
 950     Polyhedron **prev;
 951     int len = P->Dimension+2;
 952     Vector *row = Vector_Alloc(len);
 953     value_set_si(row->p[0], 1);
 954     Value g;
 955     value_init(g);
 956     Polyhedron *R = P, *N;
 957
 958     for (prev = &R; P; P = N) {
 959         Polyhedron *T;
 960         N = P->next;
 961         T = p_simplify_constraints(P, row, &g, MaxRays);
 962
 963         if (emptyQ(T) && prev != &R) {
 964             Polyhedron_Free(T);
 965             *prev = NULL;
 966             continue;
 967         }
 968
 969         if (T != P)
 970             T->next = N;
 971         *prev = T;
 972         prev = &T->next;
 973     }
 974
 975     if (R->next && emptyQ(R)) {
 976         N = R->next;
 977         Polyhedron_Free(R);
 978         R = N;
 979     }
 980
 981     value_clear(g);
 982     Vector_Free(row);
 983     return R;
 984 }
 985
 986 int line_minmax(Polyhedron *I, Value *min, Value *max)
 987 {
 988     int i;
 989
 990     if (I->NbEq >= 1) {
 991         value_oppose(I->Constraint[0][2], I->Constraint[0][2]);
 992         /* There should never be a remainder here */
 993         if (value_pos_p(I->Constraint[0][1]))
 994             mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
 995         else
 996             mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
 997         value_assign(*max, *min);
 998     } else for (i = 0; i < I->NbConstraints; ++i) {
 999         if (value_zero_p(I->Constraint[i][1])) {
1000             Polyhedron_Free(I);
1001             return 0;
1002         }
1003
1004         value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
1005         if (value_pos_p(I->Constraint[i][1]))
1006             mpz_cdiv_q(*min, I->Constraint[i][2], I->Constraint[i][1]);
1007         else
1008             mpz_fdiv_q(*max, I->Constraint[i][2], I->Constraint[i][1]);
1009     }
1010     Polyhedron_Free(I);
1011     return 1;
1012 }
1013
1014 /**
1015
1016 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1017 each imbriquation
1018
1019 @param pos index position of current loop index (1..hdim-1)
1020 @param P loop domain
1021 @param context context values for fixed indices
1022 @param exist    number of existential variables
1023 @return the number of integer points in this
1024 polyhedron
1025
1026 */
1027 void count_points_e (int pos, Polyhedron *P, int exist, int nparam,
1028                      Value *context, Value *res)
1029 {
1030     Value LB, UB, k, c;
1031
1032     if (emptyQ(P)) {
1033         value_set_si(*res, 0);
1034         return;
1035     }
1036
1037     value_init(LB); value_init(UB); value_init(k);
1038     value_set_si(LB,0);
1039     value_set_si(UB,0);
1040
1041     if (lower_upper_bounds(pos,P,context,&LB,&UB) !=0) {
1042         /* Problem if UB or LB is INFINITY */
1043         value_clear(LB); value_clear(UB); value_clear(k);
1044         if (pos > P->Dimension - nparam - exist)
1045             value_set_si(*res, 1);
1046         else
1047             value_set_si(*res, -1);
1048         return;
1049     }
1050
1051 #ifdef EDEBUG1
1052     if (!P->next) {
1053         int i;
1054         for (value_assign(k,LB); value_le(k,UB); value_increment(k,k)) {
1055             fprintf(stderr, "(");
1056             for (i=1; i<pos; i++) {
1057                 value_print(stderr,P_VALUE_FMT,context[i]);
1058                 fprintf(stderr,",");
1059             }
1060             value_print(stderr,P_VALUE_FMT,k);
1061             fprintf(stderr,")\n");
1062         }
1063     }
1064 #endif
1065
1066     value_set_si(context[pos],0);
1067     if (value_lt(UB,LB)) {
1068         value_clear(LB); value_clear(UB); value_clear(k);
1069         value_set_si(*res, 0);
1070         return;
1071     }
1072     if (!P->next) {
1073         if (exist)
1074             value_set_si(*res, 1);
1075         else {
1076             value_substract(k,UB,LB);
1077             value_add_int(k,k,1);
1078             value_assign(*res, k);
1079         }
1080         value_clear(LB); value_clear(UB); value_clear(k);
1081         return;
1082     }
1083
1084     /*-----------------------------------------------------------------*/
1085     /* Optimization idea                                               */
1086     /*   If inner loops are not a function of k (the current index)    */
1087     /*   i.e. if P->Constraint[i][pos]==0 for all P following this and */
1088     /*        for all i,                                               */
1089     /*   Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context)    */
1090     /*   (skip the for loop)                                           */
1091     /*-----------------------------------------------------------------*/
1092
1093     value_init(c);
1094     value_set_si(*res, 0);
1095     for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
1096         /* Insert k in context */
1097         value_assign(context[pos],k);
1098         count_points_e(pos+1, P->next, exist, nparam, context, &c);
1099         if(value_notmone_p(c))
1100             value_addto(*res, *res, c);
1101         else {
1102             value_set_si(*res, -1);
1103             break;
1104         }
1105         if (pos > P->Dimension - nparam - exist &&
1106                 value_pos_p(*res))
1107             break;
1108     }
1109     value_clear(c);
1110
1111 #ifdef EDEBUG11
1112     fprintf(stderr,"%d\n",CNT);
1113 #endif
1114
1115     /* Reset context */
1116     value_set_si(context[pos],0);
1117     value_clear(LB); value_clear(UB); value_clear(k);
1118     return;
1119 } /* count_points_e */