util.c

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