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