2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
12 #include "isl_map_private.h"
16 * The transitive closure implementation is based on the paper
17 * "Computing the Transitive Closure of a Union of Affine Integer
18 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
22 /* Given a set of n offsets v_i (the rows of "steps"), construct a relation
23 * of the given dimension specification (Z^{n+1} -> Z^{n+1})
24 * that maps an element x to any element that can be reached
25 * by taking a non-negative number of steps along any of
26 * the extended offsets v'_i = [v_i 1].
29 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
31 * For any element in this relation, the number of steps taken
32 * is equal to the difference in the final coordinates.
34 static __isl_give isl_map
*path_along_steps(__isl_take isl_dim
*dim
,
35 __isl_keep isl_mat
*steps
)
38 struct isl_basic_map
*path
= NULL
;
46 d
= isl_dim_size(dim
, isl_dim_in
);
48 nparam
= isl_dim_size(dim
, isl_dim_param
);
50 path
= isl_basic_map_alloc_dim(isl_dim_copy(dim
), n
, d
, n
);
52 for (i
= 0; i
< n
; ++i
) {
53 k
= isl_basic_map_alloc_div(path
);
56 isl_assert(steps
->ctx
, i
== k
, goto error
);
57 isl_int_set_si(path
->div
[k
][0], 0);
60 for (i
= 0; i
< d
; ++i
) {
61 k
= isl_basic_map_alloc_equality(path
);
64 isl_seq_clr(path
->eq
[k
], 1 + isl_basic_map_total_dim(path
));
65 isl_int_set_si(path
->eq
[k
][1 + nparam
+ i
], 1);
66 isl_int_set_si(path
->eq
[k
][1 + nparam
+ d
+ i
], -1);
68 for (j
= 0; j
< n
; ++j
)
69 isl_int_set_si(path
->eq
[k
][1 + nparam
+ 2 * d
+ j
], 1);
71 for (j
= 0; j
< n
; ++j
)
72 isl_int_set(path
->eq
[k
][1 + nparam
+ 2 * d
+ j
],
76 for (i
= 0; i
< n
; ++i
) {
77 k
= isl_basic_map_alloc_inequality(path
);
80 isl_seq_clr(path
->ineq
[k
], 1 + isl_basic_map_total_dim(path
));
81 isl_int_set_si(path
->ineq
[k
][1 + nparam
+ 2 * d
+ i
], 1);
86 path
= isl_basic_map_simplify(path
);
87 path
= isl_basic_map_finalize(path
);
88 return isl_map_from_basic_map(path
);
91 isl_basic_map_free(path
);
99 /* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
100 * Return PURE_VAR if only the coefficients of the set variables are non-zero.
101 * Return IMPURE otherwise.
103 static int purity(__isl_keep isl_basic_set
*bset
, isl_int
*c
)
109 n_div
= isl_basic_set_dim(bset
, isl_dim_div
);
110 d
= isl_basic_set_dim(bset
, isl_dim_set
);
111 nparam
= isl_basic_set_dim(bset
, isl_dim_param
);
113 if (isl_seq_first_non_zero(c
+ 1 + nparam
+ d
, n_div
) != -1)
115 if (isl_seq_first_non_zero(c
+ 1, nparam
) == -1)
117 if (isl_seq_first_non_zero(c
+ 1 + nparam
, d
) == -1)
122 /* Given a set of offsets "delta", construct a relation of the
123 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
124 * is an overapproximation of the relations that
125 * maps an element x to any element that can be reached
126 * by taking a non-negative number of steps along any of
127 * the elements in "delta".
128 * That is, construct an approximation of
130 * { [x] -> [y] : exists f \in \delta, k \in Z :
131 * y = x + k [f, 1] and k >= 0 }
133 * For any element in this relation, the number of steps taken
134 * is equal to the difference in the final coordinates.
136 * In particular, let delta be defined as
138 * \delta = [p] -> { [x] : A x + a >= and B p + b >= 0 and
139 * C x + C'p + c >= 0 }
141 * then the relation is constructed as
143 * { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
144 * A f + k a >= 0 and B p + b >= 0 and k >= 1 }
145 * union { [x] -> [x] }
147 * Existentially quantified variables in \delta are currently ignored.
148 * This is safe, but leads to an additional overapproximation.
150 static __isl_give isl_map
*path_along_delta(__isl_take isl_dim
*dim
,
151 __isl_take isl_basic_set
*delta
)
153 isl_basic_map
*path
= NULL
;
162 n_div
= isl_basic_set_dim(delta
, isl_dim_div
);
163 d
= isl_basic_set_dim(delta
, isl_dim_set
);
164 nparam
= isl_basic_set_dim(delta
, isl_dim_param
);
165 path
= isl_basic_map_alloc_dim(isl_dim_copy(dim
), n_div
+ d
+ 1,
166 d
+ 1 + delta
->n_eq
, delta
->n_ineq
+ 1);
167 off
= 1 + nparam
+ 2 * (d
+ 1) + n_div
;
169 for (i
= 0; i
< n_div
+ d
+ 1; ++i
) {
170 k
= isl_basic_map_alloc_div(path
);
173 isl_int_set_si(path
->div
[k
][0], 0);
176 for (i
= 0; i
< d
+ 1; ++i
) {
177 k
= isl_basic_map_alloc_equality(path
);
180 isl_seq_clr(path
->eq
[k
], 1 + isl_basic_map_total_dim(path
));
181 isl_int_set_si(path
->eq
[k
][1 + nparam
+ i
], 1);
182 isl_int_set_si(path
->eq
[k
][1 + nparam
+ d
+ 1 + i
], -1);
183 isl_int_set_si(path
->eq
[k
][off
+ i
], 1);
186 for (i
= 0; i
< delta
->n_eq
; ++i
) {
187 int p
= purity(delta
, delta
->eq
[i
]);
190 k
= isl_basic_map_alloc_equality(path
);
193 isl_seq_clr(path
->eq
[k
], 1 + isl_basic_map_total_dim(path
));
195 isl_seq_cpy(path
->eq
[k
] + off
,
196 delta
->eq
[i
] + 1 + nparam
, d
);
197 isl_int_set(path
->eq
[k
][off
+ d
], delta
->eq
[i
][0]);
199 isl_seq_cpy(path
->eq
[k
], delta
->eq
[i
], 1 + nparam
);
202 for (i
= 0; i
< delta
->n_ineq
; ++i
) {
203 int p
= purity(delta
, delta
->ineq
[i
]);
206 k
= isl_basic_map_alloc_inequality(path
);
209 isl_seq_clr(path
->ineq
[k
], 1 + isl_basic_map_total_dim(path
));
211 isl_seq_cpy(path
->ineq
[k
] + off
,
212 delta
->ineq
[i
] + 1 + nparam
, d
);
213 isl_int_set(path
->ineq
[k
][off
+ d
], delta
->ineq
[i
][0]);
215 isl_seq_cpy(path
->ineq
[k
], delta
->ineq
[i
], 1 + nparam
);
218 k
= isl_basic_map_alloc_inequality(path
);
221 isl_seq_clr(path
->ineq
[k
], 1 + isl_basic_map_total_dim(path
));
222 isl_int_set_si(path
->ineq
[k
][0], -1);
223 isl_int_set_si(path
->ineq
[k
][off
+ d
], 1);
225 isl_basic_set_free(delta
);
226 path
= isl_basic_map_finalize(path
);
227 return isl_basic_map_union(path
,
228 isl_basic_map_identity(isl_dim_domain(dim
)));
231 isl_basic_set_free(delta
);
232 isl_basic_map_free(path
);
236 /* Given a dimenion specification Z^{n+1} -> Z^{n+1} and a parameter "param",
237 * construct a map that equates the parameter to the difference
238 * in the final coordinates and imposes that this difference is positive.
241 * { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
243 static __isl_give isl_map
*equate_parameter_to_length(__isl_take isl_dim
*dim
,
246 struct isl_basic_map
*bmap
;
251 d
= isl_dim_size(dim
, isl_dim_in
);
252 nparam
= isl_dim_size(dim
, isl_dim_param
);
253 bmap
= isl_basic_map_alloc_dim(dim
, 0, 1, 1);
254 k
= isl_basic_map_alloc_equality(bmap
);
257 isl_seq_clr(bmap
->eq
[k
], 1 + isl_basic_map_total_dim(bmap
));
258 isl_int_set_si(bmap
->eq
[k
][1 + param
], -1);
259 isl_int_set_si(bmap
->eq
[k
][1 + nparam
+ d
- 1], -1);
260 isl_int_set_si(bmap
->eq
[k
][1 + nparam
+ d
+ d
- 1], 1);
262 k
= isl_basic_map_alloc_inequality(bmap
);
265 isl_seq_clr(bmap
->ineq
[k
], 1 + isl_basic_map_total_dim(bmap
));
266 isl_int_set_si(bmap
->ineq
[k
][1 + param
], 1);
267 isl_int_set_si(bmap
->ineq
[k
][0], -1);
269 bmap
= isl_basic_map_finalize(bmap
);
270 return isl_map_from_basic_map(bmap
);
272 isl_basic_map_free(bmap
);
276 /* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
277 * construct a map that is an overapproximation of the map
278 * that takes an element from the space D to another
279 * element from the same space, such that the difference between
280 * them is a strictly positive sum of differences between images
281 * and pre-images in one of the R_i.
282 * The number of differences in the sum is equated to parameter "param".
285 * \Delta_i = { y - x | (x, y) in R_i }
287 * then the constructed map is an overapproximation of
289 * { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
290 * d = \sum_i k_i and k = \sum_i k_i > 0 }
292 * We first construct an extended mapping with an extra coordinate
293 * that indicates the number of steps taken. In particular,
294 * the difference in the last coordinate is equal to the number
295 * of steps taken to move from a domain element to the corresponding
297 * In the final step, this difference is equated to the parameter "param"
298 * and made positive. The extra coordinates are subsequently projected out.
300 * The elements of the singleton \Delta_i's are collected as the
301 * rows of the steps matrix. For all these \Delta_i's together,
302 * a single path is constructed.
303 * For each of the other \Delta_i's, we compute an overapproximation
304 * of the paths along elements of \Delta_i.
305 * Since each of these paths performs an addition, composition is
306 * symmetric and we can simply compose all resulting paths in any order.
308 static __isl_give isl_map
*construct_path(__isl_keep isl_map
*map
,
311 struct isl_mat
*steps
= NULL
;
312 struct isl_map
*path
= NULL
;
313 struct isl_map
*diff
;
314 struct isl_dim
*dim
= NULL
;
321 dim
= isl_map_get_dim(map
);
323 d
= isl_dim_size(dim
, isl_dim_in
);
324 dim
= isl_dim_add(dim
, isl_dim_in
, 1);
325 dim
= isl_dim_add(dim
, isl_dim_out
, 1);
327 path
= isl_map_identity(isl_dim_domain(isl_dim_copy(dim
)));
329 steps
= isl_mat_alloc(map
->ctx
, map
->n
, d
);
334 for (i
= 0; i
< map
->n
; ++i
) {
335 struct isl_basic_set
*delta
;
337 delta
= isl_basic_map_deltas(isl_basic_map_copy(map
->p
[i
]));
339 for (j
= 0; j
< d
; ++j
) {
342 fixed
= isl_basic_set_fast_dim_is_fixed(delta
, j
,
345 isl_basic_set_free(delta
);
354 path
= isl_map_apply_range(path
,
355 path_along_delta(isl_dim_copy(dim
), delta
));
357 isl_basic_set_free(delta
);
364 path
= isl_map_apply_range(path
,
365 path_along_steps(isl_dim_copy(dim
), steps
));
368 diff
= equate_parameter_to_length(dim
, param
);
369 path
= isl_map_intersect(path
, diff
);
370 path
= isl_map_project_out(path
, isl_dim_in
, d
, 1);
371 path
= isl_map_project_out(path
, isl_dim_out
, d
, 1);
381 /* Check whether "path" is acyclic.
382 * That is, check whether
384 * { d | d = y - x and (x,y) in path }
386 * does not contain the origin.
388 static int is_acyclic(__isl_take isl_map
*path
)
393 struct isl_set
*delta
;
395 delta
= isl_map_deltas(path
);
396 dim
= isl_set_dim(delta
, isl_dim_set
);
397 for (i
= 0; i
< dim
; ++i
)
398 delta
= isl_set_fix_si(delta
, isl_dim_set
, i
, 0);
400 acyclic
= isl_set_is_empty(delta
);
406 /* Shift variable at position "pos" up by one.
407 * That is, replace the corresponding variable v by v - 1.
409 static __isl_give isl_basic_map
*basic_map_shift_pos(
410 __isl_take isl_basic_map
*bmap
, unsigned pos
)
414 bmap
= isl_basic_map_cow(bmap
);
418 for (i
= 0; i
< bmap
->n_eq
; ++i
)
419 isl_int_sub(bmap
->eq
[i
][0], bmap
->eq
[i
][0], bmap
->eq
[i
][pos
]);
421 for (i
= 0; i
< bmap
->n_ineq
; ++i
)
422 isl_int_sub(bmap
->ineq
[i
][0],
423 bmap
->ineq
[i
][0], bmap
->ineq
[i
][pos
]);
425 for (i
= 0; i
< bmap
->n_div
; ++i
) {
426 if (isl_int_is_zero(bmap
->div
[i
][0]))
428 isl_int_sub(bmap
->div
[i
][1],
429 bmap
->div
[i
][1], bmap
->div
[i
][1 + pos
]);
435 /* Shift variable at position "pos" up by one.
436 * That is, replace the corresponding variable v by v - 1.
438 static __isl_give isl_map
*map_shift_pos(__isl_take isl_map
*map
, unsigned pos
)
442 map
= isl_map_cow(map
);
446 for (i
= 0; i
< map
->n
; ++i
) {
447 map
->p
[i
] = basic_map_shift_pos(map
->p
[i
], pos
);
451 ISL_F_CLR(map
, ISL_MAP_NORMALIZED
);
458 /* Check whether the overapproximation of the power of "map" is exactly
459 * the power of "map". Let R be "map" and A_k the overapproximation.
460 * The approximation is exact if
463 * A_k = A_{k-1} \circ R k >= 2
465 * Since A_k is known to be an overapproximation, we only need to check
468 * A_k \subset A_{k-1} \circ R k >= 2
471 static int check_power_exactness(__isl_take isl_map
*map
,
472 __isl_take isl_map
*app
, unsigned param
)
478 app_1
= isl_map_fix_si(isl_map_copy(app
), isl_dim_param
, param
, 1);
480 exact
= isl_map_is_subset(app_1
, map
);
483 if (!exact
|| exact
< 0) {
489 app_2
= isl_map_lower_bound_si(isl_map_copy(app
),
490 isl_dim_param
, param
, 2);
491 app_1
= map_shift_pos(app
, 1 + param
);
492 app_1
= isl_map_apply_range(map
, app_1
);
494 exact
= isl_map_is_subset(app_2
, app_1
);
502 /* Check whether the overapproximation of the power of "map" is exactly
503 * the power of "map", possibly after projecting out the power (if "project"
506 * If "project" is set and if "steps" can only result in acyclic paths,
509 * A = R \cup (A \circ R)
511 * where A is the overapproximation with the power projected out, i.e.,
512 * an overapproximation of the transitive closure.
513 * More specifically, since A is known to be an overapproximation, we check
515 * A \subset R \cup (A \circ R)
517 * Otherwise, we check if the power is exact.
519 static int check_exactness(__isl_take isl_map
*map
, __isl_take isl_map
*app
,
520 __isl_take isl_map
*path
, unsigned param
, int project
)
526 project
= is_acyclic(path
);
533 return check_power_exactness(map
, app
, param
);
535 map
= isl_map_project_out(map
, isl_dim_param
, param
, 1);
536 app
= isl_map_project_out(app
, isl_dim_param
, param
, 1);
538 test
= isl_map_apply_range(isl_map_copy(map
), isl_map_copy(app
));
539 test
= isl_map_union(test
, isl_map_copy(map
));
541 exact
= isl_map_is_subset(app
, test
);
555 /* Compute the positive powers of "map", or an overapproximation.
556 * The power is given by parameter "param". If the result is exact,
557 * then *exact is set to 1.
558 * If project is set, then we are actually interested in the transitive
559 * closure, so we can use a more relaxed exactness check.
561 static __isl_give isl_map
*map_power(__isl_take isl_map
*map
, unsigned param
,
562 int *exact
, int project
)
564 struct isl_set
*domain
= NULL
;
565 struct isl_set
*range
= NULL
;
566 struct isl_map
*app
= NULL
;
567 struct isl_map
*path
= NULL
;
572 map
= isl_map_remove_empty_parts(map
);
576 if (isl_map_fast_is_empty(map
))
579 isl_assert(map
->ctx
, param
< isl_map_dim(map
, isl_dim_param
), goto error
);
581 isl_map_dim(map
, isl_dim_in
) == isl_map_dim(map
, isl_dim_out
),
584 domain
= isl_map_domain(isl_map_copy(map
));
585 domain
= isl_set_coalesce(domain
);
586 range
= isl_map_range(isl_map_copy(map
));
587 range
= isl_set_coalesce(range
);
588 app
= isl_map_from_domain_and_range(isl_set_copy(domain
),
589 isl_set_copy(range
));
591 path
= construct_path(map
, param
);
592 app
= isl_map_intersect(app
, isl_map_copy(path
));
595 (*exact
= check_exactness(isl_map_copy(map
), isl_map_copy(app
),
596 isl_map_copy(path
), param
, project
)) < 0)
599 isl_set_free(domain
);
605 isl_set_free(domain
);
613 /* Compute the positive powers of "map", or an overapproximation.
614 * The power is given by parameter "param". If the result is exact,
615 * then *exact is set to 1.
617 __isl_give isl_map
*isl_map_power(__isl_take isl_map
*map
, unsigned param
,
620 return map_power(map
, param
, exact
, 0);
623 /* Compute the transitive closure of "map", or an overapproximation.
624 * If the result is exact, then *exact is set to 1.
625 * Simply compute the powers of map and then project out the parameter
626 * describing the power.
628 __isl_give isl_map
*isl_map_transitive_closure(__isl_take isl_map
*map
,
636 param
= isl_map_dim(map
, isl_dim_param
);
637 map
= isl_map_add(map
, isl_dim_param
, 1);
638 map
= map_power(map
, param
, exact
, 1);
639 map
= isl_map_project_out(map
, isl_dim_param
, param
, 1);