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