iscc: add inverse operation
[barvinok.git] / lattice_width.c
blob0a193a2f1c2fbbdff18e93cf4e3ff1a2cf5fd630
1 #include <assert.h>
2 #include <stdlib.h>
3 #include <barvinok/options.h>
4 #include <barvinok/util.h>
5 #include "hilbert.h"
6 #include "hull.h"
7 #include "lattice_width.h"
8 #include "param_util.h"
9 #include "reduce_domain.h"
11 #define ALLOC(type) (type*)malloc(sizeof(type))
12 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
14 static void clear_width_direction(struct width_direction *wd)
16 Vector_Free(wd->width);
17 Vector_Free(wd->dir);
18 if (wd->domain)
19 Polyhedron_Free(wd->domain);
22 static struct width_direction_array *new_width_direction_array(void)
24 struct width_direction_array *dirs = ALLOC(struct width_direction_array);
25 assert(dirs);
26 dirs->n = 0;
27 dirs->alloc = 32;
28 dirs->wd = ALLOCN(struct width_direction, dirs->alloc);
29 assert(dirs->wd);
30 return dirs;
33 static void grow_width_direction_array(struct width_direction_array *dirs,
34 int extra)
36 if (dirs->n + extra <= dirs->alloc)
37 return;
38 dirs->alloc = (5*(dirs->n+extra))/4;
39 dirs->wd = (struct width_direction*)realloc(dirs->wd,
40 dirs->alloc * sizeof(struct width_direction));
41 assert(dirs->wd);
44 void free_width_direction_array(struct width_direction_array *dirs)
46 int i;
48 for (i = 0; i < dirs->n; ++i)
49 clear_width_direction(&dirs->wd[i]);
50 free(dirs->wd);
51 free(dirs);
54 #define INT_BITS (sizeof(unsigned) * 8)
56 /* For each parametric vertex, compute cone of directions
57 * for which this vertex attains the minimal value.
59 static Matrix **compute_vertex_dirs(Param_Polyhedron *PP)
61 int i;
62 unsigned nvar = PP->V->Vertex->NbRows;
63 Param_Vertices *V;
64 Matrix **vertex_dirs = ALLOCN(Matrix *, PP->nbV);
66 for (i = 0, V = PP->V; V; ++i, V = V->next) {
67 int kx;
68 unsigned bx;
69 int j, k;
70 unsigned *facets;
71 int n;
72 Matrix *M;
73 Polyhedron *P;
75 if (V->Facets) {
76 int len = (PP->Constraints->NbRows+INT_BITS-1)/INT_BITS;
77 facets = V->Facets;
78 n = bit_vector_count(facets, len);
79 } else
80 facets = supporting_constraints(PP->Constraints, V, &n);
81 M = Matrix_Alloc(n, 1+nvar+1);
82 for (k = 0, j = 0, kx = 0, bx = MSB; j < n; ++k) {
83 if (facets[kx] & bx) {
84 value_set_si(M->p[j][0], 1);
85 Vector_Copy(PP->Constraints->p[k]+1, M->p[j++]+1, nvar);
87 NEXT(kx, bx);
89 P = Constraints2Polyhedron(M, 0);
90 Matrix_Free(M);
91 vertex_dirs[i] = Matrix_Alloc(P->NbRays-1, nvar);
92 for (k = 0, j = 0; k < P->NbRays; ++k) {
93 if (value_notzero_p(P->Ray[k][1+nvar]))
94 continue;
95 Vector_Copy(P->Ray[k]+1, vertex_dirs[i]->p[j++], nvar);
97 Polyhedron_Free(P);
98 if (!V->Facets)
99 free(facets);
101 return vertex_dirs;
104 /* Compute
106 * c a b
107 * --- = --- - ---
108 * c_d a_d b_d
110 static void Vector_Subtract(Value *a, Value a_d,
111 Value *b, Value b_d,
112 Value *c, Value *c_d, int len)
114 Value ma, mb;
115 value_init(ma);
116 value_init(mb);
117 value_lcm(*c_d, a_d, b_d);
118 value_divexact(ma, *c_d, a_d);
119 value_divexact(mb, *c_d, b_d);
120 value_oppose(mb, mb);
121 Vector_Combine(a, b, c, ma, mb, len);
122 value_clear(ma);
123 value_clear(mb);
126 /* Compute width for a given direction dir and initialize width_direction
127 * structure.
129 static void compute_width_direction(Matrix *V_min, Matrix *V_max,
130 Value *dir, struct width_direction *wd)
132 Vector *max = Vector_Alloc(V_min->NbColumns);
133 unsigned nvar = V_min->NbRows;
134 unsigned nparam = V_min->NbColumns-2;
136 wd->width = Vector_Alloc(V_min->NbColumns);
137 wd->dir = Vector_Alloc(nvar);
138 Vector_Copy(dir, wd->dir->p, nvar);
139 wd->domain = NULL;
141 V_min->NbColumns--;
142 V_max->NbColumns--;
144 Vector_Matrix_Product(dir, V_max, max->p);
145 Vector_Matrix_Product(dir, V_min, wd->width->p);
146 Vector_Subtract(max->p, V_max->p[0][V_max->NbColumns],
147 wd->width->p, V_min->p[0][V_min->NbColumns],
148 wd->width->p, &wd->width->p[nparam+1],
149 nparam+1);
151 V_min->NbColumns++;
152 V_max->NbColumns++;
154 Vector_Normalize(wd->width->p, nparam+2);
156 Vector_Free(max);
159 static int Vector_Compare(Value *p1, Value *p2, unsigned len)
161 int i;
163 for (i = 0; i < len; ++i) {
164 int sign = mpz_cmp(p1[i], p2[i]);
165 if (sign)
166 return sign;
168 return 0;
171 static int wd_width_lex_cmp(const void *va, const void *vb)
173 const struct width_direction *a = (const struct width_direction *)va;
174 const struct width_direction *b = (const struct width_direction *)vb;
176 return Vector_Compare(a->width->p, b->width->p, a->width->Size);
179 static int wd_dir_lex_cmp(const void *va, const void *vb)
181 const struct width_direction *a = (const struct width_direction *)va;
182 const struct width_direction *b = (const struct width_direction *)vb;
184 return Vector_Compare(a->dir->p, b->dir->p, a->dir->Size);
187 static int add_vertex(Matrix *M, int n, Value *v)
189 if (n >= M->NbRows)
190 Matrix_Extend(M, 3*(M->NbRows+10)/2);
191 value_set_si(M->p[n][0], 1);
192 Vector_Copy(v, M->p[n]+1, M->NbColumns-2);
193 value_set_si(M->p[n][M->NbColumns-1], 1);
194 return n+1;
197 /* Puts the points in v that lie in P in front of the list
198 * and returns their number.
200 static int valid_vertices(Polyhedron *P, Matrix *v, int n_v)
202 int i, j, k;
203 Value tmp;
205 assert(v->NbColumns == P->Dimension+2);
206 value_init(tmp);
208 for (j = 0, k = 0; j < n_v; ++j) {
209 for (i = 0; i < P->NbConstraints; ++i) {
210 Inner_Product(v->p[j]+1, P->Constraint[i]+1, P->Dimension+1, &tmp);
211 if (value_neg_p(tmp))
212 break;
214 if (i < P->NbConstraints)
215 continue;
216 if (j != k)
217 Vector_Exchange(v->p[j]+1, v->p[k]+1, P->Dimension);
218 ++k;
221 value_clear(tmp);
222 return k;
225 static struct width_direction_array *
226 compute_width_directions(Param_Polyhedron *PP, struct barvinok_options *options)
228 Matrix **vertex_dirs;
229 Param_Vertices *V_max, *V_min;
230 int i, V_max_i, V_min_i;
231 unsigned nvar = PP->V->Vertex->NbRows;
232 struct width_direction_array *width_dirs = new_width_direction_array();
233 Matrix *all_vertices = Matrix_Alloc(nvar, 1+nvar+1);
234 int n_vertices = 0;
236 vertex_dirs = compute_vertex_dirs(PP);
238 for (V_max = PP->V; V_max; V_max = V_max->next)
239 Param_Vertex_Common_Denominator(V_max);
241 for (V_max = PP->V, V_max_i = 0; V_max; V_max = V_max->next, V_max_i++) {
242 for (V_min = V_max->next, V_min_i = V_max_i+1;
243 V_min;
244 V_min = V_min->next, V_min_i++) {
245 Matrix *M;
246 Matrix *basis;
247 Polyhedron *C;
248 unsigned V_max_n = vertex_dirs[V_max_i]->NbRows;
249 unsigned V_min_n = vertex_dirs[V_min_i]->NbRows;
250 int sorted_n;
251 int n_valid;
253 if (options->verbose)
254 fprintf(stderr, "%d/%d %d/%d %d \r",
255 V_max_i, PP->nbV,
256 V_min_i, PP->nbV,
257 width_dirs->n);
259 M = Matrix_Alloc(V_max_n+V_min_n, 1+nvar+1);
260 for (i = 0; i < V_max_n; ++i) {
261 value_set_si(M->p[i][0], 1);
262 Vector_Oppose(vertex_dirs[V_max_i]->p[i], M->p[i]+1, nvar);
264 for (i = 0; i < V_min_n; ++i) {
265 value_set_si(M->p[V_max_n+i][0], 1);
266 Vector_Copy(vertex_dirs[V_min_i]->p[i], M->p[V_max_n+i]+1, nvar);
268 C = Constraints2Polyhedron(M, options->MaxRays);
269 Matrix_Free(M);
270 n_valid = valid_vertices(C, all_vertices, n_vertices);
271 basis = Cone_Integer_Hull(C, all_vertices, n_valid, options);
272 grow_width_direction_array(width_dirs, basis->NbRows);
273 qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
274 wd_dir_lex_cmp);
275 sorted_n = width_dirs->n;
276 for (i = 0; i < basis->NbRows; ++i) {
277 Vector v;
278 struct width_direction wd;
279 Matrix *VM_min, *VM_max;
280 int pos;
282 VM_min = V_min->Vertex;
283 VM_max = V_max->Vertex;
284 pos = First_Non_Zero(basis->p[i], nvar);
285 if (value_neg_p(basis->p[i][pos])) {
286 Vector_Oppose(basis->p[i], basis->p[i], nvar);
287 VM_min = V_max->Vertex;
288 VM_max = V_min->Vertex;
291 v.Size = nvar;
292 v.p = basis->p[i];
293 wd.dir = &v;
294 if (bsearch(&wd, width_dirs->wd, sorted_n,
295 sizeof(struct width_direction),
296 wd_dir_lex_cmp))
297 continue;
299 n_vertices = add_vertex(all_vertices, n_vertices, basis->p[i]);
300 compute_width_direction(VM_min, VM_max, basis->p[i],
301 &width_dirs->wd[width_dirs->n++]);
303 Matrix_Free(basis);
304 Polyhedron_Free(C);
307 Matrix_Free(all_vertices);
309 for (i = 0; i < PP->nbV; ++i)
310 Matrix_Free(vertex_dirs[i]);
311 free(vertex_dirs);
313 return width_dirs;
316 /* Computes the lattice width direction of a parametric polytope.
317 * The parameter space is allowed to be unbounded.
318 * Currently, the parametric polytope and the parameter space
319 * are assumed to be full-dimensional.
321 * First, we compute the parametric vertices.
322 * Then, for each pair of vertices, we construct a (rational) cone
323 * of directions for which one vertex attains the minimal value
324 * and the other vertex attians the maximal value.
325 * The candidate directions are the elements of the integer hulls
326 * of these cones.
327 * The minimal direction is then obtained by computing the
328 * region in the parameter space where each direction yields
329 * a smaller (or equal) width than all the other directions.
331 * In principle, we can avoid computing candidate directions
332 * for vertices with no overlapping activity domains (possibly
333 * after opening some facets of the activity domains in the
334 * familiar way).
336 * The output is a list of triples, consisting of a direction,
337 * the corresponding width and the chamber in the parameter
338 * space where this direction leads to the minimal width.
340 * The algorithm is described in "Integer points in a parameterised
341 * polyhedron" by Friedrich Eisenbrand and Gennady Shmonin.
343 struct width_direction_array *
344 Polyhedron_Lattice_Width_Directions(Polyhedron *P, Polyhedron *C,
345 struct barvinok_options *options)
347 Param_Polyhedron *PP;
348 unsigned nparam = C->Dimension;
349 int i, j, k;
350 struct width_direction_array *width_dirs;
351 Polyhedron *TC;
352 Vector *inner;
354 assert(P->NbEq == 0);
355 assert(C->NbEq == 0);
357 /* Use true context since the algorithm assumes P is non-empty
358 * for every point in the context.
360 TC = true_context(P, C, options->MaxRays);
361 inner = inner_point(TC);
363 /* This is overkill, as we discard the computed chambers. */
364 PP = Polyhedron2Param_Polyhedron(P, TC, options);
366 width_dirs = compute_width_directions(PP, options);
367 Param_Polyhedron_Free(PP);
369 qsort(width_dirs->wd, width_dirs->n, sizeof(struct width_direction),
370 wd_width_lex_cmp);
372 for (i = 1, j = 1; i < width_dirs->n; ++i) {
373 /* We could also weed out width_directions that differ by a
374 * (positive) constant from another width_direction, but then
375 * we'd have to put the two width_directions on a common
376 * denominator first.
378 if (Vector_Equal(width_dirs->wd[j-1].width->p,
379 width_dirs->wd[i].width->p, nparam+2))
380 clear_width_direction(&width_dirs->wd[i]);
381 else
382 width_dirs->wd[j++] = width_dirs->wd[i];
384 width_dirs->n = j;
386 for (i = 0, k = 0; i < width_dirs->n; ++i) {
387 Matrix *M = Matrix_Alloc(TC->NbConstraints+width_dirs->n-(i-k)-1, nparam+2);
388 for (j = 0; j < TC->NbConstraints; ++j)
389 Vector_Copy(TC->Constraint[j], M->p[j], nparam+2);
390 for (j = 0; j < width_dirs->n; ++j) {
391 int pos;
392 if (k <= j && j <= i)
393 continue;
394 if (j < k)
395 pos = TC->NbConstraints + j;
396 else
397 pos = TC->NbConstraints + j - (i-k) - 1;
398 Vector_Subtract(width_dirs->wd[j].width->p,
399 width_dirs->wd[j].width->p[nparam+1],
400 width_dirs->wd[i].width->p,
401 width_dirs->wd[i].width->p[nparam+1],
402 M->p[pos]+1, M->p[pos], nparam+1);
403 value_set_si(M->p[pos][0], 1);
404 Vector_Normalize(M->p[pos]+1, nparam+1);
405 if (!is_internal(inner, M->p[pos]))
406 value_decrement(M->p[pos][nparam+1], M->p[pos][nparam+1]);
408 width_dirs->wd[i].domain = Constraints2Polyhedron(M, options->MaxRays);
409 if (emptyQ(width_dirs->wd[i].domain))
410 clear_width_direction(&width_dirs->wd[i]);
411 else
412 width_dirs->wd[k++] = width_dirs->wd[i];
413 Matrix_Free(M);
415 width_dirs->n = k;
416 Vector_Free(inner);
417 Polyhedron_Free(TC);
419 return width_dirs;
422 /* Construct evalue of chambers with their associated widths */
423 evalue *Polyhedron_Lattice_Width(Polyhedron *P, Polyhedron *C,
424 struct barvinok_options *options)
426 evalue *width;
427 struct evalue_section *s;
428 struct width_direction_array *width_dirs;
429 int i;
430 unsigned nparam = C->Dimension;
432 width_dirs = Polyhedron_Lattice_Width_Directions(P, C, options);
433 s = ALLOCN(struct evalue_section, width_dirs->n);
434 for (i = 0; i < width_dirs->n; ++i) {
435 s[i].D = width_dirs->wd[i].domain;
436 width_dirs->wd[i].domain = NULL;
437 s[i].E = affine2evalue(width_dirs->wd[i].width->p,
438 width_dirs->wd[i].width->p[nparam+1],
439 nparam);
441 free_width_direction_array(width_dirs);
443 width = evalue_from_section_array(s, i);
444 free(s);
446 return width;