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"
10 #define ALLOC(type) (type*)malloc(sizeof(type))
11 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
13 /* Computes an evalue representation of a coordinate
16 static evalue
*vertex2evalue(Value
*vertex
, int nparam
)
18 return affine2evalue(vertex
, vertex
[nparam
+1], nparam
);
21 static void matrix_print(evalue
***matrix
, int dim
, int *cols
,
22 const char * const *param_names
)
26 for (i
= 0; i
< dim
; ++i
)
27 for (j
= 0; j
< dim
; ++j
) {
28 int k
= cols
? cols
[j
] : j
;
29 fprintf(stderr
, "%d %d: ", i
, j
);
30 print_evalue(stderr
, matrix
[i
][k
], param_names
);
31 fprintf(stderr
, "\n");
35 /* Compute determinant using Laplace's formula.
36 * In particular, the determinant is expanded along the last row.
37 * The cols array is a list of columns that remain in the currect submatrix.
39 static evalue
*determinant_cols(evalue
***matrix
, int dim
, int *cols
)
49 evalue_copy(det
, matrix
[0][cols
[0]]);
54 evalue_set_si(&mone
, -1, 1);
56 newcols
= ALLOCN(int, dim
-1);
57 for (j
= 1; j
< dim
; ++j
)
58 newcols
[j
-1] = cols
[j
];
59 for (j
= 0; j
< dim
; ++j
) {
61 newcols
[j
-1] = cols
[j
-1];
62 tmp
= determinant_cols(matrix
, dim
-1, newcols
);
63 emul(matrix
[dim
-1][cols
[j
]], tmp
);
70 free_evalue_refs(tmp
);
75 free_evalue_refs(&mone
);
80 static evalue
*determinant(evalue
***matrix
, int dim
)
83 int *cols
= ALLOCN(int, dim
);
86 for (i
= 0; i
< dim
; ++i
)
89 det
= determinant_cols(matrix
, dim
, cols
);
96 /* Compute the facet of P that saturates constraint c.
98 static Polyhedron
*facet(Polyhedron
*P
, int c
, unsigned MaxRays
)
101 Vector
*row
= Vector_Alloc(1+P
->Dimension
+1);
102 Vector_Copy(P
->Constraint
[c
]+1, row
->p
+1, P
->Dimension
+1);
103 F
= AddConstraints(row
->p
, 1, P
, MaxRays
);
108 /* Compute a dummy Param_Domain that contains all vertices of Param_Domain D
109 * (which contains the vertices of P) that lie on the facet obtained by
110 * saturating constraint c of P
112 static Param_Domain
*face_vertices(Param_Polyhedron
*PP
, Param_Domain
*D
,
113 Polyhedron
*P
, int c
)
117 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
118 Vector
*row
= Vector_Alloc(1+nparam
+1);
119 Param_Domain
*FD
= ALLOC(Param_Domain
);
123 nv
= (PP
->nbV
- 1)/(8*sizeof(int)) + 1;
124 FD
->F
= ALLOCN(unsigned, nv
);
125 memset(FD
->F
, 0, nv
* sizeof(unsigned));
127 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
) /* _i, _ix, _bx internal counters */
129 Param_Inner_Product(P
->Constraint
[c
], V
->Vertex
, row
->p
);
130 if (First_Non_Zero(row
->p
+1, nparam
+1) == -1)
132 END_FORALL_PVertex_in_ParamPolyhedron
;
139 /* Substitute parameters by the corresponding element in subs
141 static evalue
*evalue_substitute_new(evalue
*e
, evalue
**subs
)
147 if (value_notzero_p(e
->d
)) {
153 assert(e
->x
.p
->type
== polynomial
);
156 for (i
= e
->x
.p
->size
-1; i
> 0; --i
) {
157 c
= evalue_substitute_new(&e
->x
.p
->arr
[i
], subs
);
161 emul(subs
[e
->x
.p
->pos
-1], res
);
163 c
= evalue_substitute_new(&e
->x
.p
->arr
[0], subs
);
171 struct parameter_point
{
176 struct parameter_point
*parameter_point_new(unsigned nparam
)
178 struct parameter_point
*point
= ALLOC(struct parameter_point
);
179 point
->coord
= Vector_Alloc(nparam
+1);
184 evalue
**parameter_point_evalue(struct parameter_point
*point
)
187 unsigned nparam
= point
->coord
->Size
-1;
192 point
->e
= ALLOCN(evalue
*, nparam
);
193 for (j
= 0; j
< nparam
; ++j
) {
194 point
->e
[j
] = ALLOC(evalue
);
195 value_init(point
->e
[j
]->d
);
196 evalue_set(point
->e
[j
], point
->coord
->p
[j
], point
->coord
->p
[nparam
]);
202 void parameter_point_free(struct parameter_point
*point
)
205 unsigned nparam
= point
->coord
->Size
-1;
207 Vector_Free(point
->coord
);
210 for (i
= 0; i
< nparam
; ++i
) {
211 free_evalue_refs(point
->e
[i
]);
219 /* Computes point in pameter space where polyhedron is non-empty.
221 static struct parameter_point
*non_empty_point(Param_Domain
*D
)
223 unsigned nparam
= D
->Domain
->Dimension
;
224 struct parameter_point
*point
;
227 v
= inner_point(D
->Domain
);
228 point
= parameter_point_new(nparam
);
229 Vector_Copy(v
->p
+1, point
->coord
->p
, nparam
+1);
235 static Matrix
*barycenter(Param_Polyhedron
*PP
, Param_Domain
*D
)
238 Matrix
*center
= NULL
;
248 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
251 center
= Matrix_Copy(V
->Vertex
);
252 nparam
= center
->NbColumns
- 2;
254 for (i
= 0; i
< center
->NbRows
; ++i
) {
255 value_assign(fc
, center
->p
[i
][nparam
+1]);
256 value_lcm(fc
, V
->Vertex
->p
[i
][nparam
+1],
257 ¢er
->p
[i
][nparam
+1]);
258 value_division(fc
, center
->p
[i
][nparam
+1], fc
);
259 value_division(fv
, center
->p
[i
][nparam
+1],
260 V
->Vertex
->p
[i
][nparam
+1]);
261 Vector_Combine(center
->p
[i
], V
->Vertex
->p
[i
], center
->p
[i
],
265 END_FORALL_PVertex_in_ParamPolyhedron
;
270 value_set_si(denom
, nbV
);
271 for (i
= 0; i
< center
->NbRows
; ++i
) {
272 value_multiply(center
->p
[i
][nparam
+1], center
->p
[i
][nparam
+1], denom
);
273 Vector_Normalize(center
->p
[i
], nparam
+2);
280 static Matrix
*triangulation_vertex(Param_Polyhedron
*PP
, Param_Domain
*D
,
285 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
287 END_FORALL_PVertex_in_ParamPolyhedron
;
293 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
294 * If F is a simplex, then the volume is computed of a recursive pyramid
295 * over F with the points already in matrix.
296 * Otherwise, the barycenter of F is added to matrix and the function
297 * is called recursively on the facets of F.
299 * The first row of matrix contain the _negative_ of the first point.
300 * The remaining rows of matrix contain the distance of the corresponding
301 * point to the first point.
303 static evalue
*volume_in_domain(Param_Polyhedron
*PP
, Param_Domain
*D
,
304 unsigned dim
, evalue
***matrix
,
305 struct parameter_point
*point
,
306 int row
, Polyhedron
*F
,
307 struct barvinok_options
*options
);
309 static evalue
*volume_triangulate(Param_Polyhedron
*PP
, Param_Domain
*D
,
310 unsigned dim
, evalue
***matrix
,
311 struct parameter_point
*point
,
312 int row
, Polyhedron
*F
,
313 struct barvinok_options
*options
)
320 unsigned cut_MaxRays
= options
->MaxRays
;
321 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
324 POL_UNSET(cut_MaxRays
, POL_INTEGER
);
327 evalue_set_si(&mone
, -1, 1);
329 if (options
->volume_triangulate
== BV_VOL_BARYCENTER
)
330 center
= barycenter(PP
, D
);
332 center
= triangulation_vertex(PP
, D
, F
);
333 for (j
= 0; j
< dim
; ++j
)
334 matrix
[row
][j
] = vertex2evalue(center
->p
[j
], center
->NbColumns
- 2);
335 if (options
->volume_triangulate
== BV_VOL_BARYCENTER
)
338 v
= Vector_Alloc(1+nparam
+1);
341 for (j
= 0; j
< dim
; ++j
)
342 emul(&mone
, matrix
[row
][j
]);
344 for (j
= 0; j
< dim
; ++j
)
345 eadd(matrix
[0][j
], matrix
[row
][j
]);
349 POL_ENSURE_FACETS(F
);
350 for (j
= F
->NbEq
; j
< F
->NbConstraints
; ++j
) {
353 if (First_Non_Zero(F
->Constraint
[j
]+1, dim
) == -1)
355 if (options
->volume_triangulate
!= BV_VOL_BARYCENTER
) {
356 Param_Inner_Product(F
->Constraint
[j
], center
, v
->p
);
357 if (First_Non_Zero(v
->p
+1, nparam
+1) == -1)
360 FF
= facet(F
, j
, options
->MaxRays
);
361 FD
= face_vertices(PP
, D
, F
, j
);
362 tmp
= volume_in_domain(PP
, FD
, dim
, matrix
, point
,
368 free_evalue_refs(tmp
);
372 Param_Domain_Free(FD
);
375 if (options
->volume_triangulate
!= BV_VOL_BARYCENTER
)
378 for (j
= 0; j
< dim
; ++j
) {
379 free_evalue_refs(matrix
[row
][j
]);
380 free(matrix
[row
][j
]);
383 free_evalue_refs(&mone
);
387 static evalue
*volume_simplex(Param_Polyhedron
*PP
, Param_Domain
*D
,
388 unsigned dim
, evalue
***matrix
,
389 struct parameter_point
*point
,
390 int row
, struct barvinok_options
*options
)
397 options
->stats
->volume_simplices
++;
400 evalue_set_si(&mone
, -1, 1);
403 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
) /* _ix, _bx internal counters */
404 for (j
= 0; j
< dim
; ++j
) {
405 matrix
[i
][j
] = vertex2evalue(V
->Vertex
->p
[j
],
406 V
->Vertex
->NbColumns
- 2);
408 emul(&mone
, matrix
[i
][j
]);
410 eadd(matrix
[0][j
], matrix
[i
][j
]);
413 END_FORALL_PVertex_in_ParamPolyhedron
;
415 vol
= determinant(matrix
+1, dim
);
417 val
= evalue_substitute_new(vol
, parameter_point_evalue(point
));
419 assert(value_notzero_p(val
->d
));
420 assert(value_notzero_p(val
->x
.n
));
421 if (value_neg_p(val
->x
.n
))
424 free_evalue_refs(val
);
427 for (i
= row
; i
< dim
+1; ++i
) {
428 for (j
= 0; j
< dim
; ++j
) {
429 free_evalue_refs(matrix
[i
][j
]);
434 free_evalue_refs(&mone
);
439 static evalue
*volume_triangulate_lift(Param_Polyhedron
*PP
, Param_Domain
*D
,
440 unsigned dim
, evalue
***matrix
,
441 struct parameter_point
*point
,
442 int row
, struct barvinok_options
*options
)
444 const static int MAX_TRY
=10;
449 Matrix
*FixedRays
, *M
;
457 nv
= (PP
->nbV
- 1)/(8*sizeof(int)) + 1;
458 SD
.F
= ALLOCN(unsigned, nv
);
460 FixedRays
= Matrix_Alloc(PP
->nbV
+1, 1+dim
+2);
462 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
463 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
464 Param_Vertex_Common_Denominator(V
);
465 for (i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
466 value_multiply(FixedRays
->p
[nbV
][1+i
], V
->Vertex
->p
[i
][nparam
],
467 point
->coord
->p
[nparam
]);
468 Inner_Product(V
->Vertex
->p
[i
], point
->coord
->p
, nparam
,
469 &FixedRays
->p
[nbV
][1+dim
]);
470 value_addto(FixedRays
->p
[nbV
][1+i
], FixedRays
->p
[nbV
][1+i
],
471 FixedRays
->p
[nbV
][1+dim
]);
473 value_multiply(FixedRays
->p
[nbV
][1+dim
+1], V
->Vertex
->p
[0][nparam
+1],
474 point
->coord
->p
[nparam
]);
475 value_set_si(FixedRays
->p
[nbV
][0], 1);
477 END_FORALL_PVertex_in_ParamPolyhedron
;
478 value_set_si(FixedRays
->p
[nbV
][0], 1);
479 value_set_si(FixedRays
->p
[nbV
][1+dim
], 1);
480 FixedRays
->NbRows
= nbV
+1;
485 /* Usually vol should still be NULL */
487 free_evalue_refs(vol
);
493 assert(t
<= MAX_TRY
);
496 for (i
= 0; i
< nbV
; ++i
)
497 value_set_si(FixedRays
->p
[i
][1+dim
], random_int((t
+1)*dim
*nbV
)+1);
499 M
= Matrix_Copy(FixedRays
);
500 L
= Rays2Polyhedron(M
, options
->MaxRays
);
503 POL_ENSURE_FACETS(L
);
504 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
506 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
507 if (value_negz_p(L
->Constraint
[i
][1+dim
]))
510 memset(SD
.F
, 0, nv
* sizeof(unsigned));
513 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
) /* _ix, _bx internal */
514 Inner_Product(FixedRays
->p
[nbV
]+1, L
->Constraint
[i
]+1, dim
+2, &tmp
);
515 if (value_zero_p(tmp
)) {
522 END_FORALL_PVertex_in_ParamPolyhedron
;
523 assert(r
== (dim
-row
)+1);
525 s
= volume_simplex(PP
, &SD
, dim
, matrix
, point
, row
, options
);
535 Matrix_Free(FixedRays
);
542 static evalue
*volume_in_domain(Param_Polyhedron
*PP
, Param_Domain
*D
,
543 unsigned dim
, evalue
***matrix
,
544 struct parameter_point
*point
,
545 int row
, Polyhedron
*F
,
546 struct barvinok_options
*options
)
555 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
557 END_FORALL_PVertex_in_ParamPolyhedron
;
559 if (nbV
> (dim
-row
) + 1) {
560 if (options
->volume_triangulate
== BV_VOL_LIFT
)
561 vol
= volume_triangulate_lift(PP
, D
, dim
, matrix
, point
,
564 vol
= volume_triangulate(PP
, D
, dim
, matrix
, point
,
567 assert(nbV
== (dim
-row
) + 1);
568 vol
= volume_simplex(PP
, D
, dim
, matrix
, point
, row
, options
);
574 evalue
* Param_Polyhedron_Volume(Polyhedron
*P
, Polyhedron
* C
,
575 struct barvinok_options
*options
)
578 unsigned nparam
= C
->Dimension
;
579 unsigned nvar
= P
->Dimension
- C
->Dimension
;
580 Param_Polyhedron
*PP
;
586 struct section
{ Polyhedron
*D
; evalue
*E
; } *s
;
590 if (options
->polynomial_approximation
== BV_APPROX_SIGN_NONE
)
591 options
->polynomial_approximation
= BV_APPROX_SIGN_APPROX
;
593 if (options
->polynomial_approximation
!= BV_APPROX_SIGN_APPROX
) {
594 int pa
= options
->polynomial_approximation
;
595 assert(pa
== BV_APPROX_SIGN_UPPER
|| pa
== BV_APPROX_SIGN_LOWER
);
597 P
= Polyhedron_Flate(P
, nparam
, pa
== BV_APPROX_SIGN_UPPER
,
600 /* Don't deflate/inflate again (on this polytope) */
601 options
->polynomial_approximation
= BV_APPROX_SIGN_APPROX
;
602 vol
= barvinok_enumerate_with_options(P
, C
, options
);
603 options
->polynomial_approximation
= pa
;
609 TC
= true_context(P
, C
, options
->MaxRays
);
611 MaxRays
= options
->MaxRays
;
612 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
615 Factorial(nvar
, &fact
);
617 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
619 for (nd
= 0, D
= PP
->D
; D
; ++nd
, D
= D
->next
);
620 s
= ALLOCN(struct section
, nd
);
622 matrix
= ALLOCN(evalue
**, nvar
+1);
623 for (i
= 0; i
< nvar
+1; ++i
)
624 matrix
[i
] = ALLOCN(evalue
*, nvar
);
626 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
628 struct parameter_point
*point
;
630 CA
= align_context(D
->Domain
, P
->Dimension
, MaxRays
);
631 F
= DomainIntersection(P
, CA
, options
->MaxRays
);
634 point
= non_empty_point(D
);
636 s
[i
].E
= volume_in_domain(PP
, D
, nvar
, matrix
, point
, 0, F
, options
);
638 parameter_point_free(point
);
639 evalue_div(s
[i
].E
, fact
);
640 END_FORALL_REDUCED_DOMAIN
641 options
->MaxRays
= MaxRays
;
646 value_set_si(vol
->d
, 0);
649 evalue_set_si(vol
, 0, 1);
651 vol
->x
.p
= new_enode(partition
, 2*nd
, C
->Dimension
);
652 for (i
= 0; i
< nd
; ++i
) {
653 EVALUE_SET_DOMAIN(vol
->x
.p
->arr
[2*i
], s
[i
].D
);
654 value_clear(vol
->x
.p
->arr
[2*i
+1].d
);
655 vol
->x
.p
->arr
[2*i
+1] = *s
[i
].E
;
661 for (i
= 0; i
< nvar
+1; ++i
)
664 Param_Polyhedron_Free(PP
);