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