volume.c

   1 #include <barvinok/polylib.h>
   2 #include <barvinok/barvinok.h>
   3 #include <barvinok/options.h>
   4 #include <barvinok/util.h>
   5 #include "reduce_domain.h"
   6 #include "scale.h"
   7 #include "volume.h"
   8
   9 #define ALLOC(type) (type*)malloc(sizeof(type))
  10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  11
  12 /* Computes an evalue representation of a coordinate
  13  * in a ParamVertices.
  14  */
  15 static evalue *vertex2evalue(Value *vertex, int nparam)
  16 {
  17     return affine2evalue(vertex, vertex[nparam+1], nparam);
  18 }
  19
  20 static void matrix_print(evalue ***matrix, int dim, int *cols,
  21                          char **param_names)
  22 {
  23     int i, j;
  24
  25     for (i = 0; i < dim; ++i)
  26         for (j = 0; j < dim; ++j) {
  27             int k = cols ? cols[j] : j;
  28             fprintf(stderr, "%d %d: ", i, j);
  29             print_evalue(stderr, matrix[i][k], param_names);
  30             fprintf(stderr, "\n");
  31         }
  32 }
  33
  34 /* Compute determinant using Laplace's formula.
  35  * In particular, the determinant is expanded along the last row.
  36  * The cols array is a list of columns that remain in the currect submatrix.
  37  */
  38 static evalue *determinant_cols(evalue ***matrix, int dim, int *cols)
  39 {
  40     evalue *det, *tmp;
  41     evalue mone;
  42
  43     if (dim == 1) {
  44         det = ALLOC(evalue);
  45         value_init(det->d);
  46         evalue_copy(det, matrix[0][cols[0]]);
  47         return det;
  48     }
  49
  50     value_init(mone.d);
  51     evalue_set_si(&mone, -1, 1);
  52     int j;
  53     det = NULL;
  54     int *newcols = ALLOCN(int, dim-1);
  55     for (j = 1; j < dim; ++j)
  56         newcols[j-1] = cols[j];
  57     for (j = 0; j < dim; ++j) {
  58         if (j != 0)
  59             newcols[j-1] = cols[j-1];
  60         tmp = determinant_cols(matrix, dim-1, newcols);
  61         emul(matrix[dim-1][cols[j]], tmp);
  62         if ((dim+j)%2 == 0)
  63             emul(&mone, tmp);
  64         if (!det)
  65             det = tmp;
  66         else {
  67             eadd(tmp, det);
  68             free_evalue_refs(tmp);
  69             free(tmp);
  70         }
  71     }
  72     free(newcols);
  73     free_evalue_refs(&mone);
  74
  75     return det;
  76 }
  77
  78 static evalue *determinant(evalue ***matrix, int dim)
  79 {
  80     int i;
  81     int *cols = ALLOCN(int, dim);
  82     evalue *det;
  83
  84     for (i = 0; i < dim; ++i)
  85         cols[i] = i;
  86
  87     det = determinant_cols(matrix, dim, cols);
  88
  89     free(cols);
  90
  91     return det;
  92 }
  93
  94 /* Compute the facet of P that saturates constraint c.
  95  */
  96 static Polyhedron *facet(Polyhedron *P, int c, unsigned MaxRays)
  97 {
  98     Polyhedron *F;
  99     Vector *row = Vector_Alloc(1+P->Dimension+1);
 100     Vector_Copy(P->Constraint[c]+1, row->p+1, P->Dimension+1);
 101     F = AddConstraints(row->p, 1, P, MaxRays);
 102     Vector_Free(row);
 103     return F;
 104 }
 105
 106 /* Compute a dummy Param_Domain that contains all vertices of Param_Domain D
 107  * (which contains the vertices of P) that lie on the facet obtain by
 108  * saturating constraint c of P
 109  */
 110 static Param_Domain *face_vertices(Param_Polyhedron *PP, Param_Domain *D,
 111                                    Polyhedron *P, int c)
 112 {
 113     int nv;
 114     Param_Vertices *V;
 115     Param_Domain *FD = ALLOC(Param_Domain);
 116     FD->Domain = 0;
 117     FD->next = 0;
 118
 119     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 120     FD->F = ALLOCN(unsigned, nv);
 121     memset(FD->F, 0, nv * sizeof(unsigned));
 122
 123     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _i, _ix, _bx internal counters */
 124         int n;
 125         unsigned char *supporting = supporting_constraints(P, V, &n);
 126         if (supporting[c])
 127             FD->F[_ix] |= _bx;
 128         free(supporting);
 129     END_FORALL_PVertex_in_ParamPolyhedron;
 130
 131     return FD;
 132 }
 133
 134 /* Substitute parameters by the corresponding element in subs
 135  */
 136 static evalue *evalue_substitute_new(evalue *e, evalue **subs)
 137 {
 138     evalue *res = NULL;
 139     evalue *c;
 140     int i;
 141
 142     if (value_notzero_p(e->d)) {
 143         res = ALLOC(evalue);
 144         value_init(res->d);
 145         evalue_copy(res, e);
 146         return res;
 147     }
 148     assert(e->x.p->type == polynomial);
 149
 150     res = evalue_zero();
 151     for (i = e->x.p->size-1; i > 0; --i) {
 152         c = evalue_substitute_new(&e->x.p->arr[i], subs);
 153         eadd(c, res);
 154         free_evalue_refs(c);
 155         free(c);
 156         emul(subs[e->x.p->pos-1], res);
 157     }
 158     c = evalue_substitute_new(&e->x.p->arr[0], subs);
 159     eadd(c, res);
 160     free_evalue_refs(c);
 161     free(c);
 162
 163     return res;
 164 }
 165
 166 /* Plug in the parametric vertex V in the constraint constraint.
 167  * The result is stored in row, with the denominator in position 0.
 168  */
 169 static void Param_Inner_Product(Value *constraint, Matrix *Vertex,
 170                                 Value *row)
 171 {
 172     unsigned nparam = Vertex->NbColumns - 2;
 173     unsigned nvar = Vertex->NbRows;
 174     int j;
 175     Value tmp, tmp2;
 176
 177     value_set_si(row[0], 1);
 178     Vector_Set(row+1, 0, nparam+1);
 179
 180     value_init(tmp);
 181     value_init(tmp2);
 182
 183     for (j = 0 ; j < nvar; ++j) {
 184         value_set_si(tmp, 1);
 185         value_assign(tmp2,  constraint[1+j]);
 186         if (value_ne(row[0], Vertex->p[j][nparam+1])) {
 187             value_assign(tmp, row[0]);
 188             value_lcm(row[0], Vertex->p[j][nparam+1], &row[0]);
 189             value_division(tmp, row[0], tmp);
 190             value_multiply(tmp2, tmp2, row[0]);
 191             value_division(tmp2, tmp2, Vertex->p[j][nparam+1]);
 192         }
 193         Vector_Combine(row+1, Vertex->p[j], row+1, tmp, tmp2, nparam+1);
 194     }
 195     value_set_si(tmp, 1);
 196     Vector_Combine(row+1, constraint+1+nvar, row+1, tmp, row[0], nparam+1);
 197
 198     value_clear(tmp);
 199     value_clear(tmp2);
 200 }
 201
 202 struct parameter_point {
 203     Vector *coord;
 204     evalue **e;
 205 };
 206
 207 struct parameter_point *parameter_point_new(unsigned nparam)
 208 {
 209     struct parameter_point *point = ALLOC(struct parameter_point);
 210     point->coord = Vector_Alloc(nparam+1);
 211     point->e = NULL;
 212     return point;
 213 }
 214
 215 evalue **parameter_point_evalue(struct parameter_point *point)
 216 {
 217     int j;
 218     unsigned nparam = point->coord->Size-1;
 219
 220     if (point->e)
 221         return point->e;
 222
 223     point->e = ALLOCN(evalue *, nparam);
 224     for (j = 0; j < nparam; ++j) {
 225         point->e[j] = ALLOC(evalue);
 226         value_init(point->e[j]->d);
 227         evalue_set(point->e[j], point->coord->p[j], point->coord->p[nparam]);
 228     }
 229
 230     return point->e;
 231 }
 232
 233 void parameter_point_free(struct parameter_point *point)
 234 {
 235     int i;
 236     unsigned nparam = point->coord->Size-1;
 237
 238     Vector_Free(point->coord);
 239
 240     if (point->e) {
 241         for (i = 0; i < nparam; ++i) {
 242             free_evalue_refs(point->e[i]);
 243             free(point->e[i]);
 244         }
 245         free(point->e);
 246     }
 247     free(point);
 248 }
 249
 250 /* Computes point in pameter space where polyhedron is non-empty.
 251  */
 252 static struct parameter_point *non_empty_point(Param_Domain *D)
 253 {
 254     unsigned nparam = D->Domain->Dimension;
 255     struct parameter_point *point;
 256     Vector *v;
 257
 258     v = inner_point(D->Domain);
 259     point = parameter_point_new(nparam);
 260     Vector_Copy(v->p+1, point->coord->p, nparam+1);
 261     Vector_Free(v);
 262
 263     return point;
 264 }
 265
 266 static Matrix *barycenter(Param_Polyhedron *PP, Param_Domain *D)
 267 {
 268     int nbV;
 269     Matrix *center = NULL;
 270     Value denom;
 271     Value fc, fv;
 272     unsigned nparam;
 273     int i;
 274     Param_Vertices *V;
 275
 276     value_init(fc);
 277     value_init(fv);
 278     nbV = 0;
 279     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 280         ++nbV;
 281         if (!center) {
 282             center = Matrix_Copy(V->Vertex);
 283             nparam = center->NbColumns - 2;
 284         } else {
 285             for (i = 0; i < center->NbRows; ++i) {
 286                 value_assign(fc, center->p[i][nparam+1]);
 287                 value_lcm(fc, V->Vertex->p[i][nparam+1],
 288                             &center->p[i][nparam+1]);
 289                 value_division(fc, center->p[i][nparam+1], fc);
 290                 value_division(fv, center->p[i][nparam+1],
 291                                 V->Vertex->p[i][nparam+1]);
 292                 Vector_Combine(center->p[i], V->Vertex->p[i], center->p[i],
 293                                fc, fv, nparam+1);
 294             }
 295         }
 296     END_FORALL_PVertex_in_ParamPolyhedron;
 297     value_clear(fc);
 298     value_clear(fv);
 299
 300     value_init(denom);
 301     value_set_si(denom, nbV);
 302     for (i = 0; i < center->NbRows; ++i) {
 303         value_multiply(center->p[i][nparam+1], center->p[i][nparam+1], denom);
 304         Vector_Normalize(center->p[i], nparam+2);
 305     }
 306     value_clear(denom);
 307
 308     return center;
 309 }
 310
 311 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
 312  * If F is a simplex, then the volume is computed of a recursive pyramid
 313  * over F with the points already in matrix.
 314  * Otherwise, the barycenter of F is added to matrix and the function
 315  * is called recursively on the facets of F.
 316  *
 317  * The first row of matrix contain the _negative_ of the first point.
 318  * The remaining rows of matrix contain the distance of the corresponding
 319  * point to the first point.
 320  */
 321 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 322                                 unsigned dim, evalue ***matrix,
 323                                 struct parameter_point *point,
 324                                 int row, Polyhedron *F,
 325                                 struct barvinok_options *options);
 326
 327 static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 328                                   unsigned dim, evalue ***matrix,
 329                                   struct parameter_point *point,
 330                                   int row, Polyhedron *F,
 331                                   struct barvinok_options *options)
 332 {
 333     int j;
 334     evalue *tmp;
 335     evalue *vol;
 336     evalue mone;
 337     Matrix *center;
 338     unsigned cut_MaxRays = options->MaxRays;
 339     unsigned nparam = D->Domain->Dimension;
 340
 341     POL_UNSET(cut_MaxRays, POL_INTEGER);
 342
 343     value_init(mone.d);
 344     evalue_set_si(&mone, -1, 1);
 345
 346     center = barycenter(PP, D);
 347     for (j = 0; j < dim; ++j)
 348         matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
 349
 350     if (row == 0) {
 351         for (j = 0; j < dim; ++j)
 352             emul(&mone, matrix[row][j]);
 353     } else {
 354         for (j = 0; j < dim; ++j)
 355             eadd(matrix[0][j], matrix[row][j]);
 356     }
 357
 358     vol = NULL;
 359     POL_ENSURE_FACETS(F);
 360     for (j = F->NbEq; j < F->NbConstraints; ++j) {
 361         Polyhedron *FF;
 362         Param_Domain *FD;
 363         if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 364             continue;
 365         FF = facet(F, j, options->MaxRays);
 366         FD = face_vertices(PP, D, F, j);
 367         tmp = volume_in_domain(PP, FD, dim, matrix, point,
 368                                row+1, FF, options);
 369         if (!vol)
 370             vol = tmp;
 371         else {
 372             eadd(tmp, vol);
 373             free_evalue_refs(tmp);
 374             free(tmp);
 375         }
 376         Polyhedron_Free(FF);
 377         Param_Domain_Free(FD);
 378     }
 379
 380     Matrix_Free(center);
 381
 382     for (j = 0; j < dim; ++j) {
 383         free_evalue_refs(matrix[row][j]);
 384         free(matrix[row][j]);
 385     }
 386
 387     free_evalue_refs(&mone);
 388     return vol;
 389 }
 390
 391 static evalue *volume_simplex(Param_Polyhedron *PP, Param_Domain *D,
 392                                 unsigned dim, evalue ***matrix,
 393                                 struct parameter_point *point,
 394                                 int row, unsigned MaxRays)
 395 {
 396     evalue mone;
 397     Param_Vertices *V;
 398     evalue *vol, *val;
 399     int i, j;
 400
 401     value_init(mone.d);
 402     evalue_set_si(&mone, -1, 1);
 403
 404     i = row;
 405     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal counters */
 406         for (j = 0; j < dim; ++j) {
 407             matrix[i][j] = vertex2evalue(V->Vertex->p[j],
 408                                            V->Vertex->NbColumns - 2);
 409             if (i == 0)
 410                 emul(&mone, matrix[i][j]);
 411             else
 412                 eadd(matrix[0][j], matrix[i][j]);
 413         }
 414         ++i;
 415     END_FORALL_PVertex_in_ParamPolyhedron;
 416
 417     vol = determinant(matrix+1, dim);
 418
 419     val = evalue_substitute_new(vol, parameter_point_evalue(point));
 420
 421     assert(value_notzero_p(val->d));
 422     assert(value_notzero_p(val->x.n));
 423     if (value_neg_p(val->x.n))
 424         emul(&mone, vol);
 425
 426     free_evalue_refs(val);
 427     free(val);
 428
 429     for (i = row; i < dim+1; ++i) {
 430         for (j = 0; j < dim; ++j) {
 431             free_evalue_refs(matrix[i][j]);
 432             free(matrix[i][j]);
 433         }
 434     }
 435
 436     free_evalue_refs(&mone);
 437
 438     return vol;
 439 }
 440
 441 static evalue *volume_triangulate_lift(Param_Polyhedron *PP, Param_Domain *D,
 442                                         unsigned dim, evalue ***matrix,
 443                                         struct parameter_point *point,
 444                                         int row, unsigned MaxRays)
 445 {
 446     const static int MAX_TRY=10;
 447     Param_Vertices *V;
 448     int nbV, nv;
 449     int i;
 450     int t = 0;
 451     Matrix *FixedRays, *M;
 452     Polyhedron *L;
 453     Param_Domain SD;
 454     Value tmp;
 455     evalue *s, *vol;
 456
 457     SD.Domain = 0;
 458     SD.next = 0;
 459     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 460     SD.F = ALLOCN(unsigned, nv);
 461
 462     FixedRays = Matrix_Alloc(PP->nbV+1, 1+dim+2);
 463     nbV = 0;
 464     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 465         unsigned nparam = V->Vertex->NbColumns - 2;
 466         Param_Vertex_Common_Denominator(V);
 467         for (i = 0; i < V->Vertex->NbRows; ++i) {
 468             value_multiply(FixedRays->p[nbV][1+i], V->Vertex->p[i][nparam],
 469                            point->coord->p[nparam]);
 470             Inner_Product(V->Vertex->p[i], point->coord->p, nparam,
 471                           &FixedRays->p[nbV][1+dim]);
 472             value_addto(FixedRays->p[nbV][1+i], FixedRays->p[nbV][1+i],
 473                         FixedRays->p[nbV][1+dim]);
 474         }
 475         value_multiply(FixedRays->p[nbV][1+dim+1], V->Vertex->p[0][nparam+1],
 476                        point->coord->p[nparam]);
 477         value_set_si(FixedRays->p[nbV][0], 1);
 478         ++nbV;
 479     END_FORALL_PVertex_in_ParamPolyhedron;
 480     value_set_si(FixedRays->p[nbV][0], 1);
 481     value_set_si(FixedRays->p[nbV][1+dim], 1);
 482     FixedRays->NbRows = nbV+1;
 483
 484     value_init(tmp);
 485     if (0) {
 486 try_again:
 487         /* Usually vol should still be NULL */
 488         if (vol) {
 489             free_evalue_refs(vol);
 490             free(vol);
 491         }
 492         Polyhedron_Free(L);
 493         ++t;
 494     }
 495     assert(t <= MAX_TRY);
 496     vol = NULL;
 497
 498     for (i = 0; i < nbV; ++i)
 499         value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
 500
 501     M = Matrix_Copy(FixedRays);
 502     L = Rays2Polyhedron(M, MaxRays);
 503     Matrix_Free(M);
 504
 505     POL_ENSURE_FACETS(L);
 506     for (i = 0; i < L->NbConstraints; ++i) {
 507         int r;
 508         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 509         if (value_negz_p(L->Constraint[i][1+dim]))
 510             continue;
 511
 512         memset(SD.F, 0, nv * sizeof(unsigned));
 513         nbV = 0;
 514         r = 0;
 515         FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal */
 516             Inner_Product(FixedRays->p[nbV]+1, L->Constraint[i]+1, dim+2, &tmp);
 517             if (value_zero_p(tmp)) {
 518                 if (r > dim-row)
 519                     goto try_again;
 520                 SD.F[_ix] |= _bx;
 521                 ++r;
 522             }
 523             ++nbV;
 524         END_FORALL_PVertex_in_ParamPolyhedron;
 525         assert(r == (dim-row)+1);
 526
 527         s = volume_simplex(PP, &SD, dim, matrix, point, row, MaxRays);
 528         if (!vol)
 529             vol = s;
 530         else {
 531             eadd(s, vol);
 532             free_evalue_refs(s);
 533             free(s);
 534         }
 535     }
 536     Polyhedron_Free(L);
 537     Matrix_Free(FixedRays);
 538     value_clear(tmp);
 539     free(SD.F);
 540
 541     return vol;
 542 }
 543
 544 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 545                                 unsigned dim, evalue ***matrix,
 546                                 struct parameter_point *point,
 547                                 int row, Polyhedron *F,
 548                                 struct barvinok_options *options)
 549 {
 550     int nbV;
 551     Param_Vertices *V;
 552     evalue *vol;
 553
 554     assert(point);
 555
 556     nbV = 0;
 557     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 558         ++nbV;
 559     END_FORALL_PVertex_in_ParamPolyhedron;
 560
 561     if (nbV > (dim-row) + 1) {
 562         if (options->volume_triangulate_lift)
 563             vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 564                                           row, options->MaxRays);
 565         else
 566             vol = volume_triangulate(PP, D, dim, matrix, point,
 567                                      row, F, options);
 568     } else {
 569         assert(nbV == (dim-row) + 1);
 570         vol = volume_simplex(PP, D, dim, matrix, point, row, options->MaxRays);
 571     }
 572
 573     return vol;
 574 }
 575
 576 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
 577                                 struct barvinok_options *options)
 578 {
 579     evalue ***matrix;
 580     unsigned nparam = C->Dimension;
 581     unsigned nvar = P->Dimension - C->Dimension;
 582     Param_Polyhedron *PP;
 583     unsigned PP_MaxRays = options->MaxRays;
 584     unsigned MaxRays;
 585     int i, j;
 586     Value fact;
 587     evalue *vol;
 588     int nd;
 589     struct section { Polyhedron *D; evalue *E; } *s;
 590     Param_Domain *D;
 591     Polyhedron *TC;
 592
 593     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE)
 594         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 595
 596     if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
 597         int pa = options->polynomial_approximation;
 598         assert(pa == BV_APPROX_SIGN_UPPER || pa == BV_APPROX_SIGN_LOWER);
 599
 600         P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
 601                              options->MaxRays);
 602
 603         /* Don't deflate/inflate again (on this polytope) */
 604         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 605         vol = barvinok_enumerate_with_options(P, C, options);
 606         options->polynomial_approximation = pa;
 607
 608         Polyhedron_Free(P);
 609         return vol;
 610     }
 611
 612     TC = true_context(P, NULL, C, options->MaxRays);
 613
 614     if (PP_MaxRays & POL_NO_DUAL)
 615         PP_MaxRays = 0;
 616
 617     MaxRays = options->MaxRays;
 618     POL_UNSET(options->MaxRays, POL_INTEGER);
 619
 620     value_init(fact);
 621     Factorial(nvar, &fact);
 622
 623     PP = Polyhedron2Param_Domain(P, C, PP_MaxRays);
 624
 625     for (nd = 0, D = PP->D; D; ++nd, D = D->next);
 626     s = ALLOCN(struct section, nd);
 627
 628     matrix = ALLOCN(evalue **, nvar+1);
 629     for (i = 0; i < nvar+1; ++i)
 630         matrix[i] = ALLOCN(evalue *, nvar);
 631
 632     FORALL_REDUCED_DOMAIN(PP, TC, NULL, NULL, nd, options, i, D, rVD)
 633         Polyhedron *CA, *F;
 634         struct parameter_point *point;
 635
 636         CA = align_context(D->Domain, P->Dimension, MaxRays);
 637         F = DomainIntersection(P, CA, options->MaxRays);
 638         Domain_Free(CA);
 639
 640         point = non_empty_point(D);
 641         s[i].D = rVD;
 642         s[i].E = volume_in_domain(PP, D, nvar, matrix, point, 0, F, options);
 643         Domain_Free(F);
 644         parameter_point_free(point);
 645         evalue_div(s[i].E, fact);
 646     END_FORALL_REDUCED_DOMAIN
 647     options->MaxRays = MaxRays;
 648     Polyhedron_Free(TC);
 649
 650     vol = ALLOC(evalue);
 651     value_init(vol->d);
 652     value_set_si(vol->d, 0);
 653
 654     if (nd == 0)
 655         evalue_set_si(vol, 0, 1);
 656     else {
 657         vol->x.p = new_enode(partition, 2*nd, C->Dimension);
 658         for (i = 0; i < nd; ++i) {
 659             EVALUE_SET_DOMAIN(vol->x.p->arr[2*i], s[i].D);
 660             value_clear(vol->x.p->arr[2*i+1].d);
 661             vol->x.p->arr[2*i+1] = *s[i].E;
 662             free(s[i].E);
 663         }
 664     }
 665     free(s);
 666
 667     for (i = 0; i < nvar+1; ++i)
 668         free(matrix[i]);
 669     free(matrix);
 670     Param_Polyhedron_Free(PP);
 671     value_clear(fact);
 672
 673     return vol;
 674 }