scale.c

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