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