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