util.c

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