1 #include <assert.h>
2 #include <stdlib.h>
3 #include <isl/space.h>
4 #include <isl/set.h>
5 #include <isl/union_set.h>
6 #include <isl/polynomial.h>
7 #include <isl_set_polylib.h>
8 #include <barvinok/options.h>
9 #include <barvinok/util.h>
10 #include "hilbert.h"
11 #include "hull.h"
12 #include "lattice_width.h"
13 #include "param_util.h"
14 #include "reduce_domain.h"
15
16 #define ALLOC(type) (type*)malloc(sizeof(type))
17 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
18
clear_width_direction(struct width_direction * wd)19 static void clear_width_direction(struct width_direction *wd)
20 {
21 Vector_Free(wd->width);
22 Vector_Free(wd->dir);
23 if (wd->domain)
24 Polyhedron_Free(wd->domain);
25 }
26
new_width_direction_array(void)27 static struct width_direction_array *new_width_direction_array(void)
28 {
29 struct width_direction_array *dirs = ALLOC(struct width_direction_array);
30 assert(dirs);
31 dirs->n = 0;
32 dirs->alloc = 32;
33 dirs->wd = ALLOCN(struct width_direction, dirs->alloc);
34 assert(dirs->wd);
35 return dirs;
36 }
37
grow_width_direction_array(struct width_direction_array * dirs,int extra)38 static void grow_width_direction_array(struct width_direction_array *dirs,
39 int extra)
40 {
41 if (dirs->n + extra <= dirs->alloc)
42 return;
43 dirs->alloc = (5*(dirs->n+extra))/4;
44 dirs->wd = (struct width_direction*)realloc(dirs->wd,
45 dirs->alloc * sizeof(struct width_direction));
46 assert(dirs->wd);
47 }
48
free_width_direction_array(struct width_direction_array * dirs)49 void free_width_direction_array(struct width_direction_array *dirs)
50 {
51 int i;
52
53 for (i = 0; i < dirs->n; ++i)
54 clear_width_direction(&dirs->wd[i]);
55 free(dirs->wd);
56 free(dirs);
57 }
58
59 #define INT_BITS (sizeof(unsigned) * 8)
60
61 /* For each parametric vertex, compute cone of directions
62 * for which this vertex attains the minimal value.
63 */
compute_vertex_dirs(Param_Polyhedron * PP)64 static Matrix **compute_vertex_dirs(Param_Polyhedron *PP)
65 {
66 int i;
67 unsigned nvar = PP->V->Vertex->NbRows;
68 Param_Vertices *V;
69 Matrix **vertex_dirs = ALLOCN(Matrix *, PP->nbV);
70
71 for (i = 0, V = PP->V; V; ++i, V = V->next) {
72 int kx;
73 unsigned bx;
74 int j, k;
75 unsigned *facets;
76 int n;
77 Matrix *M;
78 Polyhedron *P;
79
80 if (V->Facets) {
81 int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
82 facets = V->Facets;
83 n = bit_vector_count(facets, len);
84 } else
85 facets = supporting_constraints(PP->Constraints, V, &n);
86 M = Matrix_Alloc(n, 1+nvar+1);
87 for (k = 0, j = 0, kx = 0, bx = MSB; j < n; ++k) {
88 if (facets[kx] & bx) {
89 value_set_si(M->p[j][0], 1);
90 Vector_Copy(PP->Constraints->p[k]+1, M->p[j++]+1, nvar);
91 }
92 NEXT(kx, bx);
93 }
94 P = Constraints2Polyhedron(M, 0);
95 Matrix_Free(M);
96 vertex_dirs[i] = Matrix_Alloc(P->NbRays-1, nvar);
97 for (k = 0, j = 0; k < P->NbRays; ++k) {
98 if (value_notzero_p(P->Ray[k][1+nvar]))
99 continue;
100 Vector_Copy(P->Ray[k]+1, vertex_dirs[i]->p[j++], nvar);
101 }
102 Polyhedron_Free(P);
103 if (!V->Facets)
104 free(facets);
105 }
106 return vertex_dirs;
107 }
108
109 /* Compute
110 *
111 * c a b
112 * --- = --- - ---
113 * c_d a_d b_d
114 */
Vector_Subtract(Value * a,Value a_d,Value * b,Value b_d,Value * c,Value * c_d,int len)115 static void Vector_Subtract(Value *a, Value a_d,
116 Value *b, Value b_d,
117 Value *c, Value *c_d, int len)
118 {
119 Value ma, mb;
120 value_init(ma);
121 value_init(mb);
122 value_lcm(*c_d, a_d, b_d);
123 value_divexact(ma, *c_d, a_d);
124 value_divexact(mb, *c_d, b_d);
125 value_oppose(mb, mb);
126 Vector_Combine(a, b, c, ma, mb, len);
127 value_clear(ma);
128 value_clear(mb);
129 }
130
131 /* Compute width for a given direction dir and initialize width_direction
132 * structure.
133 */
compute_width_direction(Matrix * V_min,Matrix * V_max,Value * dir,struct width_direction * wd)134 static void compute_width_direction(Matrix *V_min, Matrix *V_max,
135 Value *dir, struct width_direction *wd)
136 {
137 Vector *max = Vector_Alloc(V_min->NbColumns);
138 unsigned nvar = V_min->NbRows;
139 unsigned nparam = V_min->NbColumns-2;
140
141 wd->width = Vector_Alloc(V_min->NbColumns);
142 wd->dir = Vector_Alloc(nvar);
143 Vector_Copy(dir, wd->dir->p, nvar);
144 wd->domain = NULL;
145
146 V_min->NbColumns--;
147 V_max->NbColumns--;
148
149 Vector_Matrix_Product(dir, V_max, max->p);
150 Vector_Matrix_Product(dir, V_min, wd->width->p);
151 Vector_Subtract(max->p, V_max->p[0][V_max->NbColumns],
152 wd->width->p, V_min->p[0][V_min->NbColumns],
153 wd->width->p, &wd->width->p[nparam+1],
154 nparam+1);
155
156 V_min->NbColumns++;
157 V_max->NbColumns++;
158
159 Vector_Normalize(wd->width->p, nparam+2);
160
161 Vector_Free(max);
162 }
163
Vector_Compare(Value * p1,Value * p2,unsigned len)164 static int Vector_Compare(Value *p1, Value *p2, unsigned len)
165 {
166 int i;
167
168 for (i = 0; i < len; ++i) {
169 int sign = mpz_cmp(p1[i], p2[i]);
170 if (sign)
171 return sign;
172 }
173 return 0;
174 }
175
wd_width_lex_cmp(const void * va,const void * vb)176 static int wd_width_lex_cmp(const void *va, const void *vb)
177 {
178 const struct width_direction *a = (const struct width_direction *)va;
179 const struct width_direction *b = (const struct width_direction *)vb;
180
181 return Vector_Compare(a->width->p, b->width->p, a->width->Size);
182 }
183
add_vertex(Matrix * M,int n,Value * v)184 static int add_vertex(Matrix *M, int n, Value *v)
185 {
186 if (n >= M->NbRows)
187 Matrix_Extend(M, 3*(M->NbRows+10)/2);
188 value_set_si(M->p[n][0], 1);
189 Vector_Copy(v, M->p[n]+1, M->NbColumns-2);
190 value_set_si(M->p[n][M->NbColumns-1], 1);
191 return n+1;
192 }
193
194 /* Puts the points in v that lie in P in front of the list
195 * and returns their number.
196 */
valid_vertices(Polyhedron * P,Matrix * v,int n_v)197 static int valid_vertices(Polyhedron *P, Matrix *v, int n_v)
198 {
199 int i, j, k;
200 Value tmp;
201
202 assert(v->NbColumns == P->Dimension+2);
203 value_init(tmp);
204
205 for (j = 0, k = 0; j < n_v; ++j) {
206 for (i = 0; i < P->NbConstraints; ++i) {
207 Inner_Product(v->p[j]+1, P->Constraint[i]+1, P->Dimension+1, &tmp);
208 if (value_neg_p(tmp))
209 break;
210 }
211 if (i < P->NbConstraints)
212 continue;
213 if (j != k)
214 Vector_Exchange(v->p[j]+1, v->p[k]+1, P->Dimension);
215 ++k;
216 }
217
218 value_clear(tmp);
219 return k;
220 }
221
222 static struct width_direction_array *
compute_width_directions(Param_Polyhedron * PP,struct barvinok_options * options)223 compute_width_directions(Param_Polyhedron *PP, struct barvinok_options *options)
224 {
225 Matrix **vertex_dirs;
226 Param_Vertices *V_max, *V_min;
227 int i, V_max_i, V_min_i;
228 unsigned nvar = PP->V->Vertex->NbRows;
229 struct width_direction_array *width_dirs = new_width_direction_array();
230 Matrix *all_vertices = Matrix_Alloc(nvar, 1+nvar+1);
231 int n_vertices = 0;
232
233 vertex_dirs = compute_vertex_dirs(PP);
234
235 for (V_max = PP->V; V_max; V_max = V_max->next)
236 Param_Vertex_Common_Denominator(V_max);
237
238 for (V_max = PP->V, V_max_i = 0; V_max; V_max = V_max->next, V_max_i++) {
239 for (V_min = V_max->next, V_min_i = V_max_i+1;
240 V_min;
241 V_min = V_min->next, V_min_i++) {
242 Matrix *M;
243 Matrix *basis;
244 Polyhedron *C;
245 unsigned V_max_n = vertex_dirs[V_max_i]->NbRows;
246 unsigned V_min_n = vertex_dirs[V_min_i]->NbRows;
247 int n_valid;
248
249 if (options->verbose)
250 fprintf(stderr, "%d/%d %d/%d %d \r",
251 V_max_i, PP->nbV,
252 V_min_i, PP->nbV,
253 width_dirs->n);
254
255 M = Matrix_Alloc(V_max_n+V_min_n, 1+nvar+1);
256 for (i = 0; i < V_max_n; ++i) {
257 value_set_si(M->p[i][0], 1);
258 Vector_Oppose(vertex_dirs[V_max_i]->p[i], M->p[i]+1, nvar);
259 }
260 for (i = 0; i < V_min_n; ++i) {
261 value_set_si(M->p[V_max_n+i][0], 1);
262 Vector_Copy(vertex_dirs[V_min_i]->p[i], M->p[V_max_n+i]+1, nvar);
263 }
264 C = Constraints2Polyhedron(M, options->MaxRays);
265 Matrix_Free(M);
266 n_valid = valid_vertices(C, all_vertices, n_vertices);
267 basis = Cone_Integer_Hull(C, all_vertices, n_valid, options);
268 grow_width_direction_array(width_dirs, basis->NbRows);
269 for (i = 0; i < basis->NbRows; ++i) {
270 Matrix *VM_min, *VM_max;
271 int pos;
272
273 VM_min = V_min->Vertex;
274 VM_max = V_max->Vertex;
275 pos = First_Non_Zero(basis->p[i], nvar);
276 if (value_neg_p(basis->p[i][pos])) {
277 Vector_Oppose(basis->p[i], basis->p[i], nvar);
278 VM_min = V_max->Vertex;
279 VM_max = V_min->Vertex;
280 }
281
282 n_vertices = add_vertex(all_vertices, n_vertices, basis->p[i]);
283 compute_width_direction(VM_min, VM_max, basis->p[i],
284 &width_dirs->wd[width_dirs->n++]);
285 }
286 Matrix_Free(basis);
287 Polyhedron_Free(C);
288 }
289 }
290 Matrix_Free(all_vertices);
291
292 for (i = 0; i < PP->nbV; ++i)
293 Matrix_Free(vertex_dirs[i]);
294 free(vertex_dirs);
295
296 return width_dirs;
297 }
298
299 /* Computes the lattice width direction of a parametric polytope.
300 * The parameter space is allowed to be unbounded.
301 * Currently, the parametric polytope and the parameter space
302 * are assumed to be full-dimensional.
303 *
304 * First, we compute the parametric vertices.
305 * Then, for each pair of vertices, we construct a (rational) cone
306 * of directions for which one vertex attains the minimal value
307 * and the other vertex attains the maximal value.
308 * The candidate directions are the elements of the integer hulls
309 * of these cones.
310 * The minimal direction is then obtained by computing the
311 * region in the parameter space where each direction yields
312 * a smaller (or equal) width than all the other directions.
313 *
314 * In principle, we can avoid computing candidate directions
315 * for vertices with no overlapping activity domains (possibly
316 * after opening some facets of the activity domains in the
317 * familiar way).
318 *
319 * The output is a list of triples, consisting of a direction,
320 * the corresponding width and the chamber in the parameter
321 * space where this direction leads to the minimal width.
322 *
323 * The algorithm is described in "Integer points in a parameterised
324 * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
325 */
326 struct width_direction_array *
Polyhedron_Lattice_Width_Directions(Polyhedron * P,Polyhedron * C,struct barvinok_options * options)327 Polyhedron_Lattice_Width_Directions(Polyhedron *P, Polyhedron *C,
328 struct barvinok_options *options)
329 {
330 Param_Polyhedron *PP;
331 unsigned nparam = C->Dimension;
332 int i, j, k;
333 struct width_direction_array *width_dirs;
334 Polyhedron *TC;
335 Vector *inner;
336
337 assert(P->NbEq == 0);
338 assert(C->NbEq == 0);
339
340 /* Use true context since the algorithm assumes P is non-empty
341 * for every point in the context.
342 */
343 TC = true_context(P, C, options->MaxRays);
344 inner = inner_point(TC);
345
346 /* This is overkill, as we discard the computed chambers. */
347 PP = Polyhedron2Param_Polyhedron(P, TC, options);
348
349 width_dirs = compute_width_directions(PP, options);
350 Param_Polyhedron_Free(PP);
351
352 qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
353 wd_width_lex_cmp);
354
355 for (i = 1, j = 1; i < width_dirs->n; ++i) {
356 /* We could also weed out width_directions that differ by a
357 * (positive) constant from another width_direction, but then
358 * we'd have to put the two width_directions on a common
359 * denominator first.
360 */
361 if (Vector_Equal(width_dirs->wd[j-1].width->p,
362 width_dirs->wd[i].width->p, nparam+2))
363 clear_width_direction(&width_dirs->wd[i]);
364 else
365 width_dirs->wd[j++] = width_dirs->wd[i];
366 }
367 width_dirs->n = j;
368
369 for (i = 0, k = 0; i < width_dirs->n; ++i) {
370 Matrix *M = Matrix_Alloc(TC->NbConstraints+width_dirs->n-(i-k)-1, nparam+2);
371 for (j = 0; j < TC->NbConstraints; ++j)
372 Vector_Copy(TC->Constraint[j], M->p[j], nparam+2);
373 for (j = 0; j < width_dirs->n; ++j) {
374 int pos;
375 if (k <= j && j <= i)
376 continue;
377 if (j < k)
378 pos = TC->NbConstraints + j;
379 else
380 pos = TC->NbConstraints + j - (i-k) - 1;
381 Vector_Subtract(width_dirs->wd[j].width->p,
382 width_dirs->wd[j].width->p[nparam+1],
383 width_dirs->wd[i].width->p,
384 width_dirs->wd[i].width->p[nparam+1],
385 M->p[pos]+1, M->p[pos], nparam+1);
386 value_set_si(M->p[pos][0], 1);
387 Vector_Normalize(M->p[pos]+1, nparam+1);
388 if (!is_internal(inner, M->p[pos]))
389 value_decrement(M->p[pos][nparam+1], M->p[pos][nparam+1]);
390 }
391 width_dirs->wd[i].domain = Constraints2Polyhedron(M, options->MaxRays);
392 if (emptyQ(width_dirs->wd[i].domain))
393 clear_width_direction(&width_dirs->wd[i]);
394 else
395 width_dirs->wd[k++] = width_dirs->wd[i];
396 Matrix_Free(M);
397 }
398 width_dirs->n = k;
399 Vector_Free(inner);
400 Polyhedron_Free(TC);
401
402 return width_dirs;
403 }
404
405 /* Construct evalue of chambers with their associated widths */
Polyhedron_Lattice_Width(Polyhedron * P,Polyhedron * C,struct barvinok_options * options)406 evalue *Polyhedron_Lattice_Width(Polyhedron *P, Polyhedron *C,
407 struct barvinok_options *options)
408 {
409 evalue *width;
410 struct evalue_section *s;
411 struct width_direction_array *width_dirs;
412 int i;
413 unsigned nparam = C->Dimension;
414
415 width_dirs = Polyhedron_Lattice_Width_Directions(P, C, options);
416 s = ALLOCN(struct evalue_section, width_dirs->n);
417 for (i = 0; i < width_dirs->n; ++i) {
418 s[i].D = width_dirs->wd[i].domain;
419 width_dirs->wd[i].domain = NULL;
420 s[i].E = affine2evalue(width_dirs->wd[i].width->p,
421 width_dirs->wd[i].width->p[nparam+1],
422 nparam);
423 }
424 free_width_direction_array(width_dirs);
425
426 width = evalue_from_section_array(s, i);
427 free(s);
428
429 return width;
430 }
431
isl_basic_set_lattice_width(__isl_take isl_basic_set * bset)432 __isl_give isl_pw_qpolynomial *isl_basic_set_lattice_width(
433 __isl_take isl_basic_set *bset)
434 {
435 isl_ctx *ctx;
436 isl_space *dim;
437 isl_pw_qpolynomial *pwqp;
438 unsigned nparam;
439 Polyhedron *U;
440 Polyhedron *P;
441 evalue *E;
442 struct barvinok_options *options;
443 int options_allocated = 0;
444
445 if (!bset)
446 return NULL;
447
448 ctx = isl_basic_set_get_ctx(bset);
449 options = isl_ctx_peek_barvinok_options(ctx);
450 if (!options) {
451 options = barvinok_options_new_with_defaults();
452 options_allocated = 1;
453 }
454
455 nparam = isl_basic_set_dim(bset, isl_dim_param);
456 dim = isl_basic_set_get_space(bset);
457 dim = isl_space_params(dim);
458
459 U = Universe_Polyhedron(nparam);
460 P = isl_basic_set_to_polylib(bset);
461
462 E = Polyhedron_Lattice_Width(P, U, options);
463
464 pwqp = isl_pw_qpolynomial_from_evalue(dim, E);
465 isl_basic_set_free(bset);
466
467 evalue_free(E);
468 Polyhedron_Free(P);
469 Polyhedron_Free(U);
470 if (options_allocated)
471 barvinok_options_free(options);
472
473 return pwqp;
474 }
475
isl_set_lattice_width(__isl_take isl_set * set)476 __isl_give isl_pw_qpolynomial *isl_set_lattice_width(__isl_take isl_set *set)
477 {
478 if (!set)
479 return NULL;
480
481 if (isl_set_plain_is_empty(set)) {
482 isl_space *dim;
483 dim = isl_set_get_space(set);
484 dim = isl_space_domain(isl_space_from_range(dim));
485 isl_set_free(set);
486 return isl_pw_qpolynomial_zero(dim);
487 }
488
489 if (isl_set_n_basic_set(set) != 1)
490 isl_die(isl_set_get_ctx(set), isl_error_unsupported,
491 "unions not supported (yet)", goto error);
492
493 return isl_basic_set_lattice_width(isl_set_simple_hull(set));
494 error:
495 isl_set_free(set);
496 return NULL;
497 }
498
set_lw(__isl_take isl_set * set,void * user)499 static isl_stat set_lw(__isl_take isl_set *set, void *user)
500 {
501 isl_union_pw_qpolynomial **res = (isl_union_pw_qpolynomial **)user;
502 isl_pw_qpolynomial *pwqp;
503
504 pwqp = isl_set_lattice_width(set);
505 *res = isl_union_pw_qpolynomial_add_pw_qpolynomial(*res, pwqp);
506
507 return isl_stat_ok;
508 }
509
isl_union_set_lattice_width(__isl_take isl_union_set * uset)510 __isl_give isl_union_pw_qpolynomial *isl_union_set_lattice_width(
511 __isl_take isl_union_set *uset)
512 {
513 isl_space *dim;
514 isl_union_pw_qpolynomial *res;
515
516 dim = isl_union_set_get_space(uset);
517 res = isl_union_pw_qpolynomial_zero(dim);
518 if (isl_union_set_n_set(uset) > 1)
519 isl_die(isl_union_set_get_ctx(uset), isl_error_unsupported,
520 "unions not supported (yet)", goto error);
521 if (isl_union_set_foreach_set(uset, &set_lw, &res) < 0)
522 goto error;
523 isl_union_set_free(uset);
524
525 return res;
526 error:
527 isl_union_set_free(uset);
528 isl_union_pw_qpolynomial_free(res);
529 return NULL;
530 }
531