util.c

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