scale.c

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