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 static Matrix *triangulation_vertex(Param_Polyhedron *PP, Param_Domain *D,
 315                                     Polyhedron *F)
 316 {
 317     Param_Vertices *V;
 318
 319     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 320         return V->Vertex;
 321     END_FORALL_PVertex_in_ParamPolyhedron;
 322
 323     assert(0);
 324     return NULL;
 325 }
 326
 327 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
 328  * If F is a simplex, then the volume is computed of a recursive pyramid
 329  * over F with the points already in matrix.
 330  * Otherwise, the barycenter of F is added to matrix and the function
 331  * is called recursively on the facets of F.
 332  *
 333  * The first row of matrix contain the _negative_ of the first point.
 334  * The remaining rows of matrix contain the distance of the corresponding
 335  * point to the first point.
 336  */
 337 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 338                                 unsigned dim, evalue ***matrix,
 339                                 struct parameter_point *point,
 340                                 int row, Polyhedron *F,
 341                                 struct barvinok_options *options);
 342
 343 static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 344                                   unsigned dim, evalue ***matrix,
 345                                   struct parameter_point *point,
 346                                   int row, Polyhedron *F,
 347                                   struct barvinok_options *options)
 348 {
 349     int j;
 350     evalue *tmp;
 351     evalue *vol;
 352     evalue mone;
 353     Matrix *center;
 354     unsigned cut_MaxRays = options->MaxRays;
 355     unsigned nparam = PP->V->Vertex->NbColumns-2;
 356     Vector *v = NULL;
 357
 358     POL_UNSET(cut_MaxRays, POL_INTEGER);
 359
 360     value_init(mone.d);
 361     evalue_set_si(&mone, -1, 1);
 362
 363     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 364         center = barycenter(PP, D);
 365     else
 366         center = triangulation_vertex(PP, D, F);
 367     for (j = 0; j < dim; ++j)
 368         matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
 369     if (options->volume_triangulate == BV_VOL_BARYCENTER)
 370         Matrix_Free(center);
 371     else
 372         v = Vector_Alloc(1+nparam+1);
 373
 374     if (row == 0) {
 375         for (j = 0; j < dim; ++j)
 376             emul(&mone, matrix[row][j]);
 377     } else {
 378         for (j = 0; j < dim; ++j)
 379             eadd(matrix[0][j], matrix[row][j]);
 380     }
 381
 382     vol = NULL;
 383     POL_ENSURE_FACETS(F);
 384     for (j = F->NbEq; j < F->NbConstraints; ++j) {
 385         Polyhedron *FF;
 386         Param_Domain *FD;
 387         if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 388             continue;
 389         if (options->volume_triangulate != BV_VOL_BARYCENTER) {
 390             Param_Inner_Product(F->Constraint[j], center, v->p);
 391             if (First_Non_Zero(v->p+1, nparam+1) == -1)
 392                 continue;
 393         }
 394         FF = facet(F, j, options->MaxRays);
 395         FD = face_vertices(PP, D, F, j);
 396         tmp = volume_in_domain(PP, FD, dim, matrix, point,
 397                                row+1, FF, options);
 398         if (!vol)
 399             vol = tmp;
 400         else {
 401             eadd(tmp, vol);
 402             free_evalue_refs(tmp);
 403             free(tmp);
 404         }
 405         Polyhedron_Free(FF);
 406         Param_Domain_Free(FD);
 407     }
 408
 409     if (options->volume_triangulate != BV_VOL_BARYCENTER)
 410         Vector_Free(v);
 411
 412     for (j = 0; j < dim; ++j) {
 413         free_evalue_refs(matrix[row][j]);
 414         free(matrix[row][j]);
 415     }
 416
 417     free_evalue_refs(&mone);
 418     return vol;
 419 }
 420
 421 static evalue *volume_simplex(Param_Polyhedron *PP, Param_Domain *D,
 422                                 unsigned dim, evalue ***matrix,
 423                                 struct parameter_point *point,
 424                                 int row, struct barvinok_options *options)
 425 {
 426     evalue mone;
 427     Param_Vertices *V;
 428     evalue *vol, *val;
 429     int i, j;
 430
 431     options->stats->volume_simplices++;
 432
 433     value_init(mone.d);
 434     evalue_set_si(&mone, -1, 1);
 435
 436     i = row;
 437     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal counters */
 438         for (j = 0; j < dim; ++j) {
 439             matrix[i][j] = vertex2evalue(V->Vertex->p[j],
 440                                            V->Vertex->NbColumns - 2);
 441             if (i == 0)
 442                 emul(&mone, matrix[i][j]);
 443             else
 444                 eadd(matrix[0][j], matrix[i][j]);
 445         }
 446         ++i;
 447     END_FORALL_PVertex_in_ParamPolyhedron;
 448
 449     vol = determinant(matrix+1, dim);
 450
 451     val = evalue_substitute_new(vol, parameter_point_evalue(point));
 452
 453     assert(value_notzero_p(val->d));
 454     assert(value_notzero_p(val->x.n));
 455     if (value_neg_p(val->x.n))
 456         emul(&mone, vol);
 457
 458     free_evalue_refs(val);
 459     free(val);
 460
 461     for (i = row; i < dim+1; ++i) {
 462         for (j = 0; j < dim; ++j) {
 463             free_evalue_refs(matrix[i][j]);
 464             free(matrix[i][j]);
 465         }
 466     }
 467
 468     free_evalue_refs(&mone);
 469
 470     return vol;
 471 }
 472
 473 static evalue *volume_triangulate_lift(Param_Polyhedron *PP, Param_Domain *D,
 474                                         unsigned dim, evalue ***matrix,
 475                                         struct parameter_point *point,
 476                                         int row, struct barvinok_options *options)
 477 {
 478     const static int MAX_TRY=10;
 479     Param_Vertices *V;
 480     int nbV, nv;
 481     int i;
 482     int t = 0;
 483     Matrix *FixedRays, *M;
 484     Polyhedron *L;
 485     Param_Domain SD;
 486     Value tmp;
 487     evalue *s, *vol;
 488
 489     SD.Domain = 0;
 490     SD.next = 0;
 491     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 492     SD.F = ALLOCN(unsigned, nv);
 493
 494     FixedRays = Matrix_Alloc(PP->nbV+1, 1+dim+2);
 495     nbV = 0;
 496     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 497         unsigned nparam = V->Vertex->NbColumns - 2;
 498         Param_Vertex_Common_Denominator(V);
 499         for (i = 0; i < V->Vertex->NbRows; ++i) {
 500             value_multiply(FixedRays->p[nbV][1+i], V->Vertex->p[i][nparam],
 501                            point->coord->p[nparam]);
 502             Inner_Product(V->Vertex->p[i], point->coord->p, nparam,
 503                           &FixedRays->p[nbV][1+dim]);
 504             value_addto(FixedRays->p[nbV][1+i], FixedRays->p[nbV][1+i],
 505                         FixedRays->p[nbV][1+dim]);
 506         }
 507         value_multiply(FixedRays->p[nbV][1+dim+1], V->Vertex->p[0][nparam+1],
 508                        point->coord->p[nparam]);
 509         value_set_si(FixedRays->p[nbV][0], 1);
 510         ++nbV;
 511     END_FORALL_PVertex_in_ParamPolyhedron;
 512     value_set_si(FixedRays->p[nbV][0], 1);
 513     value_set_si(FixedRays->p[nbV][1+dim], 1);
 514     FixedRays->NbRows = nbV+1;
 515
 516     value_init(tmp);
 517     if (0) {
 518 try_again:
 519         /* Usually vol should still be NULL */
 520         if (vol) {
 521             free_evalue_refs(vol);
 522             free(vol);
 523         }
 524         Polyhedron_Free(L);
 525         ++t;
 526     }
 527     assert(t <= MAX_TRY);
 528     vol = NULL;
 529
 530     for (i = 0; i < nbV; ++i)
 531         value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
 532
 533     M = Matrix_Copy(FixedRays);
 534     L = Rays2Polyhedron(M, options->MaxRays);
 535     Matrix_Free(M);
 536
 537     POL_ENSURE_FACETS(L);
 538     for (i = 0; i < L->NbConstraints; ++i) {
 539         int r;
 540         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 541         if (value_negz_p(L->Constraint[i][1+dim]))
 542             continue;
 543
 544         memset(SD.F, 0, nv * sizeof(unsigned));
 545         nbV = 0;
 546         r = 0;
 547         FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal */
 548             Inner_Product(FixedRays->p[nbV]+1, L->Constraint[i]+1, dim+2, &tmp);
 549             if (value_zero_p(tmp)) {
 550                 if (r > dim-row)
 551                     goto try_again;
 552                 SD.F[_ix] |= _bx;
 553                 ++r;
 554             }
 555             ++nbV;
 556         END_FORALL_PVertex_in_ParamPolyhedron;
 557         assert(r == (dim-row)+1);
 558
 559         s = volume_simplex(PP, &SD, dim, matrix, point, row, options);
 560         if (!vol)
 561             vol = s;
 562         else {
 563             eadd(s, vol);
 564             free_evalue_refs(s);
 565             free(s);
 566         }
 567     }
 568     Polyhedron_Free(L);
 569     Matrix_Free(FixedRays);
 570     value_clear(tmp);
 571     free(SD.F);
 572
 573     return vol;
 574 }
 575
 576 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 577                                 unsigned dim, evalue ***matrix,
 578                                 struct parameter_point *point,
 579                                 int row, Polyhedron *F,
 580                                 struct barvinok_options *options)
 581 {
 582     int nbV;
 583     Param_Vertices *V;
 584     evalue *vol;
 585
 586     assert(point);
 587
 588     nbV = 0;
 589     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 590         ++nbV;
 591     END_FORALL_PVertex_in_ParamPolyhedron;
 592
 593     if (nbV > (dim-row) + 1) {
 594         if (options->volume_triangulate == BV_VOL_LIFT)
 595             vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 596                                           row, options);
 597         else
 598             vol = volume_triangulate(PP, D, dim, matrix, point,
 599                                      row, F, options);
 600     } else {
 601         assert(nbV == (dim-row) + 1);
 602         vol = volume_simplex(PP, D, dim, matrix, point, row, options);
 603     }
 604
 605     return vol;
 606 }
 607
 608 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
 609                                 struct barvinok_options *options)
 610 {
 611     evalue ***matrix;
 612     unsigned nparam = C->Dimension;
 613     unsigned nvar = P->Dimension - C->Dimension;
 614     Param_Polyhedron *PP;
 615     unsigned PP_MaxRays = options->MaxRays;
 616     unsigned MaxRays;
 617     int i, j;
 618     Value fact;
 619     evalue *vol;
 620     int nd;
 621     struct section { Polyhedron *D; evalue *E; } *s;
 622     Param_Domain *D;
 623     Polyhedron *TC;
 624
 625     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE)
 626         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 627
 628     if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
 629         int pa = options->polynomial_approximation;
 630         assert(pa == BV_APPROX_SIGN_UPPER || pa == BV_APPROX_SIGN_LOWER);
 631
 632         P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
 633                              options->MaxRays);
 634
 635         /* Don't deflate/inflate again (on this polytope) */
 636         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 637         vol = barvinok_enumerate_with_options(P, C, options);
 638         options->polynomial_approximation = pa;
 639
 640         Polyhedron_Free(P);
 641         return vol;
 642     }
 643
 644     TC = true_context(P, C, options->MaxRays);
 645
 646     if (PP_MaxRays & POL_NO_DUAL)
 647         PP_MaxRays = 0;
 648
 649     MaxRays = options->MaxRays;
 650     POL_UNSET(options->MaxRays, POL_INTEGER);
 651
 652     value_init(fact);
 653     Factorial(nvar, &fact);
 654
 655     PP = Polyhedron2Param_Domain(P, C, PP_MaxRays);
 656
 657     for (nd = 0, D = PP->D; D; ++nd, D = D->next);
 658     s = ALLOCN(struct section, nd);
 659
 660     matrix = ALLOCN(evalue **, nvar+1);
 661     for (i = 0; i < nvar+1; ++i)
 662         matrix[i] = ALLOCN(evalue *, nvar);
 663
 664     FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, D, rVD)
 665         Polyhedron *CA, *F;
 666         struct parameter_point *point;
 667
 668         CA = align_context(D->Domain, P->Dimension, MaxRays);
 669         F = DomainIntersection(P, CA, options->MaxRays);
 670         Domain_Free(CA);
 671
 672         point = non_empty_point(D);
 673         s[i].D = rVD;
 674         s[i].E = volume_in_domain(PP, D, nvar, matrix, point, 0, F, options);
 675         Domain_Free(F);
 676         parameter_point_free(point);
 677         evalue_div(s[i].E, fact);
 678     END_FORALL_REDUCED_DOMAIN
 679     options->MaxRays = MaxRays;
 680     Polyhedron_Free(TC);
 681
 682     vol = ALLOC(evalue);
 683     value_init(vol->d);
 684     value_set_si(vol->d, 0);
 685
 686     if (nd == 0)
 687         evalue_set_si(vol, 0, 1);
 688     else {
 689         vol->x.p = new_enode(partition, 2*nd, C->Dimension);
 690         for (i = 0; i < nd; ++i) {
 691             EVALUE_SET_DOMAIN(vol->x.p->arr[2*i], s[i].D);
 692             value_clear(vol->x.p->arr[2*i+1].d);
 693             vol->x.p->arr[2*i+1] = *s[i].E;
 694             free(s[i].E);
 695         }
 696     }
 697     free(s);
 698
 699     for (i = 0; i < nvar+1; ++i)
 700         free(matrix[i]);
 701     free(matrix);
 702     Param_Polyhedron_Free(PP);
 703     value_clear(fact);
 704
 705     return vol;
 706 }