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