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