2 #include <barvinok/options.h>
3 #include <barvinok/util.h>
6 #include "lattice_width.h"
7 #include "param_util.h"
8 #include "reduce_domain.h"
10 #define ALLOC(type) (type*)malloc(sizeof(type))
11 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
13 static void clear_width_direction(struct width_direction
*wd
)
15 Vector_Free(wd
->width
);
18 Polyhedron_Free(wd
->domain
);
21 static struct width_direction_array
*new_width_direction_array(void)
23 struct width_direction_array
*dirs
= ALLOC(struct width_direction_array
);
27 dirs
->wd
= ALLOCN(struct width_direction
, dirs
->alloc
);
32 static void grow_width_direction_array(struct width_direction_array
*dirs
,
35 if (dirs
->n
+ extra
<= dirs
->alloc
)
37 dirs
->alloc
= (5*(dirs
->n
+extra
))/4;
38 dirs
->wd
= (struct width_direction
*)realloc(dirs
->wd
,
39 dirs
->alloc
* sizeof(struct width_direction
));
43 void free_width_direction_array(struct width_direction_array
*dirs
)
47 for (i
= 0; i
< dirs
->n
; ++i
)
48 clear_width_direction(&dirs
->wd
[i
]);
53 #define INT_BITS (sizeof(unsigned) * 8)
55 /* For each parametric vertex, compute cone of directions
56 * for which this vertex attains the minimal value.
58 static Matrix
**compute_vertex_dirs(Param_Polyhedron
*PP
)
61 unsigned nvar
= PP
->V
->Vertex
->NbRows
;
63 Matrix
**vertex_dirs
= ALLOCN(Matrix
*, PP
->nbV
);
65 for (i
= 0, V
= PP
->V
; V
; ++i
, V
= V
->next
) {
75 int len
= (PP
->Constraints
->NbRows
+INT_BITS
-1)/INT_BITS
;
77 n
= bit_vector_count(facets
, len
);
79 facets
= supporting_constraints(PP
->Constraints
, V
, &n
);
80 M
= Matrix_Alloc(n
, 1+nvar
+1);
81 for (k
= 0, j
= 0, kx
= 0, bx
= MSB
; j
< n
; ++k
) {
82 if (facets
[kx
] & bx
) {
83 value_set_si(M
->p
[j
][0], 1);
84 Vector_Copy(PP
->Constraints
->p
[k
]+1, M
->p
[j
++]+1, nvar
);
88 P
= Constraints2Polyhedron(M
, 0);
90 vertex_dirs
[i
] = Matrix_Alloc(P
->NbRays
-1, nvar
);
91 for (k
= 0, j
= 0; k
< P
->NbRays
; ++k
) {
92 if (value_notzero_p(P
->Ray
[k
][1+nvar
]))
94 Vector_Copy(P
->Ray
[k
]+1, vertex_dirs
[i
]->p
[j
++], nvar
);
109 static void Vector_Subtract(Value
*a
, Value a_d
,
111 Value
*c
, Value
*c_d
, int len
)
116 value_lcm(*c_d
, a_d
, b_d
);
117 value_divexact(ma
, *c_d
, a_d
);
118 value_divexact(mb
, *c_d
, b_d
);
119 value_oppose(mb
, mb
);
120 Vector_Combine(a
, b
, c
, ma
, mb
, len
);
125 /* Compute width for a given direction dir and initialize width_direction
128 static void compute_width_direction(Matrix
*V_min
, Matrix
*V_max
,
129 Value
*dir
, struct width_direction
*wd
)
131 Vector
*max
= Vector_Alloc(V_min
->NbColumns
);
132 unsigned nvar
= V_min
->NbRows
;
133 unsigned nparam
= V_min
->NbColumns
-2;
135 wd
->width
= Vector_Alloc(V_min
->NbColumns
);
136 wd
->dir
= Vector_Alloc(nvar
);
137 Vector_Copy(dir
, wd
->dir
->p
, nvar
);
143 Vector_Matrix_Product(dir
, V_max
, max
->p
);
144 Vector_Matrix_Product(dir
, V_min
, wd
->width
->p
);
145 Vector_Subtract(max
->p
, V_max
->p
[0][V_max
->NbColumns
],
146 wd
->width
->p
, V_min
->p
[0][V_min
->NbColumns
],
147 wd
->width
->p
, &wd
->width
->p
[nparam
+1],
153 Vector_Normalize(wd
->width
->p
, nparam
+2);
158 static int Vector_Compare(Value
*p1
, Value
*p2
, unsigned len
)
162 for (i
= 0; i
< len
; ++i
) {
163 int sign
= mpz_cmp(p1
[i
], p2
[i
]);
170 static int wd_width_lex_cmp(const void *va
, const void *vb
)
172 const struct width_direction
*a
= (const struct width_direction
*)va
;
173 const struct width_direction
*b
= (const struct width_direction
*)vb
;
175 return Vector_Compare(a
->width
->p
, b
->width
->p
, a
->width
->Size
);
178 static int wd_dir_lex_cmp(const void *va
, const void *vb
)
180 const struct width_direction
*a
= (const struct width_direction
*)va
;
181 const struct width_direction
*b
= (const struct width_direction
*)vb
;
183 return Vector_Compare(a
->dir
->p
, b
->dir
->p
, a
->dir
->Size
);
186 static struct width_direction_array
*
187 compute_width_directions(Param_Polyhedron
*PP
, struct barvinok_options
*options
)
189 Matrix
**vertex_dirs
;
190 Param_Vertices
*V_max
, *V_min
;
191 int i
, V_max_i
, V_min_i
;
192 unsigned nvar
= PP
->V
->Vertex
->NbRows
;
193 struct width_direction_array
*width_dirs
= new_width_direction_array();
195 vertex_dirs
= compute_vertex_dirs(PP
);
197 for (V_max
= PP
->V
; V_max
; V_max
= V_max
->next
)
198 Param_Vertex_Common_Denominator(V_max
);
200 for (V_max
= PP
->V
, V_max_i
= 0; V_max
; V_max
= V_max
->next
, V_max_i
++) {
201 for (V_min
= V_max
->next
, V_min_i
= V_max_i
+1;
203 V_min
= V_min
->next
, V_min_i
++) {
207 unsigned V_max_n
= vertex_dirs
[V_max_i
]->NbRows
;
208 unsigned V_min_n
= vertex_dirs
[V_min_i
]->NbRows
;
211 if (options
->verbose
)
212 fprintf(stderr
, "%d/%d %d/%d %d \r",
217 M
= Matrix_Alloc(V_max_n
+V_min_n
, 1+nvar
+1);
218 for (i
= 0; i
< V_max_n
; ++i
) {
219 value_set_si(M
->p
[i
][0], 1);
220 Vector_Oppose(vertex_dirs
[V_max_i
]->p
[i
], M
->p
[i
]+1, nvar
);
222 for (i
= 0; i
< V_min_n
; ++i
) {
223 value_set_si(M
->p
[V_max_n
+i
][0], 1);
224 Vector_Copy(vertex_dirs
[V_min_i
]->p
[i
], M
->p
[V_max_n
+i
]+1, nvar
);
226 C
= Constraints2Polyhedron(M
, options
->MaxRays
);
228 basis
= Cone_Integer_Hull(C
, options
);
229 grow_width_direction_array(width_dirs
, basis
->NbRows
);
230 qsort(width_dirs
->wd
, width_dirs
->n
, sizeof(struct width_direction
),
232 sorted_n
= width_dirs
->n
;
233 for (i
= 0; i
< basis
->NbRows
; ++i
) {
235 struct width_direction wd
;
240 if (bsearch(&wd
, width_dirs
->wd
, sorted_n
,
241 sizeof(struct width_direction
),
245 compute_width_direction(V_min
->Vertex
, V_max
->Vertex
,
247 &width_dirs
->wd
[width_dirs
->n
++]);
254 for (i
= 0; i
< PP
->nbV
; ++i
)
255 Matrix_Free(vertex_dirs
[i
]);
261 /* Computes the lattice width direction of a parametric polytope.
262 * The parameter space is allowed to be unbounded.
263 * Currently, the parametric polytope and the parameter space
264 * are assumed to be full-dimensional.
266 * First, we compute the parametric vertices.
267 * Then, for each pair of vertices, we construct a (rational) cone
268 * of directions for which one vertex attains the minimal value
269 * and the other vertex attians the maximal value.
270 * The candidate directions are the elements of the integer hulls
272 * The minimal direction is then obtained by computing the
273 * region in the parameter space where each direction yields
274 * a smaller (or equal) width than all the other directions.
276 * In principle, we can avoid computing candidate directions
277 * for vertices with no overlapping activity domains (possibly
278 * after opening some facets of the activity domains in the
281 * The output is a list of triples, consisting of a direction,
282 * the corresponding width and the chamber in the parameter
283 * space where this direction leads to the minimal width.
285 * The algorithm is described in "Integer points in a parameterised
286 * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
288 struct width_direction_array
*
289 Polyhedron_Lattice_Width_Directions(Polyhedron
*P
, Polyhedron
*C
,
290 struct barvinok_options
*options
)
292 Param_Polyhedron
*PP
;
293 unsigned nparam
= C
->Dimension
;
295 struct width_direction_array
*width_dirs
;
299 assert(P
->NbEq
== 0);
300 assert(C
->NbEq
== 0);
302 /* Use true context since the algorithm assumes P is non-empty
303 * for every point in the context.
305 TC
= true_context(P
, C
, options
->MaxRays
);
306 inner
= inner_point(TC
);
308 /* This is overkill, as we discard the computed chambers. */
309 PP
= Polyhedron2Param_Polyhedron(P
, TC
, options
);
311 width_dirs
= compute_width_directions(PP
, options
);
312 Param_Polyhedron_Free(PP
);
314 qsort(width_dirs
->wd
, width_dirs
->n
, sizeof(struct width_direction
),
317 for (i
= 1, j
= 1; i
< width_dirs
->n
; ++i
) {
318 /* We could also weed out width_directions that differ by a
319 * (positive) constant from another width_direction, but then
320 * we'd have to put the two width_directions on a common
323 if (Vector_Equal(width_dirs
->wd
[j
-1].width
->p
,
324 width_dirs
->wd
[i
].width
->p
, nparam
+2))
325 clear_width_direction(&width_dirs
->wd
[i
]);
327 width_dirs
->wd
[j
++] = width_dirs
->wd
[i
];
331 for (i
= 0, k
= 0; i
< width_dirs
->n
; ++i
) {
332 Matrix
*M
= Matrix_Alloc(TC
->NbConstraints
+width_dirs
->n
-(i
-k
)-1, nparam
+2);
333 for (j
= 0; j
< TC
->NbConstraints
; ++j
)
334 Vector_Copy(TC
->Constraint
[j
], M
->p
[j
], nparam
+2);
335 for (j
= 0; j
< width_dirs
->n
; ++j
) {
337 if (k
<= j
&& j
<= i
)
340 pos
= TC
->NbConstraints
+ j
;
342 pos
= TC
->NbConstraints
+ j
- (i
-k
) - 1;
343 Vector_Subtract(width_dirs
->wd
[j
].width
->p
,
344 width_dirs
->wd
[j
].width
->p
[nparam
+1],
345 width_dirs
->wd
[i
].width
->p
,
346 width_dirs
->wd
[i
].width
->p
[nparam
+1],
347 M
->p
[pos
]+1, M
->p
[pos
], nparam
+1);
348 value_set_si(M
->p
[pos
][0], 1);
349 Vector_Normalize(M
->p
[pos
]+1, nparam
+1);
350 if (!is_internal(inner
, M
->p
[pos
]))
351 value_decrement(M
->p
[pos
][nparam
+1], M
->p
[pos
][nparam
+1]);
353 width_dirs
->wd
[i
].domain
= Constraints2Polyhedron(M
, options
->MaxRays
);
354 if (emptyQ(width_dirs
->wd
[i
].domain
))
355 clear_width_direction(&width_dirs
->wd
[i
]);
357 width_dirs
->wd
[k
++] = width_dirs
->wd
[i
];
367 /* Construct evalue of chambers with their associated widths */
368 evalue
*Polyhedron_Lattice_Width(Polyhedron
*P
, Polyhedron
*C
,
369 struct barvinok_options
*options
)
372 struct evalue_section
*s
;
373 struct width_direction_array
*width_dirs
;
375 unsigned nparam
= C
->Dimension
;
377 width_dirs
= Polyhedron_Lattice_Width_Directions(P
, C
, options
);
378 s
= ALLOCN(struct evalue_section
, width_dirs
->n
);
379 for (i
= 0; i
< width_dirs
->n
; ++i
) {
380 s
[i
].D
= width_dirs
->wd
[i
].domain
;
381 width_dirs
->wd
[i
].domain
= NULL
;
382 s
[i
].E
= affine2evalue(width_dirs
->wd
[i
].width
->p
,
383 width_dirs
->wd
[i
].width
->p
[nparam
+1],
386 free_width_direction_array(width_dirs
);
388 width
= evalue_from_section_array(s
, i
);