util.c

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