util.c

   1 #include <barvinok/util.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 #define ALLOC(type) (type*)malloc(sizeof(type))
  11 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  12
  13 #ifdef __GNUC__
  14 #define NALLOC(p,n) p = (typeof(p))malloc((n) * sizeof(*p))
  15 #else
  16 #define NALLOC(p,n) p = (void *)malloc((n) * sizeof(*p))
  17 #endif
  18
  19 #ifndef HAVE_ENUMERATION_FREE
  20 #define Enumeration_Free(en)    /* just leak some memory */
  21 #endif
  22
  23 void manual_count(Polyhedron *P, Value* result)
  24 {
  25     Polyhedron *U = Universe_Polyhedron(0);
  26     Enumeration *en = Polyhedron_Enumerate(P,U,1024,NULL);
  27     Value *v = compute_poly(en,NULL);
  28     value_assign(*result, *v);
  29     value_clear(*v);
  30     free(v);
  31     Enumeration_Free(en);
  32     Polyhedron_Free(U);
  33 }
  34
  35 #ifndef HAVE_ENUMERATION_FREE
  36 #undef Enumeration_Free
  37 #endif
  38
  39 #include <barvinok/evalue.h>
  40 #include <barvinok/util.h>
  41 #include <barvinok/barvinok.h>
  42
  43 /* Return random value between 0 and max-1 inclusive
  44  */
  45 int random_int(int max) {
  46     return (int) (((double)(max))*rand()/(RAND_MAX+1.0));
  47 }
  48
  49 Polyhedron *Polyhedron_Read(unsigned MaxRays)
  50 {
  51     int vertices = 0;
  52     unsigned NbRows, NbColumns;
  53     Matrix *M;
  54     Polyhedron *P;
  55     char s[128];
  56
  57     while (fgets(s, sizeof(s), stdin)) {
  58         if (*s == '#')
  59             continue;
  60         if (strncasecmp(s, "vertices", sizeof("vertices")-1) == 0)
  61             vertices = 1;
  62         if (sscanf(s, "%u %u", &NbRows, &NbColumns) == 2)
  63             break;
  64     }
  65     if (feof(stdin))
  66         return NULL;
  67     M = Matrix_Alloc(NbRows,NbColumns);
  68     Matrix_Read_Input(M);
  69     if (vertices)
  70         P = Rays2Polyhedron(M, MaxRays);
  71     else
  72         P = Constraints2Polyhedron(M, MaxRays);
  73     Matrix_Free(M);
  74     return P;
  75 }
  76
  77 /* Inplace polarization
  78  */
  79 void Polyhedron_Polarize(Polyhedron *P)
  80 {
  81     unsigned NbRows = P->NbConstraints + P->NbRays;
  82     int i;
  83     Value **q;
  84
  85     q = (Value **)malloc(NbRows * sizeof(Value *));
  86     assert(q);
  87     for (i = 0; i < P->NbRays; ++i)
  88         q[i] = P->Ray[i];
  89     for (; i < NbRows; ++i)
  90         q[i] = P->Constraint[i-P->NbRays];
  91     P->NbConstraints = NbRows - P->NbConstraints;
  92     P->NbRays = NbRows - P->NbRays;
  93     free(P->Constraint);
  94     P->Constraint = q;
  95     P->Ray = q + P->NbConstraints;
  96 }
  97
  98 /*
  99  * Rather general polar
 100  * We can optimize it significantly if we assume that
 101  * P includes zero
 102  *
 103  * Also, we calculate the polar as defined in Schrijver
 104  * The opposite should probably work as well and would
 105  * eliminate the need for multiplying by -1
 106  */
 107 Polyhedron* Polyhedron_Polar(Polyhedron *P, unsigned NbMaxRays)
 108 {
 109     int i;
 110     Value mone;
 111     unsigned dim = P->Dimension + 2;
 112     Matrix *M = Matrix_Alloc(P->NbRays, dim);
 113
 114     assert(M);
 115     value_init(mone);
 116     value_set_si(mone, -1);
 117     for (i = 0; i < P->NbRays; ++i) {
 118         Vector_Scale(P->Ray[i], M->p[i], mone, dim);
 119         value_multiply(M->p[i][0], M->p[i][0], mone);
 120         value_multiply(M->p[i][dim-1], M->p[i][dim-1], mone);
 121     }
 122     P = Constraints2Polyhedron(M, NbMaxRays);
 123     assert(P);
 124     Matrix_Free(M);
 125     value_clear(mone);
 126     return P;
 127 }
 128
 129 /*
 130  * Returns the supporting cone of P at the vertex with index v
 131  */
 132 Polyhedron* supporting_cone(Polyhedron *P, int v)
 133 {
 134     Matrix *M;
 135     Value tmp;
 136     int i, n, j;
 137     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
 138     unsigned dim = P->Dimension + 2;
 139
 140     assert(v >=0 && v < P->NbRays);
 141     assert(value_pos_p(P->Ray[v][dim-1]));
 142     assert(supporting);
 143
 144     value_init(tmp);
 145     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
 146         Inner_Product(P->Constraint[i] + 1, P->Ray[v] + 1, dim - 1, &tmp);
 147         if ((supporting[i] = value_zero_p(tmp)))
 148             ++n;
 149     }
 150     assert(n >= dim - 2);
 151     value_clear(tmp);
 152     M = Matrix_Alloc(n, dim);
 153     assert(M);
 154     for (i = 0, j = 0; i < P->NbConstraints; ++i)
 155         if (supporting[i]) {
 156             value_set_si(M->p[j][dim-1], 0);
 157             Vector_Copy(P->Constraint[i], M->p[j++], dim-1);
 158         }
 159     free(supporting);
 160     P = Constraints2Polyhedron(M, P->NbRays+1);
 161     assert(P);
 162     Matrix_Free(M);
 163     return P;
 164 }
 165
 166 void value_lcm(Value i, Value j, Value* lcm)
 167 {
 168     Value aux;
 169     value_init(aux);
 170     value_multiply(aux,i,j);
 171     Gcd(i,j,lcm);
 172     value_division(*lcm,aux,*lcm);
 173     value_clear(aux);
 174 }
 175
 176 Polyhedron* supporting_cone_p(Polyhedron *P, Param_Vertices *v)
 177 {
 178     Matrix *M;
 179     Value lcm, tmp, tmp2;
 180     unsigned char *supporting = (unsigned char *)malloc(P->NbConstraints);
 181     unsigned dim = P->Dimension + 2;
 182     unsigned nparam = v->Vertex->NbColumns - 2;
 183     unsigned nvar = dim - nparam - 2;
 184     int i, n, j;
 185     Vector *row;
 186
 187     assert(supporting);
 188     row = Vector_Alloc(nparam+1);
 189     assert(row);
 190     value_init(lcm);
 191     value_init(tmp);
 192     value_init(tmp2);
 193     value_set_si(lcm, 1);
 194     for (i = 0, n = 0; i < P->NbConstraints; ++i) {
 195         Vector_Set(row->p, 0, nparam+1);
 196         for (j = 0 ; j < nvar; ++j) {
 197             value_set_si(tmp, 1);
 198             value_assign(tmp2,  P->Constraint[i][j+1]);
 199             if (value_ne(lcm, v->Vertex->p[j][nparam+1])) {
 200                 value_assign(tmp, lcm);
 201                 value_lcm(lcm, v->Vertex->p[j][nparam+1], &lcm);
 202                 value_division(tmp, lcm, tmp);
 203                 value_multiply(tmp2, tmp2, lcm);
 204                 value_division(tmp2, tmp2, v->Vertex->p[j][nparam+1]);
 205             }
 206             Vector_Combine(row->p, v->Vertex->p[j], row->p,
 207                            tmp, tmp2, nparam+1);
 208         }
 209         value_set_si(tmp, 1);
 210         Vector_Combine(row->p, P->Constraint[i]+1+nvar, row->p, tmp, lcm, nparam+1);
 211         for (j = 0; j < nparam+1; ++j)
 212             if (value_notzero_p(row->p[j]))
 213                 break;
 214         if ((supporting[i] = (j == nparam + 1)))
 215             ++n;
 216     }
 217     assert(n >= nvar);
 218     value_clear(tmp);
 219     value_clear(tmp2);
 220     value_clear(lcm);
 221     Vector_Free(row);
 222     M = Matrix_Alloc(n, nvar+2);
 223     assert(M);
 224     for (i = 0, j = 0; i < P->NbConstraints; ++i)
 225         if (supporting[i]) {
 226             value_set_si(M->p[j][nvar+1], 0);
 227             Vector_Copy(P->Constraint[i], M->p[j++], nvar+1);
 228         }
 229     free(supporting);
 230     P = Constraints2Polyhedron(M, P->NbRays+1);
 231     assert(P);
 232     Matrix_Free(M);
 233     return P;
 234 }
 235
 236 Polyhedron* triangulate_cone(Polyhedron *P, unsigned NbMaxCons)
 237 {
 238     const static int MAX_TRY=10;
 239     int i, j, r, n, t;
 240     Value tmp;
 241     unsigned dim = P->Dimension;
 242     Matrix *M = Matrix_Alloc(P->NbRays+1, dim+3);
 243     Matrix *M2, *M3;
 244     Polyhedron *L, *R, *T;
 245     assert(P->NbEq == 0);
 246
 247     R = NULL;
 248     value_init(tmp);
 249
 250     Vector_Set(M->p[0]+1, 0, dim+1);
 251     value_set_si(M->p[0][0], 1);
 252     value_set_si(M->p[0][dim+2], 1);
 253     Vector_Set(M->p[P->NbRays]+1, 0, dim+2);
 254     value_set_si(M->p[P->NbRays][0], 1);
 255     value_set_si(M->p[P->NbRays][dim+1], 1);
 256
 257     /* Delaunay triangulation */
 258     for (i = 0, r = 1; i < P->NbRays; ++i) {
 259         if (value_notzero_p(P->Ray[i][dim+1]))
 260             continue;
 261         Vector_Copy(P->Ray[i], M->p[r], dim+1);
 262         Inner_Product(M->p[r]+1, M->p[r]+1, dim, &tmp);
 263         value_assign(M->p[r][dim+1], tmp);
 264         value_set_si(M->p[r][dim+2], 0);
 265         ++r;
 266     }
 267
 268     M3 = Matrix_Copy(M);
 269     L = Rays2Polyhedron(M3, NbMaxCons);
 270     Matrix_Free(M3);
 271
 272     M2 = Matrix_Alloc(dim+1, dim+2);
 273
 274     t = 1;
 275     if (0) {
 276 try_again:
 277         /* Usually R should still be 0 */
 278         Domain_Free(R);
 279         Polyhedron_Free(L);
 280         for (r = 1; r < P->NbRays; ++r) {
 281             value_set_si(M->p[r][dim+1], random_int((t+1)*dim*P->NbRays)+1);
 282         }
 283         M3 = Matrix_Copy(M);
 284         L = Rays2Polyhedron(M3, NbMaxCons);
 285         Matrix_Free(M3);
 286         ++t;
 287     }
 288     assert(t <= MAX_TRY);
 289
 290     R = NULL;
 291     n = 0;
 292
 293     for (i = 0; i < L->NbConstraints; ++i) {
 294         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 295         if (value_negz_p(L->Constraint[i][dim+1]))
 296             continue;
 297         if (value_notzero_p(L->Constraint[i][dim+2]))
 298             continue;
 299         for (j = 1, r = 1; j < M->NbRows; ++j) {
 300             Inner_Product(M->p[j]+1, L->Constraint[i]+1, dim+1, &tmp);
 301             if (value_notzero_p(tmp))
 302                 continue;
 303             if (r > dim)
 304                 goto try_again;
 305             Vector_Copy(M->p[j]+1, M2->p[r]+1, dim);
 306             value_set_si(M2->p[r][0], 1);
 307             value_set_si(M2->p[r][dim+1], 0);
 308             ++r;
 309         }
 310         assert(r == dim+1);
 311         Vector_Set(M2->p[0]+1, 0, dim);
 312         value_set_si(M2->p[0][0], 1);
 313         value_set_si(M2->p[0][dim+1], 1);
 314         T = Rays2Polyhedron(M2, P->NbConstraints+1);
 315         T->next = R;
 316         R = T;
 317         ++n;
 318     }
 319     Matrix_Free(M2);
 320
 321     Polyhedron_Free(L);
 322     value_clear(tmp);
 323     Matrix_Free(M);
 324
 325     return R;
 326 }
 327
 328 void check_triangulization(Polyhedron *P, Polyhedron *T)
 329 {
 330     Polyhedron *C, *D, *E, *F, *G, *U;
 331     for (C = T; C; C = C->next) {
 332         if (C == T)
 333             U = C;
 334         else
 335             U = DomainConvex(DomainUnion(U, C, 100), 100);
 336         for (D = C->next; D; D = D->next) {
 337             F = C->next;
 338             G = D->next;
 339             C->next = NULL;
 340             D->next = NULL;
 341             E = DomainIntersection(C, D, 600);
 342             assert(E->NbRays == 0 || E->NbEq >= 1);
 343             Polyhedron_Free(E);
 344             C->next = F;
 345             D->next = G;
 346         }
 347     }
 348     assert(PolyhedronIncludes(U, P));
 349     assert(PolyhedronIncludes(P, U));
 350 }
 351
 352 /* Computes x, y and g such that g = gcd(a,b) and a*x+b*y = g */
 353 void Extended_Euclid(Value a, Value b, Value *x, Value *y, Value *g)
 354 {
 355     Value c, d, e, f, tmp;
 356
 357     value_init(c);
 358     value_init(d);
 359     value_init(e);
 360     value_init(f);
 361     value_init(tmp);
 362     value_absolute(c, a);
 363     value_absolute(d, b);
 364     value_set_si(e, 1);
 365     value_set_si(f, 0);
 366     while(value_pos_p(d)) {
 367         value_division(tmp, c, d);
 368         value_multiply(tmp, tmp, f);
 369         value_subtract(e, e, tmp);
 370         value_division(tmp, c, d);
 371         value_multiply(tmp, tmp, d);
 372         value_subtract(c, c, tmp);
 373         value_swap(c, d);
 374         value_swap(e, f);
 375     }
 376     value_assign(*g, c);
 377     if (value_zero_p(a))
 378         value_set_si(*x, 0);
 379     else if (value_pos_p(a))
 380         value_assign(*x, e);
 381     else value_oppose(*x, e);
 382     if (value_zero_p(b))
 383         value_set_si(*y, 0);
 384     else {
 385         value_multiply(tmp, a, *x);
 386         value_subtract(tmp, c, tmp);
 387         value_division(*y, tmp, b);
 388     }
 389     value_clear(c);
 390     value_clear(d);
 391     value_clear(e);
 392     value_clear(f);
 393     value_clear(tmp);
 394 }
 395
 396 Matrix * unimodular_complete(Vector *row)
 397 {
 398     Value g, b, c, old, tmp;
 399     Matrix *m;
 400     unsigned i, j;
 401
 402     value_init(b);
 403     value_init(c);
 404     value_init(g);
 405     value_init(old);
 406     value_init(tmp);
 407     m = Matrix_Alloc(row->Size, row->Size);
 408     for (j = 0; j < row->Size; ++j) {
 409         value_assign(m->p[0][j], row->p[j]);
 410     }
 411     value_assign(g, row->p[0]);
 412     for (i = 1; value_zero_p(g) && i < row->Size; ++i) {
 413         for (j = 0; j < row->Size; ++j) {
 414             if (j == i-1)
 415                 value_set_si(m->p[i][j], 1);
 416             else
 417                 value_set_si(m->p[i][j], 0);
 418         }
 419         value_assign(g, row->p[i]);
 420     }
 421     for (; i < row->Size; ++i) {
 422         value_assign(old, g);
 423         Extended_Euclid(old, row->p[i], &c, &b, &g);
 424         value_oppose(b, b);
 425         for (j = 0; j < row->Size; ++j) {
 426             if (j < i) {
 427                 value_multiply(tmp, row->p[j], b);
 428                 value_division(m->p[i][j], tmp, old);
 429             } else if (j == i)
 430                 value_assign(m->p[i][j], c);
 431             else
 432                 value_set_si(m->p[i][j], 0);
 433         }
 434     }
 435     value_clear(b);
 436     value_clear(c);
 437     value_clear(g);
 438     value_clear(old);
 439     value_clear(tmp);
 440     return m;
 441 }
 442
 443 /*
 444  * Returns a full-dimensional polyhedron with the same number
 445  * of integer points as P
 446  */
 447 Polyhedron *remove_equalities(Polyhedron *P)
 448 {
 449     Value g;
 450     Vector *v;
 451     Polyhedron *p = Polyhedron_Copy(P), *q;
 452     unsigned dim = p->Dimension;
 453     Matrix *m1, *m2;
 454     int i;
 455
 456     value_init(g);
 457     while (!emptyQ2(p) && p->NbEq > 0) {
 458         assert(dim > 0);
 459         Vector_Gcd(p->Constraint[0]+1, dim+1, &g);
 460         Vector_AntiScale(p->Constraint[0]+1, p->Constraint[0]+1, g, dim+1);
 461         Vector_Gcd(p->Constraint[0]+1, dim, &g);
 462         if (value_notone_p(g) && value_notmone_p(g)) {
 463             Polyhedron_Free(p);
 464             p = Empty_Polyhedron(0);
 465             break;
 466         }
 467         v = Vector_Alloc(dim);
 468         Vector_Copy(p->Constraint[0]+1, v->p, dim);
 469         m1 = unimodular_complete(v);
 470         m2 = Matrix_Alloc(dim, dim+1);
 471         for (i = 0; i < dim-1 ; ++i) {
 472             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 473             value_set_si(m2->p[i][dim], 0);
 474         }
 475         Vector_Set(m2->p[dim-1], 0, dim);
 476         value_set_si(m2->p[dim-1][dim], 1);
 477         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 478         Vector_Free(v);
 479         Matrix_Free(m1);
 480         Matrix_Free(m2);
 481         Polyhedron_Free(p);
 482         p = q;
 483         --dim;
 484     }
 485     value_clear(g);
 486     return p;
 487 }
 488
 489 /*
 490  * Returns a full-dimensional polyhedron with the same number
 491  * of integer points as P
 492  * nvar specifies the number of variables
 493  * The remaining dimensions are assumed to be parameters
 494  * Destroys P
 495  * factor is NbEq x (nparam+2) matrix, containing stride constraints
 496  * on the parameters; column nparam is the constant;
 497  * column nparam+1 is the stride
 498  *
 499  * if factor is NULL, only remove equalities that don't affect
 500  * the number of points
 501  */
 502 Polyhedron *remove_equalities_p(Polyhedron *P, unsigned nvar, Matrix **factor)
 503 {
 504     Value g;
 505     Vector *v;
 506     Polyhedron *p = P, *q;
 507     unsigned dim = p->Dimension;
 508     Matrix *m1, *m2, *f;
 509     int i, j, skip;
 510
 511     value_init(g);
 512     if (factor) {
 513         f = Matrix_Alloc(p->NbEq, dim-nvar+2);
 514         *factor = f;
 515     }
 516     j = 0;
 517     skip = 0;
 518     while (nvar > 0 && p->NbEq - skip > 0) {
 519         assert(dim > 0);
 520
 521         while (skip < p->NbEq &&
 522                First_Non_Zero(p->Constraint[skip]+1, nvar) == -1)
 523             ++skip;
 524         if (p->NbEq == skip)
 525             break;
 526
 527         Vector_Gcd(p->Constraint[skip]+1, dim+1, &g);
 528         Vector_AntiScale(p->Constraint[skip]+1, p->Constraint[skip]+1, g, dim+1);
 529         Vector_Gcd(p->Constraint[skip]+1, nvar, &g);
 530         if (!factor && value_notone_p(g) && value_notmone_p(g)) {
 531             ++skip;
 532             continue;
 533         }
 534         if (factor) {
 535             Vector_Copy(p->Constraint[skip]+1+nvar, f->p[j], dim-nvar+1);
 536             value_assign(f->p[j][dim-nvar+1], g);
 537         }
 538         v = Vector_Alloc(dim);
 539         Vector_AntiScale(p->Constraint[skip]+1, v->p, g, nvar);
 540         Vector_Set(v->p+nvar, 0, dim-nvar);
 541         m1 = unimodular_complete(v);
 542         m2 = Matrix_Alloc(dim, dim+1);
 543         for (i = 0; i < dim-1 ; ++i) {
 544             Vector_Copy(m1->p[i+1], m2->p[i], dim);
 545             value_set_si(m2->p[i][dim], 0);
 546         }
 547         Vector_Set(m2->p[dim-1], 0, dim);
 548         value_set_si(m2->p[dim-1][dim], 1);
 549         q = Polyhedron_Image(p, m2, p->NbConstraints+1+p->NbRays);
 550         Vector_Free(v);
 551         Matrix_Free(m1);
 552         Matrix_Free(m2);
 553         Polyhedron_Free(p);
 554         p = q;
 555         --dim;
 556         --nvar;
 557         ++j;
 558     }
 559     value_clear(g);
 560     return p;
 561 }
 562
 563 void Line_Length(Polyhedron *P, Value *len)
 564 {
 565     Value tmp, pos, neg;
 566     int p = 0, n = 0;
 567     int i;
 568
 569     assert(P->Dimension == 1);
 570
 571     value_init(tmp);
 572     value_init(pos);
 573     value_init(neg);
 574
 575     for (i = 0; i < P->NbConstraints; ++i) {
 576         value_oppose(tmp, P->Constraint[i][2]);
 577         if (value_pos_p(P->Constraint[i][1])) {
 578             mpz_cdiv_q(tmp, tmp, P->Constraint[i][1]);
 579             if (!p || value_gt(tmp, pos))
 580                 value_assign(pos, tmp);
 581             p = 1;
 582         } else {
 583             mpz_fdiv_q(tmp, tmp, P->Constraint[i][1]);
 584             if (!n || value_lt(tmp, neg))
 585                 value_assign(neg, tmp);
 586             n = 1;
 587         }
 588         if (n && p) {
 589             value_subtract(tmp, neg, pos);
 590             value_increment(*len, tmp);
 591         } else
 592             value_set_si(*len, -1);
 593     }
 594
 595     value_clear(tmp);
 596     value_clear(pos);
 597     value_clear(neg);
 598 }
 599
 600 /*
 601  * Factors the polyhedron P into polyhedra Q_i such that
 602  * the number of integer points in P is equal to the product
 603  * of the number of integer points in the individual Q_i
 604  *
 605  * If no factors can be found, NULL is returned.
 606  * Otherwise, a linked list of the factors is returned.
 607  *
 608  * The algorithm works by first computing the Hermite normal form
 609  * and then grouping columns linked by one or more constraints together,
 610  * where a constraints "links" two or more columns if the constraint
 611  * has nonzero coefficients in the columns.
 612  */
 613 Polyhedron* Polyhedron_Factor(Polyhedron *P, unsigned nparam,
 614                               unsigned NbMaxRays)
 615 {
 616     int i, j, k;
 617     Matrix *M, *H, *Q, *U;
 618     int *pos;           /* for each column: row position of pivot */
 619     int *group;         /* group to which a column belongs */
 620     int *cnt;           /* number of columns in the group */
 621     int *rowgroup;      /* group to which a constraint belongs */
 622     int nvar = P->Dimension - nparam;
 623     Polyhedron *F = NULL;
 624
 625     if (nvar <= 1)
 626         return NULL;
 627
 628     NALLOC(pos, nvar);
 629     NALLOC(group, nvar);
 630     NALLOC(cnt, nvar);
 631     NALLOC(rowgroup, P->NbConstraints);
 632
 633     M = Matrix_Alloc(P->NbConstraints, nvar);
 634     for (i = 0; i < P->NbConstraints; ++i)
 635         Vector_Copy(P->Constraint[i]+1, M->p[i], nvar);
 636     left_hermite(M, &H, &Q, &U);
 637     Matrix_Free(M);
 638     Matrix_Free(Q);
 639     Matrix_Free(U);
 640
 641     for (i = 0; i < P->NbConstraints; ++i)
 642         rowgroup[i] = -1;
 643     for (i = 0, j = 0; i < H->NbColumns; ++i) {
 644         for ( ; j < H->NbRows; ++j)
 645             if (value_notzero_p(H->p[j][i]))
 646                 break;
 647         assert (j < H->NbRows);
 648         pos[i] = j;
 649     }
 650     for (i = 0; i < nvar; ++i) {
 651         group[i] = i;
 652         cnt[i] = 1;
 653     }
 654     for (i = 0; i < H->NbColumns && cnt[0] < nvar; ++i) {
 655         if (rowgroup[pos[i]] == -1)
 656             rowgroup[pos[i]] = i;
 657         for (j = pos[i]+1; j <  H->NbRows; ++j) {
 658             if (value_zero_p(H->p[j][i]))
 659                 continue;
 660             if (rowgroup[j] != -1)
 661                 continue;
 662             rowgroup[j] = group[i];
 663             for (k = i+1; k < H->NbColumns && j >= pos[k]; ++k) {
 664                 int g = group[k];
 665                 while (cnt[g] == 0)
 666                     g = group[g];
 667                 group[k] = g;
 668                 if (group[k] != group[i] && value_notzero_p(H->p[j][k])) {
 669                     assert(cnt[group[k]] != 0);
 670                     assert(cnt[group[i]] != 0);
 671                     if (group[i] < group[k]) {
 672                         cnt[group[i]] += cnt[group[k]];
 673                         cnt[group[k]] = 0;
 674                         group[k] = group[i];
 675                     } else {
 676                         cnt[group[k]] += cnt[group[i]];
 677                         cnt[group[i]] = 0;
 678                         group[i] = group[k];
 679                     }
 680                 }
 681             }
 682         }
 683     }
 684
 685     if (cnt[0] != nvar) {
 686         /* Extract out pure context constraints separately */
 687         Polyhedron **next = &F;
 688         for (i = nparam ? -1 : 0; i < nvar; ++i) {
 689             int d;
 690
 691             if (i == -1) {
 692                 for (j = 0, k = 0; j < P->NbConstraints; ++j)
 693                     if (rowgroup[j] == -1) {
 694                         if (First_Non_Zero(P->Constraint[j]+1+nvar,
 695                                            nparam) == -1)
 696                             rowgroup[j] = -2;
 697                         else
 698                             ++k;
 699                     }
 700                 if (k == 0)
 701                     continue;
 702                 d = 0;
 703             } else {
 704                 if (cnt[i] == 0)
 705                     continue;
 706                 d = cnt[i];
 707                 for (j = 0, k = 0; j < P->NbConstraints; ++j)
 708                     if (rowgroup[j] >= 0 && group[rowgroup[j]] == i) {
 709                         rowgroup[j] = i;
 710                         ++k;
 711                     }
 712             }
 713
 714             M = Matrix_Alloc(k, d+nparam+2);
 715             for (j = 0, k = 0; j < P->NbConstraints; ++j) {
 716                 int l, m;
 717                 if (rowgroup[j] != i)
 718                     continue;
 719                 value_assign(M->p[k][0], P->Constraint[j][0]);
 720                 for (l = 0, m = 0; m < d; ++l) {
 721                     if (group[l] != i)
 722                         continue;
 723                     value_assign(M->p[k][1+m++], H->p[j][l]);
 724                 }
 725                 Vector_Copy(P->Constraint[j]+1+nvar, M->p[k]+1+m, nparam+1);
 726                 ++k;
 727             }
 728             *next = Constraints2Polyhedron(M, NbMaxRays);
 729             next = &(*next)->next;
 730             Matrix_Free(M);
 731         }
 732     }
 733     Matrix_Free(H);
 734     free(pos);
 735     free(group);
 736     free(cnt);
 737     free(rowgroup);
 738     return F;
 739 }
 740
 741 /*
 742  * Project on final dim dimensions
 743  */
 744 Polyhedron* Polyhedron_Project(Polyhedron *P, int dim)
 745 {
 746     int i;
 747     int remove = P->Dimension - dim;
 748     Matrix *T;
 749     Polyhedron *I;
 750
 751     if (P->Dimension == dim)
 752         return Polyhedron_Copy(P);
 753
 754     T = Matrix_Alloc(dim+1, P->Dimension+1);
 755     for (i = 0; i < dim+1; ++i)
 756         value_set_si(T->p[i][i+remove], 1);
 757     I = Polyhedron_Image(P, T, P->NbConstraints);
 758     Matrix_Free(T);
 759     return I;
 760 }
 761
 762 /* Constructs a new constraint that ensures that
 763  * the first constraint is (strictly) smaller than
 764  * the second.
 765  */
 766 static void smaller_constraint(Value *a, Value *b, Value *c, int pos, int shift,
 767                         int len, int strict, Value *tmp)
 768 {
 769     value_oppose(*tmp, b[pos+1]);
 770     value_set_si(c[0], 1);
 771     Vector_Combine(a+1+shift, b+1+shift, c+1, *tmp, a[pos+1], len-shift-1);
 772     if (strict)
 773         value_decrement(c[len-shift-1], c[len-shift-1]);
 774     ConstraintSimplify(c, c, len-shift, tmp);
 775 }
 776
 777 struct section { Polyhedron * D; evalue E; };
 778
 779 evalue * ParamLine_Length_mod(Polyhedron *P, Polyhedron *C, int MaxRays)
 780 {
 781     unsigned dim = P->Dimension;
 782     unsigned nvar = dim - C->Dimension;
 783     int *pos;
 784     int i, j, p, n, z;
 785     struct section *s;
 786     Matrix *M, *M2;
 787     int nd = 0;
 788     int k, l, k2, l2, q;
 789     evalue *L, *U;
 790     evalue *F;
 791     Value g;
 792     Polyhedron *T;
 793     evalue mone;
 794
 795     assert(nvar == 1);
 796
 797     NALLOC(pos, P->NbConstraints);
 798     value_init(g);
 799     value_init(mone.d);
 800     evalue_set_si(&mone, -1, 1);
 801
 802     for (i = 0, z = 0; i < P->NbConstraints; ++i)
 803         if (value_zero_p(P->Constraint[i][1]))
 804             ++z;
 805     /* put those with positive coefficients first; number: p */
 806     for (i = 0, p = 0, n = P->NbConstraints-z-1; i < P->NbConstraints; ++i)
 807         if (value_pos_p(P->Constraint[i][1]))
 808             pos[p++] = i;
 809         else if (value_neg_p(P->Constraint[i][1]))
 810             pos[n--] = i;
 811     n = P->NbConstraints-z-p;
 812     assert (p >= 1 && n >= 1);
 813     s = (struct section *) malloc(p * n * sizeof(struct section));
 814     M = Matrix_Alloc((p-1) + (n-1), dim-nvar+2);
 815     for (k = 0; k < p; ++k) {
 816         for (k2 = 0; k2 < p; ++k2) {
 817             if (k2 == k)
 818                 continue;
 819             q = k2 - (k2 > k);
 820             smaller_constraint(
 821                 P->Constraint[pos[k]],
 822                 P->Constraint[pos[k2]],
 823                 M->p[q], 0, nvar, dim+2, k2 > k, &g);
 824         }
 825         for (l = p; l < p+n; ++l) {
 826             for (l2 = p; l2 < p+n; ++l2) {
 827                 if (l2 == l)
 828                     continue;
 829                 q = l2-1 - (l2 > l);
 830                 smaller_constraint(
 831                     P->Constraint[pos[l2]],
 832                     P->Constraint[pos[l]],
 833                     M->p[q], 0, nvar, dim+2, l2 > l, &g);
 834             }
 835             M2 = Matrix_Copy(M);
 836             T = Constraints2Polyhedron(M2, P->NbRays);
 837             Matrix_Free(M2);
 838             s[nd].D = DomainIntersection(T, C, MaxRays);
 839             Domain_Free(T);
 840             POL_ENSURE_VERTICES(s[nd].D);
 841             if (emptyQ(s[nd].D)) {
 842                 Polyhedron_Free(s[nd].D);
 843                 continue;
 844             }
 845             L = bv_ceil3(P->Constraint[pos[k]]+1+nvar,
 846                       dim-nvar+1,
 847                       P->Constraint[pos[k]][0+1], s[nd].D);
 848             U = bv_ceil3(P->Constraint[pos[l]]+1+nvar,
 849                       dim-nvar+1,
 850                       P->Constraint[pos[l]][0+1], s[nd].D);
 851             eadd(L, U);
 852             eadd(&mone, U);
 853             emul(&mone, U);
 854             s[nd].E = *U;
 855             free_evalue_refs(L);
 856             free(L);
 857             free(U);
 858             ++nd;
 859         }
 860     }
 861
 862     Matrix_Free(M);
 863
 864     F = ALLOC(evalue);
 865     value_init(F->d);
 866     value_set_si(F->d, 0);
 867     F->x.p = new_enode(partition, 2*nd, dim-nvar);
 868     for (k = 0; k < nd; ++k) {
 869         EVALUE_SET_DOMAIN(F->x.p->arr[2*k], s[k].D);
 870         value_clear(F->x.p->arr[2*k+1].d);
 871         F->x.p->arr[2*k+1] = s[k].E;
 872     }
 873     free(s);
 874
 875     free_evalue_refs(&mone);
 876     value_clear(g);
 877     free(pos);
 878
 879     return F;
 880 }
 881
 882 #ifdef USE_MODULO
 883 evalue* ParamLine_Length(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
 884 {
 885     return ParamLine_Length_mod(P, C, MaxRays);
 886 }
 887 #else
 888 evalue* ParamLine_Length(Polyhedron *P, Polyhedron *C, unsigned MaxRays)
 889 {
 890     evalue* tmp;
 891     tmp = ParamLine_Length_mod(P, C, MaxRays);
 892     evalue_mod2table(tmp, C->Dimension);
 893     reduce_evalue(tmp);
 894     return tmp;
 895 }
 896 #endif
 897
 898 Bool isIdentity(Matrix *M)
 899 {
 900     unsigned i, j;
 901     if (M->NbRows != M->NbColumns)
 902         return False;
 903
 904     for (i = 0;i < M->NbRows; i ++)
 905         for (j = 0; j < M->NbColumns; j ++)
 906             if (i == j) {
 907                 if(value_notone_p(M->p[i][j]))
 908                     return False;
 909             } else {
 910                 if(value_notzero_p(M->p[i][j]))
 911                     return False;
 912             }
 913     return True;
 914 }
 915
 916 void Param_Polyhedron_Print(FILE* DST, Param_Polyhedron *PP, char **param_names)
 917 {
 918   Param_Domain *P;
 919   Param_Vertices *V;
 920
 921   for(P=PP->D;P;P=P->next) {
 922
 923     /* prints current val. dom. */
 924     printf( "---------------------------------------\n" );
 925     printf( "Domain :\n");
 926     Print_Domain( stdout, P->Domain, param_names );
 927
 928     /* scan the vertices */
 929     printf( "Vertices :\n");
 930     FORALL_PVertex_in_ParamPolyhedron(V,P,PP) {
 931
 932       /* prints each vertex */
 933       Print_Vertex( stdout, V->Vertex, param_names );
 934       printf( "\n" );
 935     }
 936     END_FORALL_PVertex_in_ParamPolyhedron;
 937   }
 938 }
 939
 940 void Enumeration_Print(FILE *Dst, Enumeration *en, char **params)
 941 {
 942     for (; en; en = en->next) {
 943         Print_Domain(Dst, en->ValidityDomain, params);
 944         print_evalue(Dst, &en->EP, params);
 945     }
 946 }
 947
 948 void Enumeration_Free(Enumeration *en)
 949 {
 950   Enumeration *ee;
 951
 952   while( en )
 953   {
 954           free_evalue_refs( &(en->EP) );
 955           Domain_Free( en->ValidityDomain );
 956           ee = en ->next;
 957           free( en );
 958           en = ee;
 959   }
 960 }
 961
 962 void Enumeration_mod2table(Enumeration *en, unsigned nparam)
 963 {
 964     for (; en; en = en->next) {
 965         evalue_mod2table(&en->EP, nparam);
 966         reduce_evalue(&en->EP);
 967     }
 968 }
 969
 970 size_t Enumeration_size(Enumeration *en)
 971 {
 972     size_t s = 0;
 973
 974     for (; en; en = en->next) {
 975         s += domain_size(en->ValidityDomain);
 976         s += evalue_size(&en->EP);
 977     }
 978     return s;
 979 }
 980
 981 void Free_ParamNames(char **params, int m)
 982 {
 983     while (--m >= 0)
 984         free(params[m]);
 985     free(params);
 986 }
 987
 988 int DomainIncludes(Polyhedron *Pol1, Polyhedron *Pol2)
 989 {
 990     Polyhedron *P2;
 991     for ( ; Pol1; Pol1 = Pol1->next) {
 992         for (P2 = Pol2; P2; P2 = P2->next)
 993             if (!PolyhedronIncludes(Pol1, P2))
 994                 break;
 995         if (!P2)
 996             return 1;
 997     }
 998     return 0;
 999 }
1000
1001 int line_minmax(Polyhedron *I, Value *min, Value *max)
1002 {
1003     int i;
1004
1005     if (I->NbEq >= 1) {
1006         value_oppose(I->Constraint[0][2], I->Constraint[0][2]);
1007         /* There should never be a remainder here */
1008         if (value_pos_p(I->Constraint[0][1]))
1009             mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
1010         else
1011             mpz_fdiv_q(*min, I->Constraint[0][2], I->Constraint[0][1]);
1012         value_assign(*max, *min);
1013     } else for (i = 0; i < I->NbConstraints; ++i) {
1014         if (value_zero_p(I->Constraint[i][1])) {
1015             Polyhedron_Free(I);
1016             return 0;
1017         }
1018
1019         value_oppose(I->Constraint[i][2], I->Constraint[i][2]);
1020         if (value_pos_p(I->Constraint[i][1]))
1021             mpz_cdiv_q(*min, I->Constraint[i][2], I->Constraint[i][1]);
1022         else
1023             mpz_fdiv_q(*max, I->Constraint[i][2], I->Constraint[i][1]);
1024     }
1025     Polyhedron_Free(I);
1026     return 1;
1027 }
1028
1029 /**
1030
1031 PROCEDURES TO COMPUTE ENUMERATION. recursive procedure, recurse for
1032 each imbriquation
1033
1034 @param pos index position of current loop index (1..hdim-1)
1035 @param P loop domain
1036 @param context context values for fixed indices
1037 @param exist    number of existential variables
1038 @return the number of integer points in this
1039 polyhedron
1040
1041 */
1042 void count_points_e (int pos, Polyhedron *P, int exist, int nparam,
1043                      Value *context, Value *res)
1044 {
1045     Value LB, UB, k, c;
1046
1047     if (emptyQ(P)) {
1048         value_set_si(*res, 0);
1049         return;
1050     }
1051
1052     value_init(LB); value_init(UB); value_init(k);
1053     value_set_si(LB,0);
1054     value_set_si(UB,0);
1055
1056     if (lower_upper_bounds(pos,P,context,&LB,&UB) !=0) {
1057         /* Problem if UB or LB is INFINITY */
1058         value_clear(LB); value_clear(UB); value_clear(k);
1059         if (pos > P->Dimension - nparam - exist)
1060             value_set_si(*res, 1);
1061         else
1062             value_set_si(*res, -1);
1063         return;
1064     }
1065
1066 #ifdef EDEBUG1
1067     if (!P->next) {
1068         int i;
1069         for (value_assign(k,LB); value_le(k,UB); value_increment(k,k)) {
1070             fprintf(stderr, "(");
1071             for (i=1; i<pos; i++) {
1072                 value_print(stderr,P_VALUE_FMT,context[i]);
1073                 fprintf(stderr,",");
1074             }
1075             value_print(stderr,P_VALUE_FMT,k);
1076             fprintf(stderr,")\n");
1077         }
1078     }
1079 #endif
1080
1081     value_set_si(context[pos],0);
1082     if (value_lt(UB,LB)) {
1083         value_clear(LB); value_clear(UB); value_clear(k);
1084         value_set_si(*res, 0);
1085         return;
1086     }
1087     if (!P->next) {
1088         if (exist)
1089             value_set_si(*res, 1);
1090         else {
1091             value_subtract(k,UB,LB);
1092             value_add_int(k,k,1);
1093             value_assign(*res, k);
1094         }
1095         value_clear(LB); value_clear(UB); value_clear(k);
1096         return;
1097     }
1098
1099     /*-----------------------------------------------------------------*/
1100     /* Optimization idea                                               */
1101     /*   If inner loops are not a function of k (the current index)    */
1102     /*   i.e. if P->Constraint[i][pos]==0 for all P following this and */
1103     /*        for all i,                                               */
1104     /*   Then CNT = (UB-LB+1)*count_points(pos+1, P->next, context)    */
1105     /*   (skip the for loop)                                           */
1106     /*-----------------------------------------------------------------*/
1107
1108     value_init(c);
1109     value_set_si(*res, 0);
1110     for (value_assign(k,LB);value_le(k,UB);value_increment(k,k)) {
1111         /* Insert k in context */
1112         value_assign(context[pos],k);
1113         count_points_e(pos+1, P->next, exist, nparam, context, &c);
1114         if(value_notmone_p(c))
1115             value_addto(*res, *res, c);
1116         else {
1117             value_set_si(*res, -1);
1118             break;
1119         }
1120         if (pos > P->Dimension - nparam - exist &&
1121                 value_pos_p(*res))
1122             break;
1123     }
1124     value_clear(c);
1125
1126 #ifdef EDEBUG11
1127     fprintf(stderr,"%d\n",CNT);
1128 #endif
1129
1130     /* Reset context */
1131     value_set_si(context[pos],0);
1132     value_clear(LB); value_clear(UB); value_clear(k);
1133     return;
1134 } /* count_points_e */
1135
1136 int DomainContains(Polyhedron *P, Value *list_args, int len,
1137                    unsigned MaxRays, int set)
1138 {
1139     int i;
1140     Value m;
1141
1142     if (P->Dimension == len)
1143         return in_domain(P, list_args);
1144
1145     assert(set); // assume list_args is large enough
1146     assert((P->Dimension - len) % 2 == 0);
1147     value_init(m);
1148     for (i = 0; i < P->Dimension - len; i += 2) {
1149         int j, k;
1150         for (j = 0 ; j < P->NbEq; ++j)
1151             if (value_notzero_p(P->Constraint[j][1+len+i]))
1152                 break;
1153         assert(j < P->NbEq);
1154         value_absolute(m, P->Constraint[j][1+len+i]);
1155         k = First_Non_Zero(P->Constraint[j]+1, len);
1156         assert(k != -1);
1157         assert(First_Non_Zero(P->Constraint[j]+1+k+1, len - k - 1) == -1);
1158         mpz_fdiv_q(list_args[len+i], list_args[k], m);
1159         mpz_fdiv_r(list_args[len+i+1], list_args[k], m);
1160     }
1161     value_clear(m);
1162
1163     return in_domain(P, list_args);
1164 }
1165
1166 Polyhedron *DomainConcat(Polyhedron *head, Polyhedron *tail)
1167 {
1168     Polyhedron *S;
1169     if (!head)
1170         return tail;
1171     for (S = head; S->next; S = S->next)
1172         ;
1173     S->next = tail;
1174     return head;
1175 }
1176
1177 #ifndef HAVE_LEXSMALLER
1178
1179 evalue *barvinok_lexsmaller_ev(Polyhedron *P, Polyhedron *D, unsigned dim,
1180                                Polyhedron *C, unsigned MaxRays)
1181 {
1182     assert(0);
1183 }
1184
1185 #else
1186 #include <polylib/ranking.h>
1187
1188 evalue *barvinok_lexsmaller_ev(Polyhedron *P, Polyhedron *D, unsigned dim,
1189                                Polyhedron *C, unsigned MaxRays)
1190 {
1191     evalue *ranking;
1192     Polyhedron *RC, *RD, *Q;
1193     unsigned nparam = dim + C->Dimension;
1194     unsigned exist;
1195     Polyhedron *CA;
1196
1197     RC = LexSmaller(P, D, dim, C, MaxRays);
1198     RD = RC->next;
1199     RC->next = NULL;
1200
1201     exist = RD->Dimension - nparam - dim;
1202     CA = align_context(RC, RD->Dimension, MaxRays);
1203     Q = DomainIntersection(RD, CA, MaxRays);
1204     Polyhedron_Free(CA);
1205     Domain_Free(RD);
1206     Polyhedron_Free(RC);
1207     RD = Q;
1208
1209     for (Q = RD; Q; Q = Q->next) {
1210         evalue *t;
1211         Polyhedron *next = Q->next;
1212         Q->next = 0;
1213
1214         t = barvinok_enumerate_e(Q, exist, nparam, MaxRays);
1215
1216         if (Q == RD)
1217             ranking = t;
1218         else {
1219             eadd(t, ranking);
1220             free_evalue_refs(t);
1221             free(t);
1222         }
1223
1224         Q->next = next;
1225     }
1226
1227     Domain_Free(RD);
1228
1229     return ranking;
1230 }
1231
1232 Enumeration *barvinok_lexsmaller(Polyhedron *P, Polyhedron *D, unsigned dim,
1233                                  Polyhedron *C, unsigned MaxRays)
1234 {
1235     evalue *EP = barvinok_lexsmaller_ev(P, D, dim, C, MaxRays);
1236
1237     return partition2enumeration(EP);
1238 }
1239 #endif
1240
1241 /* "align" matrix to have nrows by inserting
1242  * the necessary number of rows and an equal number of columns in front
1243  */
1244 Matrix *align_matrix(Matrix *M, int nrows)
1245 {
1246     int i;
1247     int newrows = nrows - M->NbRows;
1248     Matrix *M2 = Matrix_Alloc(nrows, newrows + M->NbColumns);
1249     for (i = 0; i < newrows; ++i)
1250         value_set_si(M2->p[i][i], 1);
1251     for (i = 0; i < M->NbRows; ++i)
1252         Vector_Copy(M->p[i], M2->p[newrows+i]+newrows, M->NbColumns);
1253     return M2;
1254 }
1255
1256 static void print_varlist(FILE *out, int n, char **names)
1257 {
1258     int i;
1259     fprintf(out, "[");
1260     for (i = 0; i < n; ++i) {
1261         if (i)
1262             fprintf(out, ",");
1263         fprintf(out, "%s", names[i]);
1264     }
1265     fprintf(out, "]");
1266 }
1267
1268 static void print_term(FILE *out, Value v, int pos, int dim, int nparam,
1269                        char **iter_names, char **param_names, int *first)
1270 {
1271     if (value_zero_p(v)) {
1272         if (first && *first && pos >= dim + nparam)
1273             fprintf(out, "0");
1274         return;
1275     }
1276
1277     if (first) {
1278         if (!*first && value_pos_p(v))
1279             fprintf(out, "+");
1280         *first = 0;
1281     }
1282     if (pos < dim + nparam) {
1283         if (value_mone_p(v))
1284             fprintf(out, "-");
1285         else if (!value_one_p(v))
1286             value_print(out, VALUE_FMT, v);
1287         if (pos < dim)
1288             fprintf(out, "%s", iter_names[pos]);
1289         else
1290             fprintf(out, "%s", param_names[pos-dim]);
1291     } else
1292         value_print(out, VALUE_FMT, v);
1293 }
1294
1295 char **util_generate_names(int n, char *prefix)
1296 {
1297     int i;
1298     int len = (prefix ? strlen(prefix) : 0) + 10;
1299     char **names = ALLOCN(char*, n);
1300     if (!names) {
1301         fprintf(stderr, "ERROR: memory overflow.\n");
1302         exit(1);
1303     }
1304     for (i = 0; i < n; ++i) {
1305         names[i] = ALLOCN(char, len);
1306         if (!names[i]) {
1307             fprintf(stderr, "ERROR: memory overflow.\n");
1308             exit(1);
1309         }
1310         if (!prefix)
1311             snprintf(names[i], len, "%d", i);
1312         else
1313             snprintf(names[i], len, "%s%d", prefix, i);
1314     }
1315
1316     return names;
1317 }
1318
1319 void util_free_names(int n, char **names)
1320 {
1321     int i;
1322     for (i = 0; i < n; ++i)
1323         free(names[i]);
1324     free(names);
1325 }
1326
1327 void Polyhedron_pprint(FILE *out, Polyhedron *P, int dim, int nparam,
1328                        char **iter_names, char **param_names)
1329 {
1330     int i, j;
1331     Value tmp;
1332
1333     assert(dim + nparam == P->Dimension);
1334
1335     value_init(tmp);
1336
1337     fprintf(out, "{ ");
1338     if (nparam) {
1339         print_varlist(out, nparam, param_names);
1340         fprintf(out, " -> ");
1341     }
1342     print_varlist(out, dim, iter_names);
1343     fprintf(out, " : ");
1344
1345     if (emptyQ2(P))
1346         fprintf(out, "FALSE");
1347     else for (i = 0; i < P->NbConstraints; ++i) {
1348         int first = 1;
1349         int v = First_Non_Zero(P->Constraint[i]+1, P->Dimension);
1350         if (v == -1 && value_pos_p(P->Constraint[i][0]))
1351             continue;
1352         if (i)
1353             fprintf(out, " && ");
1354         if (v == -1 && value_notzero_p(P->Constraint[i][1+P->Dimension]))
1355             fprintf(out, "FALSE");
1356         else if (value_pos_p(P->Constraint[i][v+1])) {
1357             print_term(out, P->Constraint[i][v+1], v, dim, nparam,
1358                        iter_names, param_names, NULL);
1359             if (value_zero_p(P->Constraint[i][0]))
1360                 fprintf(out, " = ");
1361             else
1362                 fprintf(out, " >= ");
1363             for (j = v+1; j <= dim+nparam; ++j) {
1364                 value_oppose(tmp, P->Constraint[i][1+j]);
1365                 print_term(out, tmp, j, dim, nparam,
1366                            iter_names, param_names, &first);
1367             }
1368         } else {
1369             value_oppose(tmp, P->Constraint[i][1+v]);
1370             print_term(out, tmp, v, dim, nparam,
1371                        iter_names, param_names, NULL);
1372             fprintf(out, " <= ");
1373             for (j = v+1; j <= dim+nparam; ++j)
1374                 print_term(out, P->Constraint[i][1+j], j, dim, nparam,
1375                            iter_names, param_names, &first);
1376         }
1377     }
1378
1379     fprintf(out, " }\n");
1380
1381     value_clear(tmp);
1382 }
1383
1384 /* Construct a cone over P with P placed at x_d = 1, with
1385  * x_d the coordinate of an extra dimension
1386  *
1387  * It's probably a mistake to depend so much on the internal
1388  * representation.  We should probably simply compute the
1389  * vertices/facets first.
1390  */
1391 Polyhedron *Cone_over_Polyhedron(Polyhedron *P)
1392 {
1393     unsigned NbConstraints = 0;
1394     unsigned NbRays = 0;
1395     Polyhedron *C;
1396     int i;
1397
1398     if (POL_HAS(P, POL_INEQUALITIES))
1399         NbConstraints = P->NbConstraints + 1;
1400     if (POL_HAS(P, POL_POINTS))
1401         NbRays = P->NbRays + 1;
1402
1403     C = Polyhedron_Alloc(P->Dimension+1, NbConstraints, NbRays);
1404     if (POL_HAS(P, POL_INEQUALITIES)) {
1405         C->NbEq = P->NbEq;
1406         for (i = 0; i < P->NbConstraints; ++i)
1407             Vector_Copy(P->Constraint[i], C->Constraint[i], P->Dimension+2);
1408         /* n >= 0 */
1409         value_set_si(C->Constraint[P->NbConstraints][0], 1);
1410         value_set_si(C->Constraint[P->NbConstraints][1+P->Dimension], 1);
1411     }
1412     if (POL_HAS(P, POL_POINTS)) {
1413         C->NbBid = P->NbBid;
1414         for (i = 0; i < P->NbRays; ++i)
1415             Vector_Copy(P->Ray[i], C->Ray[i], P->Dimension+2);
1416         /* vertex 0 */
1417         value_set_si(C->Ray[P->NbRays][0], 1);
1418         value_set_si(C->Ray[P->NbRays][1+C->Dimension], 1);
1419     }
1420     POL_SET(C, POL_VALID);
1421     if (POL_HAS(P, POL_INEQUALITIES))
1422         POL_SET(C, POL_INEQUALITIES);
1423     if (POL_HAS(P, POL_POINTS))
1424         POL_SET(C, POL_POINTS);
1425     if (POL_HAS(P, POL_VERTICES))
1426         POL_SET(C, POL_VERTICES);
1427     return C;
1428 }
1429
1430 /* Returns a (dim+nparam+1)x((dim-n)+nparam+1) matrix
1431  * mapping the transformed subspace back to the original space.
1432  * n is the number of equalities involving the variables
1433  * (i.e., not purely the parameters).
1434  * The remaining n coordinates in the transformed space would
1435  * have constant (parametric) values and are therefore not
1436  * included in the variables of the new space.
1437  */
1438 Matrix *compress_variables(Matrix *Equalities, unsigned nparam)
1439 {
1440     unsigned dim = (Equalities->NbColumns-2) - nparam;
1441     Matrix *M, *H, *Q, *U, *C, *ratH, *invH, *Ul, *T1, *T2, *T;
1442     Value mone;
1443     int n, i, j;
1444     int ok;
1445
1446     for (n = 0; n < Equalities->NbRows; ++n)
1447         if (First_Non_Zero(Equalities->p[n]+1, dim) == -1)
1448             break;
1449     if (n == 0)
1450         return Identity(dim+nparam+1);
1451     value_init(mone);
1452     value_set_si(mone, -1);
1453     M = Matrix_Alloc(n, dim);
1454     C = Matrix_Alloc(n+1, nparam+1);
1455     for (i = 0; i < n; ++i) {
1456         Vector_Copy(Equalities->p[i]+1, M->p[i], dim);
1457         Vector_Scale(Equalities->p[i]+1+dim, C->p[i], mone, nparam+1);
1458     }
1459     value_set_si(C->p[n][nparam], 1);
1460     left_hermite(M, &H, &Q, &U);
1461     Matrix_Free(M);
1462     Matrix_Free(Q);
1463     value_clear(mone);
1464
1465     ratH = Matrix_Alloc(n+1, n+1);
1466     invH = Matrix_Alloc(n+1, n+1);
1467     for (i = 0; i < n; ++i)
1468         Vector_Copy(H->p[i], ratH->p[i], n);
1469     value_set_si(ratH->p[n][n], 1);
1470     ok = Matrix_Inverse(ratH, invH);
1471     assert(ok);
1472     Matrix_Free(H);
1473     Matrix_Free(ratH);
1474     T1 = Matrix_Alloc(n+1, nparam+1);
1475     Matrix_Product(invH, C, T1);
1476     Matrix_Free(C);
1477     Matrix_Free(invH);
1478     if (value_notone_p(T1->p[n][nparam])) {
1479         for (i = 0; i < n; ++i) {
1480             if (!mpz_divisible_p(T1->p[i][nparam], T1->p[n][nparam])) {
1481                 Matrix_Free(T1);
1482                 Matrix_Free(U);
1483                 return NULL;
1484             }
1485             /* compress_params should have taken care of this */
1486             for (j = 0; j < nparam; ++j)
1487                 assert(mpz_divisible_p(T1->p[i][j], T1->p[n][nparam]));
1488             Vector_AntiScale(T1->p[i], T1->p[i], T1->p[n][nparam], nparam+1);
1489         }
1490         value_set_si(T1->p[n][nparam], 1);
1491     }
1492     Ul = Matrix_Alloc(dim+1, n+1);
1493     for (i = 0; i < dim; ++i)
1494         Vector_Copy(U->p[i], Ul->p[i], n);
1495     value_set_si(Ul->p[dim][n], 1);
1496     T2 = Matrix_Alloc(dim+1, nparam+1);
1497     Matrix_Product(Ul, T1, T2);
1498     Matrix_Free(Ul);
1499     Matrix_Free(T1);
1500
1501     T = Matrix_Alloc(dim+nparam+1, (dim-n)+nparam+1);
1502     for (i = 0; i < dim; ++i) {
1503         Vector_Copy(U->p[i]+n, T->p[i], dim-n);
1504         Vector_Copy(T2->p[i], T->p[i]+dim-n, nparam+1);
1505     }
1506     for (i = 0; i < nparam+1; ++i)
1507         value_set_si(T->p[dim+i][(dim-n)+i], 1);
1508     assert(value_one_p(T2->p[dim][nparam]));
1509     Matrix_Free(U);
1510     Matrix_Free(T2);
1511
1512     return T;
1513 }