volume.c

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