scale.c

   1 #include <assert.h>
   2 #include <barvinok/barvinok.h>
   3 #include <barvinok/util.h>
   4 #include <barvinok/options.h>
   5 #include "scale.h"
   6 #include "reduce_domain.h"
   7 #include "param_util.h"
   8
   9 #define ALLOC(type) (type*)malloc(sizeof(type))
  10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  11
  12 /* If a vertex is described by A x + B p + c = 0, then
  13  * M = [A B] and we want to compute a linear transformation L such
  14  * that H L = A and H \Z contains both A \Z and B \Z.
  15  * We compute
  16  *             [ A B ] = [ H 0 ] [ U_11  U_12 ]
  17  *                               [ U_21  U_22 ]
  18  *
  19  * U_11 is the required linear transformation.
  20  * Note that this also works if M has more rows than there are variables,
  21  * i.e., if some rows in M are linear combinations of other rows.
  22  * These extra rows only affect and H and not U.
  23  */
  24 static Lattice *extract_lattice(Matrix *M, unsigned nvar)
  25 {
  26     int row;
  27     Matrix *H, *Q, *U, *Li;
  28     Lattice *L;
  29     int ok;
  30
  31     left_hermite(M, &H, &Q, &U);
  32     Matrix_Free(U);
  33
  34     Li = Matrix_Alloc(nvar+1, nvar+1);
  35     L = Matrix_Alloc(nvar+1, nvar+1);
  36     value_set_si(Li->p[nvar][nvar], 1);
  37
  38     for (row = 0; row < nvar; ++row)
  39         Vector_Copy(Q->p[row], Li->p[row], nvar);
  40     Matrix_Free(H);
  41     Matrix_Free(Q);
  42
  43     ok = Matrix_Inverse(Li, L);
  44     assert(ok);
  45     Matrix_Free(Li);
  46
  47     return L;
  48 }
  49
  50 /* Returns the smallest (wrt inclusion) lattice that contains both X and Y */
  51 static Lattice *LatticeJoin(Lattice *X, Lattice *Y)
  52 {
  53     int i;
  54     int dim = X->NbRows-1;
  55     Value lcm;
  56     Value tmp;
  57     Lattice *L;
  58     Matrix *M, *H, *U, *Q;
  59
  60     assert(X->NbColumns-1 == dim);
  61     assert(Y->NbRows-1 == dim);
  62     assert(Y->NbColumns-1 == dim);
  63
  64     value_init(lcm);
  65     value_init(tmp);
  66
  67     M = Matrix_Alloc(dim, 2*dim);
  68     value_lcm(X->p[dim][dim], Y->p[dim][dim], &lcm);
  69
  70     value_division(tmp, lcm, X->p[dim][dim]);
  71     for (i = 0; i < dim; ++i)
  72         Vector_Scale(X->p[i], M->p[i], tmp, dim);
  73     value_division(tmp, lcm, Y->p[dim][dim]);
  74     for (i = 0; i < dim; ++i)
  75         Vector_Scale(Y->p[i], M->p[i]+dim, tmp, dim);
  76
  77     left_hermite(M, &H, &Q, &U);
  78     Matrix_Free(M);
  79     Matrix_Free(Q);
  80     Matrix_Free(U);
  81
  82     L = Matrix_Alloc(dim+1, dim+1);
  83     value_assign(L->p[dim][dim], lcm);
  84     for (i = 0; i < dim; ++i)
  85         Vector_Copy(H->p[i], L->p[i], dim);
  86     Matrix_Free(H);
  87
  88     value_clear(tmp);
  89     value_clear(lcm);
  90     return L;
  91 }
  92
  93 static void Param_Vertex_Image(Param_Vertices *V, Matrix *T)
  94 {
  95     unsigned nvar  = V->Vertex->NbRows;
  96     unsigned nparam = V->Vertex->NbColumns - 2;
  97     Matrix *Vertex;
  98     int i;
  99
 100     Param_Vertex_Common_Denominator(V);
 101     Vertex = Matrix_Alloc(V->Vertex->NbRows, V->Vertex->NbColumns);
 102     Matrix_Product(T, V->Vertex, Vertex);
 103     for (i = 0; i < nvar; ++i) {
 104         value_assign(Vertex->p[i][nparam+1], V->Vertex->p[i][nparam+1]);
 105         Vector_Normalize(Vertex->p[i], nparam+2);
 106     }
 107     Matrix_Free(V->Vertex);
 108     V->Vertex = Vertex;
 109 }
 110
 111 /* Scales the parametric polyhedron with constraints *P and vertices PP
 112  * such that the number of integer points can be represented by a polynomial.
 113  * Both *P and P->Vertex are adapted according to the scaling.
 114  * The scaling factor is returned in *det.
 115  * The transformation that maps the new coordinates to the original
 116  * coordinates is returned in *Lat (if Lat != NULL).
 117  * The enumerator of the scaled parametric polyhedron should be divided
 118  * by this number to obtain an approximation of the enumerator of the
 119  * original parametric polyhedron.
 120  *
 121  * The algorithm is described in "Approximating Ehrhart Polynomials using
 122  * affine transformations" by B. Meister.
 123  */
 124 void Param_Polyhedron_Scale_Integer_Slow(Param_Polyhedron *PP, Polyhedron **P,
 125                                          Lattice **Lat,
 126                                          Value *det, unsigned MaxRays)
 127 {
 128     Param_Vertices *V;
 129     unsigned dim = (*P)->Dimension;
 130     unsigned nparam;
 131     unsigned nvar;
 132     Lattice *L = NULL, *Li;
 133     Matrix *T;
 134     Matrix *expansion;
 135     int i;
 136     int ok;
 137
 138     if (!PP->nbV)
 139         return;
 140
 141     nparam = PP->V->Vertex->NbColumns - 2;
 142     nvar = dim - nparam;
 143
 144     for (V = PP->V; V; V = V->next) {
 145         Lattice *L2;
 146         Matrix *M;
 147         int i, j, n;
 148         unsigned char *supporting;
 149
 150         supporting = supporting_constraints(*P, V, &n);
 151         M = Matrix_Alloc(n, (*P)->Dimension);
 152         for (i = 0, j = 0; i < (*P)->NbConstraints; ++i)
 153             if (supporting[i])
 154                 Vector_Copy((*P)->Constraint[i]+1, M->p[j++], (*P)->Dimension);
 155         free(supporting);
 156         L2 = extract_lattice(M, nvar);
 157         Matrix_Free(M);
 158
 159         if (!L)
 160             L = L2;
 161         else {
 162             Lattice *L3 = LatticeJoin(L, L2);
 163             Matrix_Free(L);
 164             Matrix_Free(L2);
 165             L = L3;
 166         }
 167     }
 168
 169     if (Lat)
 170         *Lat = Matrix_Copy(L);
 171
 172     /* apply the variable expansion to the polyhedron (constraints) */
 173     expansion = Matrix_Alloc(nvar + nparam + 1,  nvar + nparam + 1);
 174     for (i = 0; i < nvar; ++i)
 175         Vector_Copy(L->p[i], expansion->p[i], nvar);
 176     for (i = nvar; i < nvar+nparam+1; ++i)
 177         value_assign(expansion->p[i][i], L->p[nvar][nvar]);
 178
 179     *P = Polyhedron_Preimage(*P, expansion, MaxRays);
 180     Matrix_Free(expansion);
 181
 182     /* apply the variable expansion to the parametric vertices */
 183     Li = Matrix_Alloc(nvar+1, nvar+1);
 184     ok = Matrix_Inverse(L, Li);
 185     assert(ok);
 186     Matrix_Free(L);
 187     assert(value_one_p(Li->p[nvar][nvar]));
 188     T = Matrix_Alloc(nvar, nvar);
 189     value_set_si(*det, 1);
 190     for (i = 0; i < nvar; ++i) {
 191         value_multiply(*det, *det, Li->p[i][i]);
 192         Vector_Copy(Li->p[i], T->p[i], nvar);
 193     }
 194     Matrix_Free(Li);
 195     for (V = PP->V; V; V = V->next)
 196         Param_Vertex_Image(V, T);
 197     Matrix_Free(T);
 198 }
 199
 200 /* Scales the parametric polyhedron with constraints *P and vertices PP
 201  * such that the number of integer points can be represented by a polynomial.
 202  * Both *P and P->Vertex are adapted according to the scaling.
 203  * The scaling factor is returned in *det.
 204  * The transformation that maps the new coordinates to the original
 205  * coordinates is returned in *Lat (if Lat != NULL).
 206  * The enumerator of the scaled parametric polyhedron should be divided
 207  * by this number to obtain an approximation of the enumerator of the
 208  * original parametric polyhedron.
 209  *
 210  * The algorithm is described in "Approximating Ehrhart Polynomials using
 211  * affine transformations" by B. Meister.
 212  */
 213 void Param_Polyhedron_Scale_Integer_Fast(Param_Polyhedron *PP, Polyhedron **P,
 214                                          Lattice **Lat,
 215                                          Value *det, unsigned MaxRays)
 216 {
 217   int i;
 218   int nb_param, nb_vars;
 219   Vector *denoms;
 220   Param_Vertices *V;
 221   Value global_var_lcm;
 222   Value tmp;
 223   Matrix *expansion;
 224
 225   value_set_si(*det, 1);
 226   if (!PP->nbV)
 227     return;
 228
 229   nb_param = PP->D->Domain->Dimension;
 230   nb_vars = PP->V->Vertex->NbRows;
 231
 232   /* Scan the vertices and make an orthogonal expansion of the variable
 233      space */
 234   /* a- prepare the array of common denominators */
 235   denoms = Vector_Alloc(nb_vars);
 236   value_init(global_var_lcm);
 237
 238   value_init(tmp);
 239   /* b- scan the vertices and compute the variables' global lcms */
 240   for (V = PP->V; V; V = V->next) {
 241     for (i = 0; i < nb_vars; i++) {
 242       Vector_Gcd(V->Vertex->p[i], nb_param, &tmp);
 243       Gcd(tmp, V->Vertex->p[i][nb_param+1], &tmp);
 244       value_division(tmp, V->Vertex->p[i][nb_param+1], tmp);
 245       Lcm3(denoms->p[i], tmp, &denoms->p[i]);
 246     }
 247   }
 248   value_clear(tmp);
 249
 250   value_set_si(global_var_lcm, 1);
 251   for (i = 0; i < nb_vars; i++) {
 252     value_multiply(*det, *det, denoms->p[i]);
 253     Lcm3(global_var_lcm, denoms->p[i], &global_var_lcm);
 254   }
 255
 256   /* scale vertices */
 257   for (V = PP->V; V; V = V->next)
 258     for (i = 0; i < nb_vars; i++) {
 259       Vector_Scale(V->Vertex->p[i], V->Vertex->p[i], denoms->p[i], nb_param+1);
 260       Vector_Normalize(V->Vertex->p[i], nb_param+2);
 261     }
 262
 263   /* the expansion can be actually writen as global_var_lcm.L^{-1} */
 264   /* this is equivalent to multiply the rows of P by denoms_det */
 265   for (i = 0; i < nb_vars; i++)
 266     value_division(denoms->p[i], global_var_lcm, denoms->p[i]);
 267
 268   /* OPT : we could use a vector instead of a diagonal matrix here (c- and d-).*/
 269   /* c- make the quick expansion matrix */
 270   expansion = Matrix_Alloc(nb_vars+nb_param+1, nb_vars+nb_param+1);
 271   for (i = 0; i < nb_vars; i++)
 272     value_assign(expansion->p[i][i], denoms->p[i]);
 273   for (i = nb_vars; i < nb_vars+nb_param+1; i++)
 274     value_assign(expansion->p[i][i], global_var_lcm);
 275
 276   /* d- apply the variable expansion to the polyhedron */
 277   if (P)
 278     *P = Polyhedron_Preimage(*P, expansion, MaxRays);
 279
 280   if (Lat) {
 281     Lattice *L = Matrix_Alloc(nb_vars+1, nb_vars+1);
 282     for (i = 0; i < nb_vars; ++i)
 283         value_assign(L->p[i][i], denoms->p[i]);
 284     value_assign(L->p[nb_vars][nb_vars], global_var_lcm);
 285     *Lat = L;
 286   }
 287
 288   Matrix_Free(expansion);
 289   value_clear(global_var_lcm);
 290   Vector_Free(denoms);
 291 }
 292
 293 /* Compute negated sum of all positive or negative coefficients in a row */
 294 static void negated_sum(Value *v, int len, int negative, Value *sum)
 295 {
 296     int j;
 297
 298     value_set_si(*sum, 0);
 299     for (j = 0; j < len; ++j)
 300         if (negative ? value_neg_p(v[j]) : value_pos_p(v[j]))
 301             value_subtract(*sum, *sum, v[j]);
 302 }
 303
 304 /* adapted from mpolyhedron_inflate in PolyLib */
 305 Polyhedron *Polyhedron_Flate(Polyhedron *P, unsigned nparam, int inflate,
 306                              unsigned MaxRays)
 307 {
 308     Value sum;
 309     int nvar = P->Dimension - nparam;
 310     Matrix *C = Polyhedron2Constraints(P);
 311     Polyhedron *P2;
 312     int i;
 313
 314     value_init(sum);
 315     /* subtract the sum of the negative coefficients of each inequality */
 316     for (i = 0; i < C->NbRows; ++i) {
 317         negated_sum(C->p[i]+1, nvar, inflate, &sum);
 318         value_addto(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], sum);
 319     }
 320     value_clear(sum);
 321     P2 = Constraints2Polyhedron(C, MaxRays);
 322     Matrix_Free(C);
 323
 324     if (inflate) {
 325         Polyhedron *C, *CA;
 326         C = Polyhedron_Project(P, nparam);
 327         CA = align_context(C, P->Dimension, MaxRays);
 328         P = P2;
 329         P2 = DomainIntersection(P, CA, MaxRays);
 330         Polyhedron_Free(C);
 331         Polyhedron_Free(CA);
 332         Polyhedron_Free(P);
 333     }
 334
 335     return P2;
 336 }
 337
 338 static Polyhedron *flate_narrow2(Polyhedron *P, Lattice *L,
 339                                  unsigned nparam, int inflate,
 340                                  unsigned MaxRays)
 341 {
 342     Value sum;
 343     unsigned nvar = P->Dimension - nparam;
 344     Matrix *expansion;
 345     Matrix *C;
 346     int i, j;
 347     Polyhedron *P2;
 348
 349     expansion = Matrix_Alloc(nvar + nparam + 1,  nvar + nparam + 1);
 350     for (i = 0; i < nvar; ++i)
 351         Vector_Copy(L->p[i], expansion->p[i], nvar);
 352     for (i = nvar; i < nvar+nparam+1; ++i)
 353         value_assign(expansion->p[i][i], L->p[nvar][nvar]);
 354
 355     C = Matrix_Alloc(P->NbConstraints+1, 1+P->Dimension+1);
 356     value_init(sum);
 357     for (i = 0; i < P->NbConstraints; ++i) {
 358         negated_sum(P->Constraint[i]+1, nvar, inflate, &sum);
 359         value_assign(C->p[i][0], P->Constraint[i][0]);
 360         Vector_Matrix_Product(P->Constraint[i]+1, expansion, C->p[i]+1);
 361         if (value_zero_p(sum))
 362             continue;
 363         Vector_Copy(C->p[i]+1, C->p[i+1]+1, P->Dimension+1);
 364         value_addmul(C->p[i][1+P->Dimension], sum, L->p[nvar][nvar]);
 365         ConstraintSimplify(C->p[i], C->p[i], P->Dimension+2, &sum);
 366         ConstraintSimplify(C->p[i+1], C->p[i+1], P->Dimension+2, &sum);
 367         if (value_ne(C->p[i][1+P->Dimension], C->p[i+1][1+P->Dimension])) {
 368             if (inflate)
 369                 value_decrement(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension]);
 370             else
 371                 value_increment(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension]);
 372         }
 373     }
 374     value_clear(sum);
 375     C->NbRows--;
 376     P2 = Constraints2Polyhedron(C, MaxRays);
 377     Matrix_Free(C);
 378
 379     Matrix_Free(expansion);
 380
 381     if (inflate) {
 382         Polyhedron *C, *CA;
 383         C = Polyhedron_Project(P, nparam);
 384         CA = align_context(C, P->Dimension, MaxRays);
 385         P = P2;
 386         P2 = DomainIntersection(P, CA, MaxRays);
 387         Polyhedron_Free(C);
 388         Polyhedron_Free(CA);
 389         Polyhedron_Free(P);
 390     }
 391
 392     return P2;
 393 }
 394
 395 static void linear_min(Polyhedron *D, Value *obj, Value *min)
 396 {
 397     int i;
 398     Value tmp;
 399     value_init(tmp);
 400     POL_ENSURE_VERTICES(D);
 401     for (i = 0; i < D->NbRays; ++i) {
 402         Inner_Product(obj, D->Ray[i]+1, D->Dimension, &tmp);
 403         mpz_cdiv_q(tmp, tmp, D->Ray[i][1+D->Dimension]);
 404         if (!i || value_lt(tmp, *min))
 405             value_assign(*min, tmp);
 406     }
 407     value_clear(tmp);
 408 }
 409
 410 static Polyhedron *inflate_deflate_domain(Lattice *L, unsigned MaxRays)
 411 {
 412     unsigned nvar = L->NbRows-1;
 413     int i;
 414     Matrix *M;
 415     Polyhedron *D;
 416
 417     M = Matrix_Alloc(2*nvar, 1+nvar+1);
 418     for (i = 0; i < nvar; ++i) {
 419         value_set_si(M->p[2*i][0], 1);
 420         Vector_Copy(L->p[i], M->p[2*i]+1, nvar);
 421         Vector_Normalize(M->p[2*i]+1, nvar);
 422
 423         value_set_si(M->p[2*i+1][0], 1);
 424         Vector_Oppose(L->p[i], M->p[2*i+1]+1, nvar);
 425         value_assign(M->p[2*i+1][1+nvar], L->p[nvar][nvar]);
 426         Vector_Normalize(M->p[2*i+1]+1, nvar+1);
 427         value_decrement(M->p[2*i+1][1+nvar], M->p[2*i+1][1+nvar]);
 428     }
 429     D = Constraints2Polyhedron(M, MaxRays);
 430     Matrix_Free(M);
 431
 432     return D;
 433 }
 434
 435 static Polyhedron *flate_narrow(Polyhedron *P, Lattice *L,
 436                                 unsigned nparam, int inflate, unsigned MaxRays)
 437 {
 438     int i;
 439     unsigned nvar = P->Dimension - nparam;
 440     Vector *obj;
 441     Value min;
 442     Matrix *C;
 443     Polyhedron *D;
 444     Polyhedron *P2;
 445
 446     D = inflate_deflate_domain(L, MaxRays);
 447     value_init(min);
 448     obj = Vector_Alloc(nvar);
 449     C = Polyhedron2Constraints(P);
 450
 451     for (i = 0; i < C->NbRows; ++i) {
 452         if (inflate)
 453             Vector_Copy(C->p[i]+1, obj->p, nvar);
 454         else
 455             Vector_Oppose(C->p[i]+1, obj->p, nvar);
 456         linear_min(D, obj->p, &min);
 457         if (inflate)
 458             value_subtract(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], min);
 459         else
 460             value_addto(C->p[i][1+P->Dimension], C->p[i][1+P->Dimension], min);
 461     }
 462
 463     Polyhedron_Free(D);
 464     P2 = Constraints2Polyhedron(C, MaxRays);
 465     Matrix_Free(C);
 466     Vector_Free(obj);
 467     value_clear(min);
 468
 469     if (inflate) {
 470         Polyhedron *C, *CA;
 471         C = Polyhedron_Project(P, nparam);
 472         CA = align_context(C, P->Dimension, MaxRays);
 473         P = P2;
 474         P2 = DomainIntersection(P, CA, MaxRays);
 475         Polyhedron_Free(C);
 476         Polyhedron_Free(CA);
 477         Polyhedron_Free(P);
 478     }
 479
 480     return P2;
 481 }
 482
 483 static Polyhedron *flate(Polyhedron *P, Lattice *L,
 484                          unsigned nparam, int inflate,
 485                          struct barvinok_options *options)
 486 {
 487     if (options->scale_flags & BV_APPROX_SCALE_NARROW2)
 488         return flate_narrow2(P, L, nparam, inflate, options->MaxRays);
 489     else if (options->scale_flags & BV_APPROX_SCALE_NARROW)
 490         return flate_narrow(P, L, nparam, inflate, options->MaxRays);
 491     else
 492         return Polyhedron_Flate(P, nparam, inflate, options->MaxRays);
 493 }
 494
 495 static void Param_Polyhedron_Scale(Param_Polyhedron *PP, Polyhedron **P,
 496                                    Lattice **L,
 497                                    Value *det, struct barvinok_options *options)
 498 {
 499     if (options->scale_flags & BV_APPROX_SCALE_FAST)
 500         Param_Polyhedron_Scale_Integer_Fast(PP, P, L, det, options->MaxRays);
 501     else
 502         Param_Polyhedron_Scale_Integer_Slow(PP, P, L, det, options->MaxRays);
 503 }
 504
 505 static evalue *enumerate_flated(Polyhedron *P, Polyhedron *C, Lattice *L,
 506                                 struct barvinok_options *options)
 507 {
 508     unsigned nparam = C->Dimension;
 509     evalue *eres;
 510     int save_approximation = options->polynomial_approximation;
 511
 512     if (options->polynomial_approximation == BV_APPROX_SIGN_UPPER)
 513         P = flate(P, L, nparam, 1, options);
 514     if (options->polynomial_approximation == BV_APPROX_SIGN_LOWER)
 515         P = flate(P, L, nparam, 0, options);
 516
 517     /* Don't deflate/inflate again (on this polytope) */
 518     options->polynomial_approximation = BV_APPROX_SIGN_NONE;
 519     eres = barvinok_enumerate_with_options(P, C, options);
 520     options->polynomial_approximation = save_approximation;
 521     Polyhedron_Free(P);
 522
 523     return eres;
 524 }
 525
 526 static evalue *PP_enumerate_narrow_flated(Param_Polyhedron *PP,
 527                                           Polyhedron *P, Polyhedron *C,
 528                                           struct barvinok_options *options)
 529 {
 530     Polyhedron *Porig = P;
 531     int scale_narrow2 = options->scale_flags & BV_APPROX_SCALE_NARROW2;
 532     Lattice *L = NULL;
 533     Value det;
 534     evalue *eres;
 535
 536     value_init(det);
 537     value_set_si(det, 1);
 538
 539     Param_Polyhedron_Scale(PP, &P, &L, &det, options);
 540     Param_Polyhedron_Free(PP);
 541     if (scale_narrow2) {
 542         Polyhedron_Free(P);
 543         P = Porig;
 544     }
 545     /* Don't scale again (on this polytope) */
 546     options->approximation_method = BV_APPROX_NONE;
 547     eres = enumerate_flated(P, C, L, options);
 548     options->approximation_method = BV_APPROX_SCALE;
 549     Matrix_Free(L);
 550     if (P != Porig)
 551         Polyhedron_Free(P);
 552     if (value_notone_p(det))
 553         evalue_div(eres, det);
 554     value_clear(det);
 555     return eres;
 556 }
 557
 558 static Param_Polyhedron *Param_Polyhedron_Domain(Param_Polyhedron *PP,
 559                                                  Param_Domain *D,
 560                                                  Polyhedron *rVD)
 561 {
 562     int nv;
 563     Param_Polyhedron *PP_D;
 564     int i, ix;
 565     unsigned bx;
 566     Param_Vertices **next, *V;
 567
 568     PP_D = ALLOC(Param_Polyhedron);
 569     PP_D->D = ALLOC(Param_Domain);
 570     PP_D->D->next = NULL;
 571     PP_D->D->Domain = Domain_Copy(rVD);
 572     PP_D->V = NULL;
 573
 574     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 575     PP_D->D->F = ALLOCN(unsigned, nv);
 576     memset(PP_D->D->F, 0, nv * sizeof(unsigned));
 577
 578     next = &PP_D->V;
 579     i = 0;
 580     ix = 0;
 581     bx = MSB;
 582     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 583         Param_Vertices *V2 = ALLOC(Param_Vertices);
 584         V2->Vertex = Matrix_Copy(V->Vertex);
 585         V2->Domain = NULL;
 586         V2->next = NULL;
 587         *next = V2;
 588         next = &V2->next;
 589         PP_D->D->F[ix] |= bx;
 590         NEXT(ix, bx);
 591         ++i;
 592     END_FORALL_PVertex_in_ParamPolyhedron;
 593     PP_D->nbV = i;
 594
 595     return PP_D;
 596 }
 597
 598 static evalue *enumerate_narrow_flated(Polyhedron *P, Polyhedron *C,
 599                                         struct barvinok_options *options)
 600 {
 601     unsigned Rat_MaxRays = options->MaxRays;
 602     Param_Polyhedron *PP;
 603     PP = Polyhedron2Param_Polyhedron(P, C, options);
 604     POL_UNSET(Rat_MaxRays, POL_INTEGER);
 605
 606     if ((options->scale_flags & BV_APPROX_SCALE_CHAMBER) && PP->D->next) {
 607         int nd = -1;
 608         evalue *tmp, *eres = NULL;
 609         Polyhedron *TC = true_context(P, C, options->MaxRays);
 610
 611         FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
 612             Polyhedron *P2, *CA;
 613             Param_Polyhedron *PP_D;
 614             /* Intersect with D->Domain, so we only have the relevant constraints
 615              * left.  Don't use rVD, though, since we still want to recognize
 616              * the defining constraints of the parametric vertices.
 617              */
 618             CA = align_context(D->Domain, P->Dimension, options->MaxRays);
 619             P2 = DomainIntersection(P, CA, Rat_MaxRays);
 620             POL_ENSURE_VERTICES(P2);
 621             Polyhedron_Free(CA);
 622             /* Use rVD for context, to avoid overlapping domains in
 623              * results of computations in different chambers.
 624              */
 625             PP_D = Param_Polyhedron_Domain(PP, D, rVD);
 626             tmp = PP_enumerate_narrow_flated(PP_D, P2, rVD, options);
 627             Polyhedron_Free(P2);
 628             if (!eres)
 629                 eres = tmp;
 630             else {
 631                 eadd(tmp, eres);
 632                 free_evalue_refs(tmp);
 633                 free(tmp);
 634             }
 635             Polyhedron_Free(rVD);
 636         END_FORALL_REDUCED_DOMAIN
 637         Param_Polyhedron_Free(PP);
 638         if (!eres)
 639             eres = evalue_zero();
 640         Polyhedron_Free(TC);
 641         return eres;
 642     } else
 643         return PP_enumerate_narrow_flated(PP, P, C, options);
 644 }
 645
 646 /* If scaling is to be performed in combination with deflation/inflation,
 647  * do both and return the result.
 648  * Otherwise return NULL.
 649  */
 650 evalue *scale_bound(Polyhedron *P, Polyhedron *C,
 651                     struct barvinok_options *options)
 652 {
 653     int scale_narrow = options->scale_flags & BV_APPROX_SCALE_NARROW;
 654     int scale_narrow2 = options->scale_flags & BV_APPROX_SCALE_NARROW2;
 655
 656     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE ||
 657         options->polynomial_approximation == BV_APPROX_SIGN_APPROX)
 658         return NULL;
 659
 660     if (scale_narrow || scale_narrow2)
 661         return enumerate_narrow_flated(P, C, options);
 662     else
 663         return enumerate_flated(P, C, NULL, options);
 664 }
 665
 666 evalue *scale(Param_Polyhedron *PP, Polyhedron *P, Polyhedron *C,
 667               struct barvinok_options *options)
 668 {
 669     Polyhedron *T = P;
 670     unsigned MaxRays;
 671     evalue *eres = NULL;
 672     Value det;
 673
 674     if ((options->scale_flags & BV_APPROX_SCALE_CHAMBER) && PP->D->next) {
 675         int nd = -1;
 676         evalue *tmp;
 677         Polyhedron *TC = true_context(P, C, options->MaxRays);
 678
 679         FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
 680             Param_Polyhedron *PP_D = Param_Polyhedron_Domain(PP, D, rVD);
 681             tmp = scale(PP_D, P, rVD, options);
 682             if (!eres)
 683                 eres = tmp;
 684             else {
 685                 eadd(tmp, eres);
 686                 free_evalue_refs(tmp);
 687                 free(tmp);
 688             }
 689             Param_Polyhedron_Free(PP_D);
 690             Polyhedron_Free(rVD);
 691         END_FORALL_REDUCED_DOMAIN
 692         if (!eres)
 693             eres = evalue_zero();
 694         Polyhedron_Free(TC);
 695         return eres;
 696     }
 697
 698     value_init(det);
 699     value_set_si(det, 1);
 700
 701     MaxRays = options->MaxRays;
 702     POL_UNSET(options->MaxRays, POL_INTEGER);
 703     Param_Polyhedron_Scale(PP, &T, NULL, &det, options);
 704     options->MaxRays = MaxRays;
 705
 706     eres = Param_Polyhedron_Enumerate(PP, T, C, options);
 707     if (P != T)
 708         Polyhedron_Free(T);
 709
 710     if (value_notone_p(det))
 711         evalue_div(eres, det);
 712     value_clear(det);
 713
 714     return eres;
 715 }