2 #include <barvinok/polylib.h>
3 #include <barvinok/barvinok.h>
4 #include <barvinok/options.h>
5 #include <barvinok/util.h>
6 #include "reduce_domain.h"
7 #include "param_util.h"
11 #define ALLOC(type) (type*)malloc(sizeof(type))
12 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
14 /* Computes an evalue representation of a coordinate
17 static evalue
*vertex2evalue(Value
*vertex
, int nparam
)
19 return affine2evalue(vertex
, vertex
[nparam
+1], nparam
);
22 static void matrix_print(evalue
***matrix
, int dim
, int *cols
,
23 const char * const *param_names
)
27 for (i
= 0; i
< dim
; ++i
)
28 for (j
= 0; j
< dim
; ++j
) {
29 int k
= cols
? cols
[j
] : j
;
30 fprintf(stderr
, "%d %d: ", i
, j
);
31 print_evalue(stderr
, matrix
[i
][k
], param_names
);
32 fprintf(stderr
, "\n");
36 /* Compute determinant using Laplace's formula.
37 * In particular, the determinant is expanded along the last row.
38 * The cols array is a list of columns that remain in the currect submatrix.
40 static evalue
*determinant_cols(evalue
***matrix
, int dim
, int *cols
)
50 evalue_copy(det
, matrix
[0][cols
[0]]);
55 evalue_set_si(&mone
, -1, 1);
57 newcols
= ALLOCN(int, dim
-1);
58 for (j
= 1; j
< dim
; ++j
)
59 newcols
[j
-1] = cols
[j
];
60 for (j
= 0; j
< dim
; ++j
) {
62 newcols
[j
-1] = cols
[j
-1];
63 tmp
= determinant_cols(matrix
, dim
-1, newcols
);
64 emul(matrix
[dim
-1][cols
[j
]], tmp
);
71 free_evalue_refs(tmp
);
76 free_evalue_refs(&mone
);
81 static evalue
*determinant(evalue
***matrix
, int dim
)
84 int *cols
= ALLOCN(int, dim
);
87 for (i
= 0; i
< dim
; ++i
)
90 det
= determinant_cols(matrix
, dim
, cols
);
97 /* Compute the facet of P that saturates constraint c.
99 static Polyhedron
*facet(Polyhedron
*P
, int c
, unsigned MaxRays
)
102 Vector
*row
= Vector_Alloc(1+P
->Dimension
+1);
103 Vector_Copy(P
->Constraint
[c
]+1, row
->p
+1, P
->Dimension
+1);
104 F
= AddConstraints(row
->p
, 1, P
, MaxRays
);
109 /* Substitute parameters by the corresponding element in subs
111 static evalue
*evalue_substitute_new(evalue
*e
, evalue
**subs
)
117 if (value_notzero_p(e
->d
)) {
123 assert(e
->x
.p
->type
== polynomial
);
126 for (i
= e
->x
.p
->size
-1; i
> 0; --i
) {
127 c
= evalue_substitute_new(&e
->x
.p
->arr
[i
], subs
);
131 emul(subs
[e
->x
.p
->pos
-1], res
);
133 c
= evalue_substitute_new(&e
->x
.p
->arr
[0], subs
);
141 struct parameter_point
{
146 struct parameter_point
*parameter_point_new(unsigned nparam
)
148 struct parameter_point
*point
= ALLOC(struct parameter_point
);
149 point
->coord
= Vector_Alloc(nparam
+1);
154 evalue
**parameter_point_evalue(struct parameter_point
*point
)
157 unsigned nparam
= point
->coord
->Size
-1;
162 point
->e
= ALLOCN(evalue
*, nparam
);
163 for (j
= 0; j
< nparam
; ++j
) {
164 point
->e
[j
] = ALLOC(evalue
);
165 value_init(point
->e
[j
]->d
);
166 evalue_set(point
->e
[j
], point
->coord
->p
[j
], point
->coord
->p
[nparam
]);
172 void parameter_point_free(struct parameter_point
*point
)
175 unsigned nparam
= point
->coord
->Size
-1;
177 Vector_Free(point
->coord
);
180 for (i
= 0; i
< nparam
; ++i
) {
181 free_evalue_refs(point
->e
[i
]);
189 /* Computes point in pameter space where polyhedron is non-empty.
191 static struct parameter_point
*non_empty_point(Param_Domain
*D
)
193 unsigned nparam
= D
->Domain
->Dimension
;
194 struct parameter_point
*point
;
197 v
= inner_point(D
->Domain
);
198 point
= parameter_point_new(nparam
);
199 Vector_Copy(v
->p
+1, point
->coord
->p
, nparam
+1);
205 static Matrix
*barycenter(Param_Polyhedron
*PP
, Param_Domain
*D
)
208 Matrix
*center
= NULL
;
218 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
221 center
= Matrix_Copy(V
->Vertex
);
222 nparam
= center
->NbColumns
- 2;
224 for (i
= 0; i
< center
->NbRows
; ++i
) {
225 value_assign(fc
, center
->p
[i
][nparam
+1]);
226 value_lcm(fc
, V
->Vertex
->p
[i
][nparam
+1],
227 ¢er
->p
[i
][nparam
+1]);
228 value_division(fc
, center
->p
[i
][nparam
+1], fc
);
229 value_division(fv
, center
->p
[i
][nparam
+1],
230 V
->Vertex
->p
[i
][nparam
+1]);
231 Vector_Combine(center
->p
[i
], V
->Vertex
->p
[i
], center
->p
[i
],
235 END_FORALL_PVertex_in_ParamPolyhedron
;
240 value_set_si(denom
, nbV
);
241 for (i
= 0; i
< center
->NbRows
; ++i
) {
242 value_multiply(center
->p
[i
][nparam
+1], center
->p
[i
][nparam
+1], denom
);
243 Vector_Normalize(center
->p
[i
], nparam
+2);
250 static Matrix
*triangulation_vertex(Param_Polyhedron
*PP
, Param_Domain
*D
,
255 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
257 END_FORALL_PVertex_in_ParamPolyhedron
;
263 /* Compute dim! times the volume of polyhedron F in Param_Domain D.
264 * If F is a simplex, then the volume is computed of a recursive pyramid
265 * over F with the points already in matrix.
266 * Otherwise, the barycenter of F is added to matrix and the function
267 * is called recursively on the facets of F.
269 * The first row of matrix contain the _negative_ of the first point.
270 * The remaining rows of matrix contain the distance of the corresponding
271 * point to the first point.
273 static evalue
*volume_in_domain(Param_Polyhedron
*PP
, Param_Domain
*D
,
274 unsigned dim
, evalue
***matrix
,
275 struct parameter_point
*point
,
276 int row
, Polyhedron
*F
,
277 struct barvinok_options
*options
);
279 static evalue
*volume_triangulate(Param_Polyhedron
*PP
, Param_Domain
*D
,
280 unsigned dim
, evalue
***matrix
,
281 struct parameter_point
*point
,
282 int row
, Polyhedron
*F
,
283 struct barvinok_options
*options
)
290 unsigned cut_MaxRays
= options
->MaxRays
;
291 unsigned nparam
= PP
->V
->Vertex
->NbColumns
-2;
294 POL_UNSET(cut_MaxRays
, POL_INTEGER
);
297 evalue_set_si(&mone
, -1, 1);
299 if (options
->volume_triangulate
== BV_VOL_BARYCENTER
)
300 center
= barycenter(PP
, D
);
302 center
= triangulation_vertex(PP
, D
, F
);
303 for (j
= 0; j
< dim
; ++j
)
304 matrix
[row
][j
] = vertex2evalue(center
->p
[j
], center
->NbColumns
- 2);
305 if (options
->volume_triangulate
== BV_VOL_BARYCENTER
)
308 v
= Vector_Alloc(1+nparam
+1);
311 for (j
= 0; j
< dim
; ++j
)
312 emul(&mone
, matrix
[row
][j
]);
314 for (j
= 0; j
< dim
; ++j
)
315 eadd(matrix
[0][j
], matrix
[row
][j
]);
319 POL_ENSURE_FACETS(F
);
320 for (j
= F
->NbEq
; j
< F
->NbConstraints
; ++j
) {
323 if (First_Non_Zero(F
->Constraint
[j
]+1, dim
) == -1)
325 if (options
->volume_triangulate
!= BV_VOL_BARYCENTER
) {
326 Param_Inner_Product(F
->Constraint
[j
], center
, v
->p
);
327 if (First_Non_Zero(v
->p
+1, nparam
+1) == -1)
330 FF
= facet(F
, j
, options
->MaxRays
);
331 FD
= Param_Polyhedron_Facet(PP
, D
, F
, j
);
332 tmp
= volume_in_domain(PP
, FD
, dim
, matrix
, point
,
338 free_evalue_refs(tmp
);
342 Param_Domain_Free(FD
);
345 if (options
->volume_triangulate
!= BV_VOL_BARYCENTER
)
348 for (j
= 0; j
< dim
; ++j
) {
349 free_evalue_refs(matrix
[row
][j
]);
350 free(matrix
[row
][j
]);
353 free_evalue_refs(&mone
);
357 static evalue
*volume_simplex(Param_Polyhedron
*PP
, Param_Domain
*D
,
358 unsigned dim
, evalue
***matrix
,
359 struct parameter_point
*point
,
360 int row
, struct barvinok_options
*options
)
367 options
->stats
->volume_simplices
++;
370 evalue_set_si(&mone
, -1, 1);
373 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
) /* _ix, _bx internal counters */
374 for (j
= 0; j
< dim
; ++j
) {
375 matrix
[i
][j
] = vertex2evalue(V
->Vertex
->p
[j
],
376 V
->Vertex
->NbColumns
- 2);
378 emul(&mone
, matrix
[i
][j
]);
380 eadd(matrix
[0][j
], matrix
[i
][j
]);
383 END_FORALL_PVertex_in_ParamPolyhedron
;
385 vol
= determinant(matrix
+1, dim
);
387 val
= evalue_substitute_new(vol
, parameter_point_evalue(point
));
389 assert(value_notzero_p(val
->d
));
390 assert(value_notzero_p(val
->x
.n
));
391 if (value_neg_p(val
->x
.n
))
394 free_evalue_refs(val
);
397 for (i
= row
; i
< dim
+1; ++i
) {
398 for (j
= 0; j
< dim
; ++j
) {
399 free_evalue_refs(matrix
[i
][j
]);
404 free_evalue_refs(&mone
);
409 static evalue
*volume_triangulate_lift(Param_Polyhedron
*PP
, Param_Domain
*D
,
410 unsigned dim
, evalue
***matrix
,
411 struct parameter_point
*point
,
412 int row
, struct barvinok_options
*options
)
414 const static int MAX_TRY
=10;
419 Matrix
*FixedRays
, *M
;
427 nv
= (PP
->nbV
- 1)/(8*sizeof(int)) + 1;
428 SD
.F
= ALLOCN(unsigned, nv
);
430 FixedRays
= Matrix_Alloc(PP
->nbV
+1, 1+dim
+2);
432 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
433 unsigned nparam
= V
->Vertex
->NbColumns
- 2;
434 Param_Vertex_Common_Denominator(V
);
435 for (i
= 0; i
< V
->Vertex
->NbRows
; ++i
) {
436 value_multiply(FixedRays
->p
[nbV
][1+i
], V
->Vertex
->p
[i
][nparam
],
437 point
->coord
->p
[nparam
]);
438 Inner_Product(V
->Vertex
->p
[i
], point
->coord
->p
, nparam
,
439 &FixedRays
->p
[nbV
][1+dim
]);
440 value_addto(FixedRays
->p
[nbV
][1+i
], FixedRays
->p
[nbV
][1+i
],
441 FixedRays
->p
[nbV
][1+dim
]);
443 value_multiply(FixedRays
->p
[nbV
][1+dim
+1], V
->Vertex
->p
[0][nparam
+1],
444 point
->coord
->p
[nparam
]);
445 value_set_si(FixedRays
->p
[nbV
][0], 1);
447 END_FORALL_PVertex_in_ParamPolyhedron
;
448 value_set_si(FixedRays
->p
[nbV
][0], 1);
449 value_set_si(FixedRays
->p
[nbV
][1+dim
], 1);
450 FixedRays
->NbRows
= nbV
+1;
455 /* Usually vol should still be NULL */
457 free_evalue_refs(vol
);
463 assert(t
<= MAX_TRY
);
466 for (i
= 0; i
< nbV
; ++i
)
467 value_set_si(FixedRays
->p
[i
][1+dim
], random_int((t
+1)*dim
*nbV
)+1);
469 M
= Matrix_Copy(FixedRays
);
470 L
= Rays2Polyhedron(M
, options
->MaxRays
);
473 POL_ENSURE_FACETS(L
);
474 for (i
= 0; i
< L
->NbConstraints
; ++i
) {
476 /* Ignore perpendicular facets, i.e., facets with 0 z-coordinate */
477 if (value_negz_p(L
->Constraint
[i
][1+dim
]))
480 memset(SD
.F
, 0, nv
* sizeof(unsigned));
483 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
) /* _ix, _bx internal */
484 Inner_Product(FixedRays
->p
[nbV
]+1, L
->Constraint
[i
]+1, dim
+2, &tmp
);
485 if (value_zero_p(tmp
)) {
492 END_FORALL_PVertex_in_ParamPolyhedron
;
493 assert(r
== (dim
-row
)+1);
495 s
= volume_simplex(PP
, &SD
, dim
, matrix
, point
, row
, options
);
505 Matrix_Free(FixedRays
);
512 static evalue
*volume_in_domain(Param_Polyhedron
*PP
, Param_Domain
*D
,
513 unsigned dim
, evalue
***matrix
,
514 struct parameter_point
*point
,
515 int row
, Polyhedron
*F
,
516 struct barvinok_options
*options
)
525 FORALL_PVertex_in_ParamPolyhedron(V
, D
, PP
)
527 END_FORALL_PVertex_in_ParamPolyhedron
;
529 if (nbV
> (dim
-row
) + 1) {
530 if (options
->volume_triangulate
== BV_VOL_LIFT
)
531 vol
= volume_triangulate_lift(PP
, D
, dim
, matrix
, point
,
534 vol
= volume_triangulate(PP
, D
, dim
, matrix
, point
,
537 assert(nbV
== (dim
-row
) + 1);
538 vol
= volume_simplex(PP
, D
, dim
, matrix
, point
, row
, options
);
544 evalue
* Param_Polyhedron_Volume(Polyhedron
*P
, Polyhedron
* C
,
545 struct barvinok_options
*options
)
548 unsigned nparam
= C
->Dimension
;
549 unsigned nvar
= P
->Dimension
- C
->Dimension
;
550 Param_Polyhedron
*PP
;
556 struct evalue_section
*s
;
560 if (options
->polynomial_approximation
== BV_APPROX_SIGN_NONE
)
561 options
->polynomial_approximation
= BV_APPROX_SIGN_APPROX
;
563 if (options
->polynomial_approximation
!= BV_APPROX_SIGN_APPROX
) {
564 int pa
= options
->polynomial_approximation
;
565 assert(pa
== BV_APPROX_SIGN_UPPER
|| pa
== BV_APPROX_SIGN_LOWER
);
567 P
= Polyhedron_Flate(P
, nparam
, pa
== BV_APPROX_SIGN_UPPER
,
570 /* Don't deflate/inflate again (on this polytope) */
571 options
->polynomial_approximation
= BV_APPROX_SIGN_APPROX
;
572 vol
= barvinok_enumerate_with_options(P
, C
, options
);
573 options
->polynomial_approximation
= pa
;
579 TC
= true_context(P
, C
, options
->MaxRays
);
581 MaxRays
= options
->MaxRays
;
582 POL_UNSET(options
->MaxRays
, POL_INTEGER
);
585 Factorial(nvar
, &fact
);
587 PP
= Polyhedron2Param_Polyhedron(P
, C
, options
);
589 for (nd
= 0, D
= PP
->D
; D
; ++nd
, D
= D
->next
);
590 s
= ALLOCN(struct evalue_section
, nd
);
592 matrix
= ALLOCN(evalue
**, nvar
+1);
593 for (i
= 0; i
< nvar
+1; ++i
)
594 matrix
[i
] = ALLOCN(evalue
*, nvar
);
596 FORALL_REDUCED_DOMAIN(PP
, TC
, nd
, options
, i
, D
, rVD
)
598 struct parameter_point
*point
;
600 CA
= align_context(D
->Domain
, P
->Dimension
, MaxRays
);
601 F
= DomainIntersection(P
, CA
, options
->MaxRays
);
604 point
= non_empty_point(D
);
606 s
[i
].E
= volume_in_domain(PP
, D
, nvar
, matrix
, point
, 0, F
, options
);
608 parameter_point_free(point
);
609 evalue_div(s
[i
].E
, fact
);
610 END_FORALL_REDUCED_DOMAIN
611 options
->MaxRays
= MaxRays
;
614 vol
= evalue_from_section_array(s
, nd
);
617 for (i
= 0; i
< nvar
+1; ++i
)
620 Param_Polyhedron_Free(PP
);