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