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,
 402                                 struct barvinok_options *options);
 403
 404 static evalue *volume_triangulate(Param_Polyhedron *PP, Param_Domain *D,
 405                                   unsigned dim, evalue ***matrix,
 406                                   struct parameter_point *point, Polyhedron *C,
 407                                   int row, Polyhedron *F,
 408                                   struct barvinok_options *options)
 409 {
 410     int j;
 411     evalue *tmp;
 412     evalue *vol;
 413     evalue mone;
 414     Matrix *center;
 415     unsigned cut_MaxRays = options->MaxRays;
 416     unsigned nparam = C->Dimension;
 417     Matrix *M = NULL;
 418
 419     POL_UNSET(cut_MaxRays, POL_INTEGER);
 420
 421     value_init(mone.d);
 422     evalue_set_si(&mone, -1, 1);
 423
 424     center = barycenter(PP, D);
 425     for (j = 0; j < dim; ++j)
 426         matrix[row][j] = vertex2evalue(center->p[j], center->NbColumns - 2);
 427
 428     if (row == 0) {
 429         for (j = 0; j < dim; ++j)
 430             emul(&mone, matrix[row][j]);
 431     } else {
 432         for (j = 0; j < dim; ++j)
 433             eadd(matrix[0][j], matrix[row][j]);
 434     }
 435
 436     if (!point)
 437         M = Matrix_Alloc(1, nparam+2);
 438
 439     vol = NULL;
 440     POL_ENSURE_FACETS(F);
 441     for (j = F->NbEq; j < F->NbConstraints; ++j) {
 442         Polyhedron *FC;
 443         Polyhedron *FF;
 444         Param_Domain *FD;
 445         if (First_Non_Zero(F->Constraint[j]+1, dim) == -1)
 446             continue;
 447         if (point)
 448             FC = C;
 449         else {
 450             Polyhedron *cut;
 451             int pos;
 452             Param_Inner_Product(F->Constraint[j], center, M->p[0]);
 453             pos = First_Non_Zero(M->p[0]+1, nparam+1);
 454             if (pos == -1)
 455                 /* barycenter always lies on facet */
 456                 continue;
 457             if (pos == nparam)
 458                 FC = C;
 459             else {
 460                 value_set_si(M->p[0][0], 0);
 461                 cut = Constraints2Polyhedron(M, cut_MaxRays);
 462                 FC = DomainDifference(C, cut, options->MaxRays);
 463                 Polyhedron_Free(cut);
 464                 POL_ENSURE_VERTICES(FC);
 465                 if (emptyQ(FC)) {
 466                     /* barycenter lies on facet for all parameters in C */
 467                     Polyhedron_Free(FC);
 468                     continue;
 469                 }
 470             }
 471         }
 472         FF = facet(F, j, options->MaxRays);
 473         FD = face_vertices(PP, D, F, j);
 474         tmp = volume_in_domain(PP, FD, dim, matrix, point, FC,
 475                                row+1, FF, options);
 476         if (FC != C)
 477             Domain_Free(FC);
 478         if (!vol)
 479             vol = tmp;
 480         else {
 481             eadd(tmp, vol);
 482             free_evalue_refs(tmp);
 483             free(tmp);
 484         }
 485         Polyhedron_Free(FF);
 486         Param_Domain_Free(FD);
 487     }
 488
 489     Matrix_Free(center);
 490     if (!point)
 491         Matrix_Free(M);
 492
 493     for (j = 0; j < dim; ++j) {
 494         free_evalue_refs(matrix[row][j]);
 495         free(matrix[row][j]);
 496     }
 497
 498     free_evalue_refs(&mone);
 499     return vol;
 500 }
 501
 502 static evalue *volume_simplex(Param_Polyhedron *PP, Param_Domain *D,
 503                                 unsigned dim, evalue ***matrix,
 504                                 struct parameter_point *point,
 505                                 int row, unsigned MaxRays)
 506 {
 507     evalue mone;
 508     Param_Vertices *V;
 509     evalue *vol, *val;
 510     int i, j;
 511
 512     if (!point)
 513         return evalue_zero();
 514
 515     value_init(mone.d);
 516     evalue_set_si(&mone, -1, 1);
 517
 518     i = row;
 519     FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal counters */
 520         for (j = 0; j < dim; ++j) {
 521             matrix[i][j] = vertex2evalue(V->Vertex->p[j],
 522                                            V->Vertex->NbColumns - 2);
 523             if (i == 0)
 524                 emul(&mone, matrix[i][j]);
 525             else
 526                 eadd(matrix[0][j], matrix[i][j]);
 527         }
 528         ++i;
 529     END_FORALL_PVertex_in_ParamPolyhedron;
 530
 531     vol = determinant(matrix+1, dim);
 532
 533     val = evalue_substitute(vol, parameter_point_evalue(point));
 534
 535     assert(value_notzero_p(val->d));
 536     assert(value_notzero_p(val->x.n));
 537     if (value_neg_p(val->x.n))
 538         emul(&mone, vol);
 539
 540     free_evalue_refs(val);
 541     free(val);
 542
 543     for (i = row; i < dim+1; ++i) {
 544         for (j = 0; j < dim; ++j) {
 545             free_evalue_refs(matrix[i][j]);
 546             free(matrix[i][j]);
 547         }
 548     }
 549
 550     free_evalue_refs(&mone);
 551
 552     return vol;
 553 }
 554
 555 static evalue *volume_triangulate_lift(Param_Polyhedron *PP, Param_Domain *D,
 556                                         unsigned dim, evalue ***matrix,
 557                                         struct parameter_point *point,
 558                                         int row, unsigned MaxRays)
 559 {
 560     const static int MAX_TRY=10;
 561     Param_Vertices *V;
 562     int nbV, nv;
 563     int i;
 564     int t = 0;
 565     Matrix *FixedRays, *M;
 566     Polyhedron *L;
 567     Param_Domain SD;
 568     Value tmp;
 569     evalue *s, *vol;
 570
 571     SD.Domain = 0;
 572     SD.next = 0;
 573     nv = (PP->nbV - 1)/(8*sizeof(int)) + 1;
 574     SD.F = ALLOCN(unsigned, nv);
 575
 576     FixedRays = Matrix_Alloc(PP->nbV+1, 1+dim+2);
 577     nbV = 0;
 578     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 579         unsigned nparam = V->Vertex->NbColumns - 2;
 580         Param_Vertex_Common_Denominator(V);
 581         for (i = 0; i < V->Vertex->NbRows; ++i) {
 582             value_multiply(FixedRays->p[nbV][1+i], V->Vertex->p[i][nparam],
 583                            point->coord->p[nparam]);
 584             Inner_Product(V->Vertex->p[i], point->coord->p, nparam,
 585                           &FixedRays->p[nbV][1+dim]);
 586             value_addto(FixedRays->p[nbV][1+i], FixedRays->p[nbV][1+i],
 587                         FixedRays->p[nbV][1+dim]);
 588         }
 589         value_multiply(FixedRays->p[nbV][1+dim+1], V->Vertex->p[0][nparam+1],
 590                        point->coord->p[nparam]);
 591         value_set_si(FixedRays->p[nbV][0], 1);
 592         ++nbV;
 593     END_FORALL_PVertex_in_ParamPolyhedron;
 594     value_set_si(FixedRays->p[nbV][0], 1);
 595     value_set_si(FixedRays->p[nbV][1+dim], 1);
 596     FixedRays->NbRows = nbV+1;
 597
 598     value_init(tmp);
 599     if (0) {
 600 try_again:
 601         /* Usually vol should still be NULL */
 602         if (vol) {
 603             free_evalue_refs(vol);
 604             free(vol);
 605         }
 606         Polyhedron_Free(L);
 607         ++t;
 608     }
 609     assert(t <= MAX_TRY);
 610     vol = NULL;
 611
 612     for (i = 0; i < nbV; ++i)
 613         value_set_si(FixedRays->p[i][1+dim], random_int((t+1)*dim*nbV)+1);
 614
 615     M = Matrix_Copy(FixedRays);
 616     L = Rays2Polyhedron(M, MaxRays);
 617     Matrix_Free(M);
 618
 619     POL_ENSURE_FACETS(L);
 620     for (i = 0; i < L->NbConstraints; ++i) {
 621         int r;
 622         /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
 623         if (value_negz_p(L->Constraint[i][1+dim]))
 624             continue;
 625
 626         memset(SD.F, 0, nv * sizeof(unsigned));
 627         nbV = 0;
 628         r = 0;
 629         FORALL_PVertex_in_ParamPolyhedron(V, D, PP) /* _ix, _bx internal */
 630             Inner_Product(FixedRays->p[nbV]+1, L->Constraint[i]+1, dim+2, &tmp);
 631             if (value_zero_p(tmp)) {
 632                 if (r > dim-row)
 633                     goto try_again;
 634                 SD.F[_ix] |= _bx;
 635                 ++r;
 636             }
 637             ++nbV;
 638         END_FORALL_PVertex_in_ParamPolyhedron;
 639         assert(r == (dim-row)+1);
 640
 641         s = volume_simplex(PP, &SD, dim, matrix, point, row, MaxRays);
 642         if (!vol)
 643             vol = s;
 644         else {
 645             eadd(s, vol);
 646             free_evalue_refs(s);
 647             free(s);
 648         }
 649     }
 650     Polyhedron_Free(L);
 651     Matrix_Free(FixedRays);
 652     value_clear(tmp);
 653     free(SD.F);
 654
 655     return vol;
 656 }
 657
 658 static evalue *volume_in_domain(Param_Polyhedron *PP, Param_Domain *D,
 659                                 unsigned dim, evalue ***matrix,
 660                                 struct parameter_point *point, Polyhedron *C,
 661                                 int row, Polyhedron *F,
 662                                 struct barvinok_options *options)
 663 {
 664     int nbV;
 665     Param_Vertices *V;
 666     evalue *vol;
 667     int point_computed = 0;
 668
 669     if (!point) {
 670         point = non_empty_point(PP, D, F, C, options->MaxRays);
 671         if (point)
 672             point_computed = 1;
 673     }
 674
 675     nbV = 0;
 676     FORALL_PVertex_in_ParamPolyhedron(V, D, PP)
 677         ++nbV;
 678     END_FORALL_PVertex_in_ParamPolyhedron;
 679
 680     if (nbV > (dim-row) + 1) {
 681         if (point && options->volume_triangulate_lift)
 682             vol = volume_triangulate_lift(PP, D, dim, matrix, point,
 683                                           row, options->MaxRays);
 684         else
 685             vol = volume_triangulate(PP, D, dim, matrix, point, C,
 686                                      row, F, options);
 687     } else {
 688         assert(nbV == (dim-row) + 1);
 689         vol = volume_simplex(PP, D, dim, matrix, point, row, options->MaxRays);
 690     }
 691
 692     if (point_computed)
 693         parameter_point_free(point);
 694
 695     return vol;
 696 }
 697
 698 evalue* Param_Polyhedron_Volume(Polyhedron *P, Polyhedron* C,
 699                                 struct barvinok_options *options)
 700 {
 701     evalue ***matrix;
 702     unsigned nparam = C->Dimension;
 703     unsigned nvar = P->Dimension - C->Dimension;
 704     Param_Polyhedron *PP;
 705     unsigned PP_MaxRays = options->MaxRays;
 706     unsigned MaxRays;
 707     int i, j;
 708     Value fact;
 709     evalue *vol;
 710     int nd;
 711     struct section { Polyhedron *D; evalue *E; } *s;
 712     Polyhedron **fVD;
 713     Param_Domain *D, *next;
 714     Polyhedron *CA, *F;
 715
 716     if (options->polynomial_approximation == BV_APPROX_SIGN_NONE)
 717         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 718
 719     if (options->polynomial_approximation != BV_APPROX_SIGN_APPROX) {
 720         int pa = options->polynomial_approximation;
 721         assert(pa == BV_APPROX_SIGN_UPPER || pa == BV_APPROX_SIGN_LOWER);
 722
 723         P = Polyhedron_Flate(P, nparam, pa == BV_APPROX_SIGN_UPPER,
 724                              options->MaxRays);
 725
 726         /* Don't deflate/inflate again (on this polytope) */
 727         options->polynomial_approximation = BV_APPROX_SIGN_APPROX;
 728         vol = barvinok_enumerate_with_options(P, C, options);
 729         options->polynomial_approximation = pa;
 730
 731         Polyhedron_Free(P);
 732         return vol;
 733     }
 734
 735     if (PP_MaxRays & POL_NO_DUAL)
 736         PP_MaxRays = 0;
 737
 738     MaxRays = options->MaxRays;
 739     POL_UNSET(options->MaxRays, POL_INTEGER);
 740
 741     value_init(fact);
 742     Factorial(nvar, &fact);
 743
 744     PP = Polyhedron2Param_Domain(P, C, PP_MaxRays);
 745
 746     for (nd = 0, D = PP->D; D; ++nd, D = D->next);
 747     s = ALLOCN(struct section, nd);
 748     fVD = ALLOCN(Polyhedron *, nd);
 749
 750     matrix = ALLOCN(evalue **, nvar+1);
 751     for (i = 0; i < nvar+1; ++i)
 752         matrix[i] = ALLOCN(evalue *, nvar);
 753
 754     for (nd = 0, D = PP->D; D; D = next) {
 755         Polyhedron *rVD = reduce_domain(D->Domain, NULL, NULL, fVD, nd, options);
 756
 757         next = D->next;
 758
 759         if (!rVD)
 760             continue;
 761
 762         CA = align_context(D->Domain, P->Dimension, MaxRays);
 763         F = DomainIntersection(P, CA, options->MaxRays);
 764         Domain_Free(CA);
 765
 766         s[nd].D = rVD;
 767         s[nd].E = volume_in_domain(PP, D, nvar, matrix, NULL, rVD,
 768                                    0, F, options);
 769         Domain_Free(F);
 770         evalue_div(s[nd].E, fact);
 771
 772         ++nd;
 773     }
 774     options->MaxRays = MaxRays;
 775
 776     vol = ALLOC(evalue);
 777     value_init(vol->d);
 778     value_set_si(vol->d, 0);
 779
 780     if (nd == 0)
 781         evalue_set_si(vol, 0, 1);
 782     else {
 783         vol->x.p = new_enode(partition, 2*nd, C->Dimension);
 784         for (i = 0; i < nd; ++i) {
 785             EVALUE_SET_DOMAIN(vol->x.p->arr[2*i], s[i].D);
 786             value_clear(vol->x.p->arr[2*i+1].d);
 787             vol->x.p->arr[2*i+1] = *s[i].E;
 788             free(s[i].E);
 789             Domain_Free(fVD[i]);
 790         }
 791     }
 792     free(s);
 793     free(fVD);
 794
 795     for (i = 0; i < nvar+1; ++i)
 796         free(matrix[i]);
 797     free(matrix);
 798     Param_Polyhedron_Free(PP);
 799     value_clear(fact);
 800
 801     return vol;
 802 }