util.c

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