brintos / llvm-project-archived public Read only

1/*2 * Copyright 2010      INRIA Saclay3 *4 * Use of this software is governed by the MIT license5 *6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,8 * 91893 Orsay, France 9 */10 11#include <isl_ctx_private.h>12#include <isl_map_private.h>13#include <isl/map.h>14#include <isl_seq.h>15#include <isl_space_private.h>16#include <isl_lp_private.h>17#include <isl/union_map.h>18#include <isl_mat_private.h>19#include <isl_vec_private.h>20#include <isl_options_private.h>21#include <isl_tarjan.h>22 23isl_bool isl_map_is_transitively_closed(__isl_keep isl_map *map)24{25	isl_map *map2;26	isl_bool closed;27 28	map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));29	closed = isl_map_is_subset(map2, map);30	isl_map_free(map2);31 32	return closed;33}34 35isl_bool isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)36{37	isl_union_map *umap2;38	isl_bool closed;39 40	umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),41					  isl_union_map_copy(umap));42	closed = isl_union_map_is_subset(umap2, umap);43	isl_union_map_free(umap2);44 45	return closed;46}47 48/* Given a map that represents a path with the length of the path49 * encoded as the difference between the last output coordindate50 * and the last input coordinate, set this length to either51 * exactly "length" (if "exactly" is set) or at least "length"52 * (if "exactly" is not set).53 */54static __isl_give isl_map *set_path_length(__isl_take isl_map *map,55	int exactly, int length)56{57	isl_space *space;58	struct isl_basic_map *bmap;59	isl_size d;60	isl_size nparam;61	isl_size total;62	int k;63	isl_int *c;64 65	if (!map)66		return NULL;67 68	space = isl_map_get_space(map);69	d = isl_space_dim(space, isl_dim_in);70	nparam = isl_space_dim(space, isl_dim_param);71	total = isl_space_dim(space, isl_dim_all);72	if (d < 0 || nparam < 0 || total < 0)73		space = isl_space_free(space);74	bmap = isl_basic_map_alloc_space(space, 0, 1, 1);75	if (exactly) {76		k = isl_basic_map_alloc_equality(bmap);77		if (k < 0)78			goto error;79		c = bmap->eq[k];80	} else {81		k = isl_basic_map_alloc_inequality(bmap);82		if (k < 0)83			goto error;84		c = bmap->ineq[k];85	}86	isl_seq_clr(c, 1 + total);87	isl_int_set_si(c[0], -length);88	isl_int_set_si(c[1 + nparam + d - 1], -1);89	isl_int_set_si(c[1 + nparam + d + d - 1], 1);90 91	bmap = isl_basic_map_finalize(bmap);92	map = isl_map_intersect(map, isl_map_from_basic_map(bmap));93 94	return map;95error:96	isl_basic_map_free(bmap);97	isl_map_free(map);98	return NULL;99}100 101/* Check whether the overapproximation of the power of "map" is exactly102 * the power of "map".  Let R be "map" and A_k the overapproximation.103 * The approximation is exact if104 *105 *	A_1 = R106 *	A_k = A_{k-1} \circ R			k >= 2107 *108 * Since A_k is known to be an overapproximation, we only need to check109 *110 *	A_1 \subset R111 *	A_k \subset A_{k-1} \circ R		k >= 2112 *113 * In practice, "app" has an extra input and output coordinate114 * to encode the length of the path.  So, we first need to add115 * this coordinate to "map" and set the length of the path to116 * one.117 */118static isl_bool check_power_exactness(__isl_take isl_map *map,119	__isl_take isl_map *app)120{121	isl_bool exact;122	isl_map *app_1;123	isl_map *app_2;124 125	map = isl_map_add_dims(map, isl_dim_in, 1);126	map = isl_map_add_dims(map, isl_dim_out, 1);127	map = set_path_length(map, 1, 1);128 129	app_1 = set_path_length(isl_map_copy(app), 1, 1);130 131	exact = isl_map_is_subset(app_1, map);132	isl_map_free(app_1);133 134	if (!exact || exact < 0) {135		isl_map_free(app);136		isl_map_free(map);137		return exact;138	}139 140	app_1 = set_path_length(isl_map_copy(app), 0, 1);141	app_2 = set_path_length(app, 0, 2);142	app_1 = isl_map_apply_range(map, app_1);143 144	exact = isl_map_is_subset(app_2, app_1);145 146	isl_map_free(app_1);147	isl_map_free(app_2);148 149	return exact;150}151 152/* Check whether the overapproximation of the power of "map" is exactly153 * the power of "map", possibly after projecting out the power (if "project"154 * is set).155 *156 * If "project" is set and if "steps" can only result in acyclic paths,157 * then we check158 *159 *	A = R \cup (A \circ R)160 *161 * where A is the overapproximation with the power projected out, i.e.,162 * an overapproximation of the transitive closure.163 * More specifically, since A is known to be an overapproximation, we check164 *165 *	A \subset R \cup (A \circ R)166 *167 * Otherwise, we check if the power is exact.168 *169 * Note that "app" has an extra input and output coordinate to encode170 * the length of the part.  If we are only interested in the transitive171 * closure, then we can simply project out these coordinates first.172 */173static isl_bool check_exactness(__isl_take isl_map *map,174	__isl_take isl_map *app, int project)175{176	isl_map *test;177	isl_bool exact;178	isl_size d;179 180	if (!project)181		return check_power_exactness(map, app);182 183	d = isl_map_dim(map, isl_dim_in);184	if (d < 0)185		app = isl_map_free(app);186	app = set_path_length(app, 0, 1);187	app = isl_map_project_out(app, isl_dim_in, d, 1);188	app = isl_map_project_out(app, isl_dim_out, d, 1);189 190	app = isl_map_reset_space(app, isl_map_get_space(map));191 192	test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));193	test = isl_map_union(test, isl_map_copy(map));194 195	exact = isl_map_is_subset(app, test);196 197	isl_map_free(app);198	isl_map_free(test);199 200	isl_map_free(map);201 202	return exact;203}204 205/*206 * The transitive closure implementation is based on the paper207 * "Computing the Transitive Closure of a Union of Affine Integer208 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and209 * Albert Cohen.210 */211 212/* Given a set of n offsets v_i (the rows of "steps"), construct a relation213 * of the given dimension specification (Z^{n+1} -> Z^{n+1})214 * that maps an element x to any element that can be reached215 * by taking a non-negative number of steps along any of216 * the extended offsets v'_i = [v_i 1].217 * That is, construct218 *219 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }220 *221 * For any element in this relation, the number of steps taken222 * is equal to the difference in the final coordinates.223 */224static __isl_give isl_map *path_along_steps(__isl_take isl_space *space,225	__isl_keep isl_mat *steps)226{227	int i, j, k;228	struct isl_basic_map *path = NULL;229	isl_size d;230	unsigned n;231	isl_size nparam;232	isl_size total;233 234	d = isl_space_dim(space, isl_dim_in);235	nparam = isl_space_dim(space, isl_dim_param);236	if (d < 0 || nparam < 0 || !steps)237		goto error;238 239	n = steps->n_row;240 241	path = isl_basic_map_alloc_space(isl_space_copy(space), n, d, n);242 243	for (i = 0; i < n; ++i) {244		k = isl_basic_map_alloc_div(path);245		if (k < 0)246			goto error;247		isl_assert(steps->ctx, i == k, goto error);248		isl_int_set_si(path->div[k][0], 0);249	}250 251	total = isl_basic_map_dim(path, isl_dim_all);252	if (total < 0)253		goto error;254	for (i = 0; i < d; ++i) {255		k = isl_basic_map_alloc_equality(path);256		if (k < 0)257			goto error;258		isl_seq_clr(path->eq[k], 1 + total);259		isl_int_set_si(path->eq[k][1 + nparam + i], 1);260		isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);261		if (i == d - 1)262			for (j = 0; j < n; ++j)263				isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);264		else265			for (j = 0; j < n; ++j)266				isl_int_set(path->eq[k][1 + nparam + 2 * d + j],267					    steps->row[j][i]);268	}269 270	for (i = 0; i < n; ++i) {271		k = isl_basic_map_alloc_inequality(path);272		if (k < 0)273			goto error;274		isl_seq_clr(path->ineq[k], 1 + total);275		isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);276	}277 278	isl_space_free(space);279 280	path = isl_basic_map_simplify(path);281	path = isl_basic_map_finalize(path);282	return isl_map_from_basic_map(path);283error:284	isl_space_free(space);285	isl_basic_map_free(path);286	return NULL;287}288 289#define IMPURE		0290#define PURE_PARAM	1291#define PURE_VAR	2292#define MIXED		3293 294/* Check whether the parametric constant term of constraint c is never295 * positive in "bset".296 */297static isl_bool parametric_constant_never_positive(298	__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity)299{300	isl_size d;301	isl_size n_div;302	isl_size nparam;303	isl_size total;304	int i;305	int k;306	isl_bool empty;307 308	n_div = isl_basic_set_dim(bset, isl_dim_div);309	d = isl_basic_set_dim(bset, isl_dim_set);310	nparam = isl_basic_set_dim(bset, isl_dim_param);311	total = isl_basic_set_dim(bset, isl_dim_all);312	if (n_div < 0 || d < 0 || nparam < 0 || total < 0)313		return isl_bool_error;314 315	bset = isl_basic_set_copy(bset);316	bset = isl_basic_set_cow(bset);317	bset = isl_basic_set_extend_constraints(bset, 0, 1);318	k = isl_basic_set_alloc_inequality(bset);319	if (k < 0)320		goto error;321	isl_seq_clr(bset->ineq[k], 1 + total);322	isl_seq_cpy(bset->ineq[k], c, 1 + nparam);323	for (i = 0; i < n_div; ++i) {324		if (div_purity[i] != PURE_PARAM)325			continue;326		isl_int_set(bset->ineq[k][1 + nparam + d + i],327			    c[1 + nparam + d + i]);328	}329	isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);330	empty = isl_basic_set_is_empty(bset);331	isl_basic_set_free(bset);332 333	return empty;334error:335	isl_basic_set_free(bset);336	return isl_bool_error;337}338 339/* Return PURE_PARAM if only the coefficients of the parameters are non-zero.340 * Return PURE_VAR if only the coefficients of the set variables are non-zero.341 * Return MIXED if only the coefficients of the parameters and the set342 * 	variables are non-zero and if moreover the parametric constant343 * 	can never attain positive values.344 * Return IMPURE otherwise.345 */346static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,347	int eq)348{349	isl_size d;350	isl_size n_div;351	isl_size nparam;352	isl_bool empty;353	int i;354	int p = 0, v = 0;355 356	n_div = isl_basic_set_dim(bset, isl_dim_div);357	d = isl_basic_set_dim(bset, isl_dim_set);358	nparam = isl_basic_set_dim(bset, isl_dim_param);359	if (n_div < 0 || d < 0 || nparam < 0)360		return -1;361 362	for (i = 0; i < n_div; ++i) {363		if (isl_int_is_zero(c[1 + nparam + d + i]))364			continue;365		switch (div_purity[i]) {366		case PURE_PARAM: p = 1; break;367		case PURE_VAR: v = 1; break;368		default: return IMPURE;369		}370	}371	if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1)372		return PURE_VAR;373	if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)374		return PURE_PARAM;375 376	empty = parametric_constant_never_positive(bset, c, div_purity);377	if (eq && empty >= 0 && !empty) {378		isl_seq_neg(c, c, 1 + nparam + d + n_div);379		empty = parametric_constant_never_positive(bset, c, div_purity);380	}381 382	return empty < 0 ? -1 : empty ? MIXED : IMPURE;383}384 385/* Return an array of integers indicating the type of each div in bset.386 * If the div is (recursively) defined in terms of only the parameters,387 * then the type is PURE_PARAM.388 * If the div is (recursively) defined in terms of only the set variables,389 * then the type is PURE_VAR.390 * Otherwise, the type is IMPURE.391 */392static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)393{394	int i, j;395	int *div_purity;396	isl_size d;397	isl_size n_div;398	isl_size nparam;399 400	n_div = isl_basic_set_dim(bset, isl_dim_div);401	d = isl_basic_set_dim(bset, isl_dim_set);402	nparam = isl_basic_set_dim(bset, isl_dim_param);403	if (n_div < 0 || d < 0 || nparam < 0)404		return NULL;405 406	div_purity = isl_alloc_array(bset->ctx, int, n_div);407	if (n_div && !div_purity)408		return NULL;409 410	for (i = 0; i < bset->n_div; ++i) {411		int p = 0, v = 0;412		if (isl_int_is_zero(bset->div[i][0])) {413			div_purity[i] = IMPURE;414			continue;415		}416		if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)417			p = 1;418		if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)419			v = 1;420		for (j = 0; j < i; ++j) {421			if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))422				continue;423			switch (div_purity[j]) {424			case PURE_PARAM: p = 1; break;425			case PURE_VAR: v = 1; break;426			default: p = v = 1; break;427			}428		}429		div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM;430	}431 432	return div_purity;433}434 435/* Given a path with the as yet unconstrained length at div position "pos",436 * check if setting the length to zero results in only the identity437 * mapping.438 */439static isl_bool empty_path_is_identity(__isl_keep isl_basic_map *path,440	unsigned pos)441{442	isl_basic_map *test = NULL;443	isl_basic_map *id = NULL;444	isl_bool is_id;445 446	test = isl_basic_map_copy(path);447	test = isl_basic_map_fix_si(test, isl_dim_div, pos, 0);448	id = isl_basic_map_identity(isl_basic_map_get_space(path));449	is_id = isl_basic_map_is_equal(test, id);450	isl_basic_map_free(test);451	isl_basic_map_free(id);452	return is_id;453}454 455/* If any of the constraints is found to be impure then this function456 * sets *impurity to 1.457 *458 * If impurity is NULL then we are dealing with a non-parametric set459 * and so the constraints are obviously PURE_VAR.460 */461static __isl_give isl_basic_map *add_delta_constraints(462	__isl_take isl_basic_map *path,463	__isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,464	unsigned d, int *div_purity, int eq, int *impurity)465{466	int i, k;467	int n = eq ? delta->n_eq : delta->n_ineq;468	isl_int **delta_c = eq ? delta->eq : delta->ineq;469	isl_size n_div, total;470 471	n_div = isl_basic_set_dim(delta, isl_dim_div);472	total = isl_basic_map_dim(path, isl_dim_all);473	if (n_div < 0 || total < 0)474		return isl_basic_map_free(path);475 476	for (i = 0; i < n; ++i) {477		isl_int *path_c;478		int p = PURE_VAR;479		if (impurity)480			p = purity(delta, delta_c[i], div_purity, eq);481		if (p < 0)482			goto error;483		if (p != PURE_VAR && p != PURE_PARAM && !*impurity)484			*impurity = 1;485		if (p == IMPURE)486			continue;487		if (eq && p != MIXED) {488			k = isl_basic_map_alloc_equality(path);489			if (k < 0)490				goto error;491			path_c = path->eq[k];492		} else {493			k = isl_basic_map_alloc_inequality(path);494			if (k < 0)495				goto error;496			path_c = path->ineq[k];497		}498		isl_seq_clr(path_c, 1 + total);499		if (p == PURE_VAR) {500			isl_seq_cpy(path_c + off,501				    delta_c[i] + 1 + nparam, d);502			isl_int_set(path_c[off + d], delta_c[i][0]);503		} else if (p == PURE_PARAM) {504			isl_seq_cpy(path_c, delta_c[i], 1 + nparam);505		} else {506			isl_seq_cpy(path_c + off,507				    delta_c[i] + 1 + nparam, d);508			isl_seq_cpy(path_c, delta_c[i], 1 + nparam);509		}510		isl_seq_cpy(path_c + off - n_div,511			    delta_c[i] + 1 + nparam + d, n_div);512	}513 514	return path;515error:516	isl_basic_map_free(path);517	return NULL;518}519 520/* Given a set of offsets "delta", construct a relation of the521 * given dimension specification (Z^{n+1} -> Z^{n+1}) that522 * is an overapproximation of the relations that523 * maps an element x to any element that can be reached524 * by taking a non-negative number of steps along any of525 * the elements in "delta".526 * That is, construct an approximation of527 *528 *	{ [x] -> [y] : exists f \in \delta, k \in Z :529 *					y = x + k [f, 1] and k >= 0 }530 *531 * For any element in this relation, the number of steps taken532 * is equal to the difference in the final coordinates.533 *534 * In particular, let delta be defined as535 *536 *	\delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and537 *				C x + C'p + c >= 0 and538 *				D x + D'p + d >= 0 }539 *540 * where the constraints C x + C'p + c >= 0 are such that the parametric541 * constant term of each constraint j, "C_j x + C'_j p + c_j",542 * can never attain positive values, then the relation is constructed as543 *544 *	{ [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and545 *			A f + k a >= 0 and B p + b >= 0 and546 *			C f + C'p + c >= 0 and k >= 1 }547 *	union { [x] -> [x] }548 *549 * If the zero-length paths happen to correspond exactly to the identity550 * mapping, then we return551 *552 *	{ [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and553 *			A f + k a >= 0 and B p + b >= 0 and554 *			C f + C'p + c >= 0 and k >= 0 }555 *556 * instead.557 *558 * Existentially quantified variables in \delta are handled by559 * classifying them as independent of the parameters, purely560 * parameter dependent and others.  Constraints containing561 * any of the other existentially quantified variables are removed.562 * This is safe, but leads to an additional overapproximation.563 *564 * If there are any impure constraints, then we also eliminate565 * the parameters from \delta, resulting in a set566 *567 *	\delta' = { [x] : E x + e >= 0 }568 *569 * and add the constraints570 *571 *			E f + k e >= 0572 *573 * to the constructed relation.574 */575static __isl_give isl_map *path_along_delta(__isl_take isl_space *space,576	__isl_take isl_basic_set *delta)577{578	isl_basic_map *path = NULL;579	isl_size d;580	isl_size n_div;581	isl_size nparam;582	isl_size total;583	unsigned off;584	int i, k;585	isl_bool is_id;586	int *div_purity = NULL;587	int impurity = 0;588 589	n_div = isl_basic_set_dim(delta, isl_dim_div);590	d = isl_basic_set_dim(delta, isl_dim_set);591	nparam = isl_basic_set_dim(delta, isl_dim_param);592	if (n_div < 0 || d < 0 || nparam < 0)593		goto error;594	path = isl_basic_map_alloc_space(isl_space_copy(space), n_div + d + 1,595			d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);596	off = 1 + nparam + 2 * (d + 1) + n_div;597 598	for (i = 0; i < n_div + d + 1; ++i) {599		k = isl_basic_map_alloc_div(path);600		if (k < 0)601			goto error;602		isl_int_set_si(path->div[k][0], 0);603	}604 605	total = isl_basic_map_dim(path, isl_dim_all);606	if (total < 0)607		goto error;608	for (i = 0; i < d + 1; ++i) {609		k = isl_basic_map_alloc_equality(path);610		if (k < 0)611			goto error;612		isl_seq_clr(path->eq[k], 1 + total);613		isl_int_set_si(path->eq[k][1 + nparam + i], 1);614		isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);615		isl_int_set_si(path->eq[k][off + i], 1);616	}617 618	div_purity = get_div_purity(delta);619	if (n_div && !div_purity)620		goto error;621 622	path = add_delta_constraints(path, delta, off, nparam, d,623				     div_purity, 1, &impurity);624	path = add_delta_constraints(path, delta, off, nparam, d,625				     div_purity, 0, &impurity);626	if (impurity) {627		isl_space *space = isl_basic_set_get_space(delta);628		delta = isl_basic_set_project_out(delta,629						  isl_dim_param, 0, nparam);630		delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam);631		delta = isl_basic_set_reset_space(delta, space);632		if (!delta)633			goto error;634		path = isl_basic_map_extend_constraints(path, delta->n_eq,635							delta->n_ineq + 1);636		path = add_delta_constraints(path, delta, off, nparam, d,637					     NULL, 1, NULL);638		path = add_delta_constraints(path, delta, off, nparam, d,639					     NULL, 0, NULL);640		path = isl_basic_map_gauss(path, NULL);641	}642 643	is_id = empty_path_is_identity(path, n_div + d);644	if (is_id < 0)645		goto error;646 647	k = isl_basic_map_alloc_inequality(path);648	if (k < 0)649		goto error;650	isl_seq_clr(path->ineq[k], 1 + total);651	if (!is_id)652		isl_int_set_si(path->ineq[k][0], -1);653	isl_int_set_si(path->ineq[k][off + d], 1);654			655	free(div_purity);656	isl_basic_set_free(delta);657	path = isl_basic_map_finalize(path);658	if (is_id) {659		isl_space_free(space);660		return isl_map_from_basic_map(path);661	}662	return isl_basic_map_union(path, isl_basic_map_identity(space));663error:664	free(div_purity);665	isl_space_free(space);666	isl_basic_set_free(delta);667	isl_basic_map_free(path);668	return NULL;669}670 671/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param",672 * construct a map that equates the parameter to the difference673 * in the final coordinates and imposes that this difference is positive.674 * That is, construct675 *676 *	{ [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }677 */678static __isl_give isl_map *equate_parameter_to_length(679	__isl_take isl_space *space, unsigned param)680{681	struct isl_basic_map *bmap;682	isl_size d;683	isl_size nparam;684	isl_size total;685	int k;686 687	d = isl_space_dim(space, isl_dim_in);688	nparam = isl_space_dim(space, isl_dim_param);689	total = isl_space_dim(space, isl_dim_all);690	if (d < 0 || nparam < 0 || total < 0)691		space = isl_space_free(space);692	bmap = isl_basic_map_alloc_space(space, 0, 1, 1);693	k = isl_basic_map_alloc_equality(bmap);694	if (k < 0)695		goto error;696	isl_seq_clr(bmap->eq[k], 1 + total);697	isl_int_set_si(bmap->eq[k][1 + param], -1);698	isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);699	isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);700 701	k = isl_basic_map_alloc_inequality(bmap);702	if (k < 0)703		goto error;704	isl_seq_clr(bmap->ineq[k], 1 + total);705	isl_int_set_si(bmap->ineq[k][1 + param], 1);706	isl_int_set_si(bmap->ineq[k][0], -1);707 708	bmap = isl_basic_map_finalize(bmap);709	return isl_map_from_basic_map(bmap);710error:711	isl_basic_map_free(bmap);712	return NULL;713}714 715/* Check whether "path" is acyclic, where the last coordinates of domain716 * and range of path encode the number of steps taken.717 * That is, check whether718 *719 *	{ d | d = y - x and (x,y) in path }720 *721 * does not contain any element with positive last coordinate (positive length)722 * and zero remaining coordinates (cycle).723 */724static isl_bool is_acyclic(__isl_take isl_map *path)725{726	int i;727	isl_bool acyclic;728	isl_size dim;729	struct isl_set *delta;730 731	delta = isl_map_deltas(path);732	dim = isl_set_dim(delta, isl_dim_set);733	if (dim < 0)734		delta = isl_set_free(delta);735	for (i = 0; i < dim; ++i) {736		if (i == dim -1)737			delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);738		else739			delta = isl_set_fix_si(delta, isl_dim_set, i, 0);740	}741 742	acyclic = isl_set_is_empty(delta);743	isl_set_free(delta);744 745	return acyclic;746}747 748/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D749 * and a dimension specification (Z^{n+1} -> Z^{n+1}),750 * construct a map that is an overapproximation of the map751 * that takes an element from the space D \times Z to another752 * element from the same space, such that the first n coordinates of the753 * difference between them is a sum of differences between images754 * and pre-images in one of the R_i and such that the last coordinate755 * is equal to the number of steps taken.756 * That is, let757 *758 *	\Delta_i = { y - x | (x, y) in R_i }759 *760 * then the constructed map is an overapproximation of761 *762 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :763 *				d = (\sum_i k_i \delta_i, \sum_i k_i) }764 *765 * The elements of the singleton \Delta_i's are collected as the766 * rows of the steps matrix.  For all these \Delta_i's together,767 * a single path is constructed.768 * For each of the other \Delta_i's, we compute an overapproximation769 * of the paths along elements of \Delta_i.770 * Since each of these paths performs an addition, composition is771 * symmetric and we can simply compose all resulting paths in any order.772 */773static __isl_give isl_map *construct_extended_path(__isl_take isl_space *space,774	__isl_keep isl_map *map, int *project)775{776	struct isl_mat *steps = NULL;777	struct isl_map *path = NULL;778	isl_size d;779	int i, j, n;780 781	d = isl_map_dim(map, isl_dim_in);782	if (d < 0)783		goto error;784 785	path = isl_map_identity(isl_space_copy(space));786 787	steps = isl_mat_alloc(map->ctx, map->n, d);788	if (!steps)789		goto error;790 791	n = 0;792	for (i = 0; i < map->n; ++i) {793		struct isl_basic_set *delta;794 795		delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));796 797		for (j = 0; j < d; ++j) {798			isl_bool fixed;799 800			fixed = isl_basic_set_plain_dim_is_fixed(delta, j,801							    &steps->row[n][j]);802			if (fixed < 0) {803				isl_basic_set_free(delta);804				goto error;805			}806			if (!fixed)807				break;808		}809 810 811		if (j < d) {812			path = isl_map_apply_range(path,813				path_along_delta(isl_space_copy(space), delta));814			path = isl_map_coalesce(path);815		} else {816			isl_basic_set_free(delta);817			++n;818		}819	}820 821	if (n > 0) {822		steps->n_row = n;823		path = isl_map_apply_range(path,824				path_along_steps(isl_space_copy(space), steps));825	}826 827	if (project && *project) {828		*project = is_acyclic(isl_map_copy(path));829		if (*project < 0)830			goto error;831	}832 833	isl_space_free(space);834	isl_mat_free(steps);835	return path;836error:837	isl_space_free(space);838	isl_mat_free(steps);839	isl_map_free(path);840	return NULL;841}842 843static isl_bool isl_set_overlaps(__isl_keep isl_set *set1,844	__isl_keep isl_set *set2)845{846	return isl_bool_not(isl_set_is_disjoint(set1, set2));847}848 849/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D850 * and a dimension specification (Z^{n+1} -> Z^{n+1}),851 * construct a map that is an overapproximation of the map852 * that takes an element from the dom R \times Z to an853 * element from ran R \times Z, such that the first n coordinates of the854 * difference between them is a sum of differences between images855 * and pre-images in one of the R_i and such that the last coordinate856 * is equal to the number of steps taken.857 * That is, let858 *859 *	\Delta_i = { y - x | (x, y) in R_i }860 *861 * then the constructed map is an overapproximation of862 *863 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :864 *				d = (\sum_i k_i \delta_i, \sum_i k_i) and865 *				x in dom R and x + d in ran R and866 *				\sum_i k_i >= 1 }867 */868static __isl_give isl_map *construct_component(__isl_take isl_space *space,869	__isl_keep isl_map *map, isl_bool *exact, int project)870{871	struct isl_set *domain = NULL;872	struct isl_set *range = NULL;873	struct isl_map *app = NULL;874	struct isl_map *path = NULL;875	isl_bool overlaps;876	int check;877 878	domain = isl_map_domain(isl_map_copy(map));879	domain = isl_set_coalesce(domain);880	range = isl_map_range(isl_map_copy(map));881	range = isl_set_coalesce(range);882	overlaps = isl_set_overlaps(domain, range);883	if (overlaps < 0 || !overlaps) {884		isl_set_free(domain);885		isl_set_free(range);886		isl_space_free(space);887 888		if (overlaps < 0)889			map = NULL;890		map = isl_map_copy(map);891		map = isl_map_add_dims(map, isl_dim_in, 1);892		map = isl_map_add_dims(map, isl_dim_out, 1);893		map = set_path_length(map, 1, 1);894		return map;895	}896	app = isl_map_from_domain_and_range(domain, range);897	app = isl_map_add_dims(app, isl_dim_in, 1);898	app = isl_map_add_dims(app, isl_dim_out, 1);899 900	check = exact && *exact == isl_bool_true;901	path = construct_extended_path(isl_space_copy(space), map,902					check ? &project : NULL);903	app = isl_map_intersect(app, path);904 905	if (check &&906	    (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),907				      project)) < 0)908		goto error;909 910	isl_space_free(space);911	app = set_path_length(app, 0, 1);912	return app;913error:914	isl_space_free(space);915	isl_map_free(app);916	return NULL;917}918 919/* Call construct_component and, if "project" is set, project out920 * the final coordinates.921 */922static __isl_give isl_map *construct_projected_component(923	__isl_take isl_space *space,924	__isl_keep isl_map *map, isl_bool *exact, int project)925{926	isl_map *app;927	unsigned d;928 929	if (!space)930		return NULL;931	d = isl_space_dim(space, isl_dim_in);932 933	app = construct_component(space, map, exact, project);934	if (project) {935		app = isl_map_project_out(app, isl_dim_in, d - 1, 1);936		app = isl_map_project_out(app, isl_dim_out, d - 1, 1);937	}938	return app;939}940 941/* Compute an extended version, i.e., with path lengths, of942 * an overapproximation of the transitive closure of "bmap"943 * with path lengths greater than or equal to zero and with944 * domain and range equal to "dom".945 */946static __isl_give isl_map *q_closure(__isl_take isl_space *space,947	__isl_take isl_set *dom, __isl_keep isl_basic_map *bmap,948	isl_bool *exact)949{950	int project = 1;951	isl_map *path;952	isl_map *map;953	isl_map *app;954 955	dom = isl_set_add_dims(dom, isl_dim_set, 1);956	app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));957	map = isl_map_from_basic_map(isl_basic_map_copy(bmap));958	path = construct_extended_path(space, map, &project);959	app = isl_map_intersect(app, path);960 961	if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)962		goto error;963 964	return app;965error:966	isl_map_free(app);967	return NULL;968}969 970/* Check whether qc has any elements of length at least one971 * with domain and/or range outside of dom and ran.972 */973static isl_bool has_spurious_elements(__isl_keep isl_map *qc,974	__isl_keep isl_set *dom, __isl_keep isl_set *ran)975{976	isl_set *s;977	isl_bool subset;978	isl_size d;979 980	d = isl_map_dim(qc, isl_dim_in);981	if (d < 0 || !dom || !ran)982		return isl_bool_error;983 984	qc = isl_map_copy(qc);985	qc = set_path_length(qc, 0, 1);986	qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);987	qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);988 989	s = isl_map_domain(isl_map_copy(qc));990	subset = isl_set_is_subset(s, dom);991	isl_set_free(s);992	if (subset < 0)993		goto error;994	if (!subset) {995		isl_map_free(qc);996		return isl_bool_true;997	}998 999	s = isl_map_range(qc);1000	subset = isl_set_is_subset(s, ran);1001	isl_set_free(s);1002 1003	return isl_bool_not(subset);1004error:1005	isl_map_free(qc);1006	return isl_bool_error;1007}1008 1009#define LEFT	21010#define RIGHT	11011 1012/* For each basic map in "map", except i, check whether it combines1013 * with the transitive closure that is reflexive on C combines1014 * to the left and to the right.1015 *1016 * In particular, if1017 *1018 *	dom map_j \subseteq C1019 *1020 * then right[j] is set to 1.  Otherwise, if1021 *1022 *	ran map_i \cap dom map_j = \emptyset1023 *1024 * then right[j] is set to 0.  Otherwise, composing to the right1025 * is impossible.1026 *1027 * Similar, for composing to the left, we have if1028 *1029 *	ran map_j \subseteq C1030 *1031 * then left[j] is set to 1.  Otherwise, if1032 *1033 *	dom map_i \cap ran map_j = \emptyset1034 *1035 * then left[j] is set to 0.  Otherwise, composing to the left1036 * is impossible.1037 *1038 * The return value is or'd with LEFT if composing to the left1039 * is possible and with RIGHT if composing to the right is possible.1040 */1041static int composability(__isl_keep isl_set *C, int i,1042	isl_set **dom, isl_set **ran, int *left, int *right,1043	__isl_keep isl_map *map)1044{1045	int j;1046	int ok;1047 1048	ok = LEFT | RIGHT;1049	for (j = 0; j < map->n && ok; ++j) {1050		isl_bool overlaps, subset;1051		if (j == i)1052			continue;1053 1054		if (ok & RIGHT) {1055			if (!dom[j])1056				dom[j] = isl_set_from_basic_set(1057					isl_basic_map_domain(1058						isl_basic_map_copy(map->p[j])));1059			if (!dom[j])1060				return -1;1061			overlaps = isl_set_overlaps(ran[i], dom[j]);1062			if (overlaps < 0)1063				return -1;1064			if (!overlaps)1065				right[j] = 0;1066			else {1067				subset = isl_set_is_subset(dom[j], C);1068				if (subset < 0)1069					return -1;1070				if (subset)1071					right[j] = 1;1072				else1073					ok &= ~RIGHT;1074			}1075		}1076 1077		if (ok & LEFT) {1078			if (!ran[j])1079				ran[j] = isl_set_from_basic_set(1080					isl_basic_map_range(1081						isl_basic_map_copy(map->p[j])));1082			if (!ran[j])1083				return -1;1084			overlaps = isl_set_overlaps(dom[i], ran[j]);1085			if (overlaps < 0)1086				return -1;1087			if (!overlaps)1088				left[j] = 0;1089			else {1090				subset = isl_set_is_subset(ran[j], C);1091				if (subset < 0)1092					return -1;1093				if (subset)1094					left[j] = 1;1095				else1096					ok &= ~LEFT;1097			}1098		}1099	}1100 1101	return ok;1102}1103 1104static __isl_give isl_map *anonymize(__isl_take isl_map *map)1105{1106	map = isl_map_reset(map, isl_dim_in);1107	map = isl_map_reset(map, isl_dim_out);1108	return map;1109}1110 1111/* Return a map that is a union of the basic maps in "map", except i,1112 * composed to left and right with qc based on the entries of "left"1113 * and "right".1114 */1115static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,1116	__isl_take isl_map *qc, int *left, int *right)1117{1118	int j;1119	isl_map *comp;1120 1121	comp = isl_map_empty(isl_map_get_space(map));1122	for (j = 0; j < map->n; ++j) {1123		isl_map *map_j;1124 1125		if (j == i)1126			continue;1127 1128		map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));1129		map_j = anonymize(map_j);1130		if (left && left[j])1131			map_j = isl_map_apply_range(map_j, isl_map_copy(qc));1132		if (right && right[j])1133			map_j = isl_map_apply_range(isl_map_copy(qc), map_j);1134		comp = isl_map_union(comp, map_j);1135	}1136 1137	comp = isl_map_compute_divs(comp);1138	comp = isl_map_coalesce(comp);1139 1140	isl_map_free(qc);1141 1142	return comp;1143}1144 1145/* Compute the transitive closure of "map" incrementally by1146 * computing1147 *1148 *	map_i^+ \cup qc^+1149 *1150 * or1151 *1152 *	map_i^+ \cup ((id \cup map_i^) \circ qc^+)1153 *1154 * or1155 *1156 *	map_i^+ \cup (qc^+ \circ (id \cup map_i^))1157 *1158 * depending on whether left or right are NULL.1159 */1160static __isl_give isl_map *compute_incremental(1161	__isl_take isl_space *space, __isl_keep isl_map *map,1162	int i, __isl_take isl_map *qc, int *left, int *right, isl_bool *exact)1163{1164	isl_map *map_i;1165	isl_map *tc;1166	isl_map *rtc = NULL;1167 1168	if (!map)1169		goto error;1170	isl_assert(map->ctx, left || right, goto error);1171 1172	map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));1173	tc = construct_projected_component(isl_space_copy(space), map_i,1174						exact, 1);1175	isl_map_free(map_i);1176 1177	if (*exact)1178		qc = isl_map_transitive_closure(qc, exact);1179 1180	if (!*exact) {1181		isl_space_free(space);1182		isl_map_free(tc);1183		isl_map_free(qc);1184		return isl_map_universe(isl_map_get_space(map));1185	}1186 1187	if (!left || !right)1188		rtc = isl_map_union(isl_map_copy(tc),1189				    isl_map_identity(isl_map_get_space(tc)));1190	if (!right)1191		qc = isl_map_apply_range(rtc, qc);1192	if (!left)1193		qc = isl_map_apply_range(qc, rtc);1194	qc = isl_map_union(tc, qc);1195 1196	isl_space_free(space);1197 1198	return qc;1199error:1200	isl_space_free(space);1201	isl_map_free(qc);1202	return NULL;1203}1204 1205/* Given a map "map", try to find a basic map such that1206 * map^+ can be computed as1207 *1208 * map^+ = map_i^+ \cup1209 *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+1210 *1211 * with C the simple hull of the domain and range of the input map.1212 * map_i^ \cup Id_C is computed by allowing the path lengths to be zero1213 * and by intersecting domain and range with C.1214 * Of course, we need to check that this is actually equal to map_i^ \cup Id_C.1215 * Also, we only use the incremental computation if all the transitive1216 * closures are exact and if the number of basic maps in the union,1217 * after computing the integer divisions, is smaller than the number1218 * of basic maps in the input map.1219 */1220static isl_bool incremental_on_entire_domain(__isl_keep isl_space *space,1221	__isl_keep isl_map *map,1222	isl_set **dom, isl_set **ran, int *left, int *right,1223	__isl_give isl_map **res)1224{1225	int i;1226	isl_set *C;1227	isl_size d;1228 1229	*res = NULL;1230 1231	d = isl_map_dim(map, isl_dim_in);1232	if (d < 0)1233		return isl_bool_error;1234 1235	C = isl_set_union(isl_map_domain(isl_map_copy(map)),1236			  isl_map_range(isl_map_copy(map)));1237	C = isl_set_from_basic_set(isl_set_simple_hull(C));1238	if (!C)1239		return isl_bool_error;1240	if (C->n != 1) {1241		isl_set_free(C);1242		return isl_bool_false;1243	}1244 1245	for (i = 0; i < map->n; ++i) {1246		isl_map *qc;1247		isl_bool exact_i;1248		isl_bool spurious;1249		int j;1250		dom[i] = isl_set_from_basic_set(isl_basic_map_domain(1251					isl_basic_map_copy(map->p[i])));1252		ran[i] = isl_set_from_basic_set(isl_basic_map_range(1253					isl_basic_map_copy(map->p[i])));1254		qc = q_closure(isl_space_copy(space), isl_set_copy(C),1255				map->p[i], &exact_i);1256		if (!qc)1257			goto error;1258		if (!exact_i) {1259			isl_map_free(qc);1260			continue;1261		}1262		spurious = has_spurious_elements(qc, dom[i], ran[i]);1263		if (spurious) {1264			isl_map_free(qc);1265			if (spurious < 0)1266				goto error;1267			continue;1268		}1269		qc = isl_map_project_out(qc, isl_dim_in, d, 1);1270		qc = isl_map_project_out(qc, isl_dim_out, d, 1);1271		qc = isl_map_compute_divs(qc);1272		for (j = 0; j < map->n; ++j)1273			left[j] = right[j] = 1;1274		qc = compose(map, i, qc, left, right);1275		if (!qc)1276			goto error;1277		if (qc->n >= map->n) {1278			isl_map_free(qc);1279			continue;1280		}1281		*res = compute_incremental(isl_space_copy(space), map, i, qc,1282				left, right, &exact_i);1283		if (!*res)1284			goto error;1285		if (exact_i)1286			break;1287		isl_map_free(*res);1288		*res = NULL;1289	}1290 1291	isl_set_free(C);1292 1293	return isl_bool_ok(*res != NULL);1294error:1295	isl_set_free(C);1296	return isl_bool_error;1297}1298 1299/* Try and compute the transitive closure of "map" as1300 *1301 * map^+ = map_i^+ \cup1302 *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+1303 *1304 * with C either the simple hull of the domain and range of the entire1305 * map or the simple hull of domain and range of map_i.1306 */1307static __isl_give isl_map *incremental_closure(__isl_take isl_space *space,1308	__isl_keep isl_map *map, isl_bool *exact, int project)1309{1310	int i;1311	isl_set **dom = NULL;1312	isl_set **ran = NULL;1313	int *left = NULL;1314	int *right = NULL;1315	isl_set *C;1316	isl_size d;1317	isl_map *res = NULL;1318 1319	if (!project)1320		return construct_projected_component(space, map, exact,1321							project);1322 1323	if (!map)1324		goto error;1325	if (map->n <= 1)1326		return construct_projected_component(space, map, exact,1327							project);1328 1329	d = isl_map_dim(map, isl_dim_in);1330	if (d < 0)1331		goto error;1332 1333	dom = isl_calloc_array(map->ctx, isl_set *, map->n);1334	ran = isl_calloc_array(map->ctx, isl_set *, map->n);1335	left = isl_calloc_array(map->ctx, int, map->n);1336	right = isl_calloc_array(map->ctx, int, map->n);1337	if (!ran || !dom || !left || !right)1338		goto error;1339 1340	if (incremental_on_entire_domain(space, map, dom, ran, left, right,1341					&res) < 0)1342		goto error;1343 1344	for (i = 0; !res && i < map->n; ++i) {1345		isl_map *qc;1346		int comp;1347		isl_bool exact_i, spurious;1348		if (!dom[i])1349			dom[i] = isl_set_from_basic_set(1350					isl_basic_map_domain(1351						isl_basic_map_copy(map->p[i])));1352		if (!dom[i])1353			goto error;1354		if (!ran[i])1355			ran[i] = isl_set_from_basic_set(1356					isl_basic_map_range(1357						isl_basic_map_copy(map->p[i])));1358		if (!ran[i])1359			goto error;1360		C = isl_set_union(isl_set_copy(dom[i]),1361				      isl_set_copy(ran[i]));1362		C = isl_set_from_basic_set(isl_set_simple_hull(C));1363		if (!C)1364			goto error;1365		if (C->n != 1) {1366			isl_set_free(C);1367			continue;1368		}1369		comp = composability(C, i, dom, ran, left, right, map);1370		if (!comp || comp < 0) {1371			isl_set_free(C);1372			if (comp < 0)1373				goto error;1374			continue;1375		}1376		qc = q_closure(isl_space_copy(space), C, map->p[i], &exact_i);1377		if (!qc)1378			goto error;1379		if (!exact_i) {1380			isl_map_free(qc);1381			continue;1382		}1383		spurious = has_spurious_elements(qc, dom[i], ran[i]);1384		if (spurious) {1385			isl_map_free(qc);1386			if (spurious < 0)1387				goto error;1388			continue;1389		}1390		qc = isl_map_project_out(qc, isl_dim_in, d, 1);1391		qc = isl_map_project_out(qc, isl_dim_out, d, 1);1392		qc = isl_map_compute_divs(qc);1393		qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,1394				(comp & RIGHT) ? right : NULL);1395		if (!qc)1396			goto error;1397		if (qc->n >= map->n) {1398			isl_map_free(qc);1399			continue;1400		}1401		res = compute_incremental(isl_space_copy(space), map, i, qc,1402				(comp & LEFT) ? left : NULL,1403				(comp & RIGHT) ? right : NULL, &exact_i);1404		if (!res)1405			goto error;1406		if (exact_i)1407			break;1408		isl_map_free(res);1409		res = NULL;1410	}1411 1412	for (i = 0; i < map->n; ++i) {1413		isl_set_free(dom[i]);1414		isl_set_free(ran[i]);1415	}1416	free(dom);1417	free(ran);1418	free(left);1419	free(right);1420 1421	if (res) {1422		isl_space_free(space);1423		return res;1424	}1425 1426	return construct_projected_component(space, map, exact, project);1427error:1428	if (dom)1429		for (i = 0; i < map->n; ++i)1430			isl_set_free(dom[i]);1431	free(dom);1432	if (ran)1433		for (i = 0; i < map->n; ++i)1434			isl_set_free(ran[i]);1435	free(ran);1436	free(left);1437	free(right);1438	isl_space_free(space);1439	return NULL;1440}1441 1442/* Given an array of sets "set", add "dom" at position "pos"1443 * and search for elements at earlier positions that overlap with "dom".1444 * If any can be found, then merge all of them, together with "dom", into1445 * a single set and assign the union to the first in the array,1446 * which becomes the new group leader for all groups involved in the merge.1447 * During the search, we only consider group leaders, i.e., those with1448 * group[i] = i, as the other sets have already been combined1449 * with one of the group leaders.1450 */1451static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)1452{1453	int i;1454 1455	group[pos] = pos;1456	set[pos] = isl_set_copy(dom);1457 1458	for (i = pos - 1; i >= 0; --i) {1459		isl_bool o;1460 1461		if (group[i] != i)1462			continue;1463 1464		o = isl_set_overlaps(set[i], dom);1465		if (o < 0)1466			goto error;1467		if (!o)1468			continue;1469 1470		set[i] = isl_set_union(set[i], set[group[pos]]);1471		set[group[pos]] = NULL;1472		if (!set[i])1473			goto error;1474		group[group[pos]] = i;1475		group[pos] = i;1476	}1477 1478	isl_set_free(dom);1479	return 0;1480error:1481	isl_set_free(dom);1482	return -1;1483}1484 1485/* Construct a map [x] -> [x+1], with parameters prescribed by "space".1486 */1487static __isl_give isl_map *increment(__isl_take isl_space *space)1488{1489	int k;1490	isl_basic_map *bmap;1491	isl_size total;1492 1493	space = isl_space_set_from_params(space);1494	space = isl_space_add_dims(space, isl_dim_set, 1);1495	space = isl_space_map_from_set(space);1496	bmap = isl_basic_map_alloc_space(space, 0, 1, 0);1497	total = isl_basic_map_dim(bmap, isl_dim_all);1498	k = isl_basic_map_alloc_equality(bmap);1499	if (total < 0 || k < 0)1500		goto error;1501	isl_seq_clr(bmap->eq[k], 1 + total);1502	isl_int_set_si(bmap->eq[k][0], 1);1503	isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);1504	isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);1505	return isl_map_from_basic_map(bmap);1506error:1507	isl_basic_map_free(bmap);1508	return NULL;1509}1510 1511/* Replace each entry in the n by n grid of maps by the cross product1512 * with the relation { [i] -> [i + 1] }.1513 */1514static isl_stat add_length(__isl_keep isl_map *map, isl_map ***grid, int n)1515{1516	int i, j;1517	isl_space *space;1518	isl_map *step;1519 1520	space = isl_space_params(isl_map_get_space(map));1521	step = increment(space);1522 1523	if (!step)1524		return isl_stat_error;1525 1526	for (i = 0; i < n; ++i)1527		for (j = 0; j < n; ++j)1528			grid[i][j] = isl_map_product(grid[i][j],1529						     isl_map_copy(step));1530 1531	isl_map_free(step);1532 1533	return isl_stat_ok;1534}1535 1536/* The core of the Floyd-Warshall algorithm.1537 * Updates the given n x x matrix of relations in place.1538 *1539 * The algorithm iterates over all vertices.  In each step, the whole1540 * matrix is updated to include all paths that go to the current vertex,1541 * possibly stay there a while (including passing through earlier vertices)1542 * and then come back.  At the start of each iteration, the diagonal1543 * element corresponding to the current vertex is replaced by its1544 * transitive closure to account for all indirect paths that stay1545 * in the current vertex.1546 */1547static void floyd_warshall_iterate(isl_map ***grid, int n, isl_bool *exact)1548{1549	int r, p, q;1550 1551	for (r = 0; r < n; ++r) {1552		isl_bool r_exact;1553		int check = exact && *exact == isl_bool_true;1554		grid[r][r] = isl_map_transitive_closure(grid[r][r],1555				check ? &r_exact : NULL);1556		if (check && !r_exact)1557			*exact = isl_bool_false;1558 1559		for (p = 0; p < n; ++p)1560			for (q = 0; q < n; ++q) {1561				isl_map *loop;1562				if (p == r && q == r)1563					continue;1564				loop = isl_map_apply_range(1565						isl_map_copy(grid[p][r]),1566						isl_map_copy(grid[r][q]));1567				grid[p][q] = isl_map_union(grid[p][q], loop);1568				loop = isl_map_apply_range(1569						isl_map_copy(grid[p][r]),1570					isl_map_apply_range(1571						isl_map_copy(grid[r][r]),1572						isl_map_copy(grid[r][q])));1573				grid[p][q] = isl_map_union(grid[p][q], loop);1574				grid[p][q] = isl_map_coalesce(grid[p][q]);1575			}1576	}1577}1578 1579/* Given a partition of the domains and ranges of the basic maps in "map",1580 * apply the Floyd-Warshall algorithm with the elements in the partition1581 * as vertices.1582 *1583 * In particular, there are "n" elements in the partition and "group" is1584 * an array of length 2 * map->n with entries in [0,n-1].1585 *1586 * We first construct a matrix of relations based on the partition information,1587 * apply Floyd-Warshall on this matrix of relations and then take the1588 * union of all entries in the matrix as the final result.1589 *1590 * If we are actually computing the power instead of the transitive closure,1591 * i.e., when "project" is not set, then the result should have the1592 * path lengths encoded as the difference between an extra pair of1593 * coordinates.  We therefore apply the nested transitive closures1594 * to relations that include these lengths.  In particular, we replace1595 * the input relation by the cross product with the unit length relation1596 * { [i] -> [i + 1] }.1597 */1598static __isl_give isl_map *floyd_warshall_with_groups(1599	__isl_take isl_space *space, __isl_keep isl_map *map,1600	isl_bool *exact, int project, int *group, int n)1601{1602	int i, j, k;1603	isl_map ***grid = NULL;1604	isl_map *app;1605 1606	if (!map)1607		goto error;1608 1609	if (n == 1) {1610		free(group);1611		return incremental_closure(space, map, exact, project);1612	}1613 1614	grid = isl_calloc_array(map->ctx, isl_map **, n);1615	if (!grid)1616		goto error;1617	for (i = 0; i < n; ++i) {1618		grid[i] = isl_calloc_array(map->ctx, isl_map *, n);1619		if (!grid[i])1620			goto error;1621		for (j = 0; j < n; ++j)1622			grid[i][j] = isl_map_empty(isl_map_get_space(map));1623	}1624 1625	for (k = 0; k < map->n; ++k) {1626		i = group[2 * k];1627		j = group[2 * k + 1];1628		grid[i][j] = isl_map_union(grid[i][j],1629				isl_map_from_basic_map(1630					isl_basic_map_copy(map->p[k])));1631	}1632 1633	if (!project && add_length(map, grid, n) < 0)1634		goto error;1635 1636	floyd_warshall_iterate(grid, n, exact);1637 1638	app = isl_map_empty(isl_map_get_space(grid[0][0]));1639 1640	for (i = 0; i < n; ++i) {1641		for (j = 0; j < n; ++j)1642			app = isl_map_union(app, grid[i][j]);1643		free(grid[i]);1644	}1645	free(grid);1646 1647	free(group);1648	isl_space_free(space);1649 1650	return app;1651error:1652	if (grid)1653		for (i = 0; i < n; ++i) {1654			if (!grid[i])1655				continue;1656			for (j = 0; j < n; ++j)1657				isl_map_free(grid[i][j]);1658			free(grid[i]);1659		}1660	free(grid);1661	free(group);1662	isl_space_free(space);1663	return NULL;1664}1665 1666/* Partition the domains and ranges of the n basic relations in list1667 * into disjoint cells.1668 *1669 * To find the partition, we simply consider all of the domains1670 * and ranges in turn and combine those that overlap.1671 * "set" contains the partition elements and "group" indicates1672 * to which partition element a given domain or range belongs.1673 * The domain of basic map i corresponds to element 2 * i in these arrays,1674 * while the domain corresponds to element 2 * i + 1.1675 * During the construction group[k] is either equal to k,1676 * in which case set[k] contains the union of all the domains and1677 * ranges in the corresponding group, or is equal to some l < k,1678 * with l another domain or range in the same group.1679 */1680static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,1681	isl_set ***set, int *n_group)1682{1683	int i;1684	int *group = NULL;1685	int g;1686 1687	*set = isl_calloc_array(ctx, isl_set *, 2 * n);1688	group = isl_alloc_array(ctx, int, 2 * n);1689 1690	if (!*set || !group)1691		goto error;1692 1693	for (i = 0; i < n; ++i) {1694		isl_set *dom;1695		dom = isl_set_from_basic_set(isl_basic_map_domain(1696				isl_basic_map_copy(list[i])));1697		if (merge(*set, group, dom, 2 * i) < 0)1698			goto error;1699		dom = isl_set_from_basic_set(isl_basic_map_range(1700				isl_basic_map_copy(list[i])));1701		if (merge(*set, group, dom, 2 * i + 1) < 0)1702			goto error;1703	}1704 1705	g = 0;1706	for (i = 0; i < 2 * n; ++i)1707		if (group[i] == i) {1708			if (g != i) {1709				(*set)[g] = (*set)[i];1710				(*set)[i] = NULL;1711			}1712			group[i] = g++;1713		} else1714			group[i] = group[group[i]];1715 1716	*n_group = g;1717 1718	return group;1719error:1720	if (*set) {1721		for (i = 0; i < 2 * n; ++i)1722			isl_set_free((*set)[i]);1723		free(*set);1724		*set = NULL;1725	}1726	free(group);1727	return NULL;1728}1729 1730/* Check if the domains and ranges of the basic maps in "map" can1731 * be partitioned, and if so, apply Floyd-Warshall on the elements1732 * of the partition.  Note that we also apply this algorithm1733 * if we want to compute the power, i.e., when "project" is not set.1734 * However, the results are unlikely to be exact since the recursive1735 * calls inside the Floyd-Warshall algorithm typically result in1736 * non-linear path lengths quite quickly.1737 */1738static __isl_give isl_map *floyd_warshall(__isl_take isl_space *space,1739	__isl_keep isl_map *map, isl_bool *exact, int project)1740{1741	int i;1742	isl_set **set = NULL;1743	int *group = NULL;1744	int n;1745 1746	if (!map)1747		goto error;1748	if (map->n <= 1)1749		return incremental_closure(space, map, exact, project);1750 1751	group = setup_groups(map->ctx, map->p, map->n, &set, &n);1752	if (!group)1753		goto error;1754 1755	for (i = 0; i < 2 * map->n; ++i)1756		isl_set_free(set[i]);1757 1758	free(set);1759 1760	return floyd_warshall_with_groups(space, map, exact, project, group, n);1761error:1762	isl_space_free(space);1763	return NULL;1764}1765 1766/* Structure for representing the nodes of the graph of which1767 * strongly connected components are being computed.1768 *1769 * list contains the actual nodes1770 * check_closed is set if we may have used the fact that1771 * a pair of basic maps can be interchanged1772 */1773struct isl_tc_follows_data {1774	isl_basic_map **list;1775	int check_closed;1776};1777 1778/* Check whether in the computation of the transitive closure1779 * "list[i]" (R_1) should follow (or be part of the same component as)1780 * "list[j]" (R_2).1781 *1782 * That is check whether1783 *1784 *	R_1 \circ R_21785 *1786 * is a subset of1787 *1788 *	R_2 \circ R_11789 *1790 * If so, then there is no reason for R_1 to immediately follow R_21791 * in any path.1792 *1793 * *check_closed is set if the subset relation holds while1794 * R_1 \circ R_2 is not empty.1795 */1796static isl_bool basic_map_follows(int i, int j, void *user)1797{1798	struct isl_tc_follows_data *data = user;1799	struct isl_map *map12 = NULL;1800	struct isl_map *map21 = NULL;1801	isl_bool applies, subset;1802 1803	applies = isl_basic_map_applies_range(data->list[j], data->list[i]);1804	if (applies < 0)1805		return isl_bool_error;1806	if (!applies)1807		return isl_bool_false;1808 1809	map21 = isl_map_from_basic_map(1810			isl_basic_map_apply_range(1811				isl_basic_map_copy(data->list[j]),1812				isl_basic_map_copy(data->list[i])));1813	subset = isl_map_is_empty(map21);1814	if (subset < 0)1815		goto error;1816	if (subset) {1817		isl_map_free(map21);1818		return isl_bool_false;1819	}1820 1821	if (!isl_basic_map_is_transformation(data->list[i]) ||1822	    !isl_basic_map_is_transformation(data->list[j])) {1823		isl_map_free(map21);1824		return isl_bool_true;1825	}1826 1827	map12 = isl_map_from_basic_map(1828			isl_basic_map_apply_range(1829				isl_basic_map_copy(data->list[i]),1830				isl_basic_map_copy(data->list[j])));1831 1832	subset = isl_map_is_subset(map21, map12);1833 1834	isl_map_free(map12);1835	isl_map_free(map21);1836 1837	if (subset)1838		data->check_closed = 1;1839 1840	return isl_bool_not(subset);1841error:1842	isl_map_free(map21);1843	return isl_bool_error;1844}1845 1846/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D1847 * and a dimension specification (Z^{n+1} -> Z^{n+1}),1848 * construct a map that is an overapproximation of the map1849 * that takes an element from the dom R \times Z to an1850 * element from ran R \times Z, such that the first n coordinates of the1851 * difference between them is a sum of differences between images1852 * and pre-images in one of the R_i and such that the last coordinate1853 * is equal to the number of steps taken.1854 * If "project" is set, then these final coordinates are not included,1855 * i.e., a relation of type Z^n -> Z^n is returned.1856 * That is, let1857 *1858 *	\Delta_i = { y - x | (x, y) in R_i }1859 *1860 * then the constructed map is an overapproximation of1861 *1862 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :1863 *				d = (\sum_i k_i \delta_i, \sum_i k_i) and1864 *				x in dom R and x + d in ran R }1865 *1866 * or1867 *1868 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :1869 *				d = (\sum_i k_i \delta_i) and1870 *				x in dom R and x + d in ran R }1871 *1872 * if "project" is set.1873 *1874 * We first split the map into strongly connected components, perform1875 * the above on each component and then join the results in the correct1876 * order, at each join also taking in the union of both arguments1877 * to allow for paths that do not go through one of the two arguments.1878 */1879static __isl_give isl_map *construct_power_components(1880	__isl_take isl_space *space, __isl_keep isl_map *map, isl_bool *exact,1881	int project)1882{1883	int i, n, c;1884	struct isl_map *path = NULL;1885	struct isl_tc_follows_data data;1886	struct isl_tarjan_graph *g = NULL;1887	isl_bool *orig_exact;1888	isl_bool local_exact;1889 1890	if (!map)1891		goto error;1892	if (map->n <= 1)1893		return floyd_warshall(space, map, exact, project);1894 1895	data.list = map->p;1896	data.check_closed = 0;1897	g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data);1898	if (!g)1899		goto error;1900 1901	orig_exact = exact;1902	if (data.check_closed && !exact)1903		exact = &local_exact;1904 1905	c = 0;1906	i = 0;1907	n = map->n;1908	if (project)1909		path = isl_map_empty(isl_map_get_space(map));1910	else1911		path = isl_map_empty(isl_space_copy(space));1912	path = anonymize(path);1913	while (n) {1914		struct isl_map *comp;1915		isl_map *path_comp, *path_comb;1916		comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);1917		while (g->order[i] != -1) {1918			comp = isl_map_add_basic_map(comp,1919				    isl_basic_map_copy(map->p[g->order[i]]));1920			--n;1921			++i;1922		}1923		path_comp = floyd_warshall(isl_space_copy(space),1924						comp, exact, project);1925		path_comp = anonymize(path_comp);1926		path_comb = isl_map_apply_range(isl_map_copy(path),1927						isl_map_copy(path_comp));1928		path = isl_map_union(path, path_comp);1929		path = isl_map_union(path, path_comb);1930		isl_map_free(comp);1931		++i;1932		++c;1933	}1934 1935	if (c > 1 && data.check_closed && !*exact) {1936		isl_bool closed;1937 1938		closed = isl_map_is_transitively_closed(path);1939		if (closed < 0)1940			goto error;1941		if (!closed) {1942			isl_tarjan_graph_free(g);1943			isl_map_free(path);1944			return floyd_warshall(space, map, orig_exact, project);1945		}1946	}1947 1948	isl_tarjan_graph_free(g);1949	isl_space_free(space);1950 1951	return path;1952error:1953	isl_tarjan_graph_free(g);1954	isl_space_free(space);1955	isl_map_free(path);1956	return NULL;1957}1958 1959/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,1960 * construct a map that is an overapproximation of the map1961 * that takes an element from the space D to another1962 * element from the same space, such that the difference between1963 * them is a strictly positive sum of differences between images1964 * and pre-images in one of the R_i.1965 * The number of differences in the sum is equated to parameter "param".1966 * That is, let1967 *1968 *	\Delta_i = { y - x | (x, y) in R_i }1969 *1970 * then the constructed map is an overapproximation of1971 *1972 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :1973 *				d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }1974 * or1975 *1976 *	{ (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :1977 *				d = \sum_i k_i \delta_i and \sum_i k_i > 0 }1978 *1979 * if "project" is set.1980 *1981 * If "project" is not set, then1982 * we construct an extended mapping with an extra coordinate1983 * that indicates the number of steps taken.  In particular,1984 * the difference in the last coordinate is equal to the number1985 * of steps taken to move from a domain element to the corresponding1986 * image element(s).1987 */1988static __isl_give isl_map *construct_power(__isl_keep isl_map *map,1989	isl_bool *exact, int project)1990{1991	struct isl_map *app = NULL;1992	isl_space *space = NULL;1993 1994	if (!map)1995		return NULL;1996 1997	space = isl_map_get_space(map);1998 1999	space = isl_space_add_dims(space, isl_dim_in, 1);2000	space = isl_space_add_dims(space, isl_dim_out, 1);2001 2002	app = construct_power_components(isl_space_copy(space), map,2003					exact, project);2004 2005	isl_space_free(space);2006 2007	return app;2008}2009 2010/* Compute the positive powers of "map", or an overapproximation.2011 * If the result is exact, then *exact is set to 1.2012 *2013 * If project is set, then we are actually interested in the transitive2014 * closure, so we can use a more relaxed exactness check.2015 * The lengths of the paths are also projected out instead of being2016 * encoded as the difference between an extra pair of final coordinates.2017 */2018static __isl_give isl_map *map_power(__isl_take isl_map *map,2019	isl_bool *exact, int project)2020{2021	struct isl_map *app = NULL;2022 2023	if (exact)2024		*exact = isl_bool_true;2025 2026	if (isl_map_check_transformation(map) < 0)2027		return isl_map_free(map);2028 2029	app = construct_power(map, exact, project);2030 2031	isl_map_free(map);2032	return app;2033}2034 2035/* Compute the positive powers of "map", or an overapproximation.2036 * The result maps the exponent to a nested copy of the corresponding power.2037 * If the result is exact, then *exact is set to 1.2038 * map_power constructs an extended relation with the path lengths2039 * encoded as the difference between the final coordinates.2040 * In the final step, this difference is equated to an extra parameter2041 * and made positive.  The extra coordinates are subsequently projected out2042 * and the parameter is turned into the domain of the result.2043 */2044__isl_give isl_map *isl_map_power(__isl_take isl_map *map, isl_bool *exact)2045{2046	isl_space *target_space;2047	isl_space *space;2048	isl_map *diff;2049	isl_size d;2050	isl_size param;2051 2052	d = isl_map_dim(map, isl_dim_in);2053	param = isl_map_dim(map, isl_dim_param);2054	if (d < 0 || param < 0)2055		return isl_map_free(map);2056 2057	map = isl_map_compute_divs(map);2058	map = isl_map_coalesce(map);2059 2060	if (isl_map_plain_is_empty(map)) {2061		map = isl_map_from_range(isl_map_wrap(map));2062		map = isl_map_add_dims(map, isl_dim_in, 1);2063		map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");2064		return map;2065	}2066 2067	target_space = isl_map_get_space(map);2068	target_space = isl_space_from_range(isl_space_wrap(target_space));2069	target_space = isl_space_add_dims(target_space, isl_dim_in, 1);2070	target_space = isl_space_set_dim_name(target_space, isl_dim_in, 0, "k");2071 2072	map = map_power(map, exact, 0);2073 2074	map = isl_map_add_dims(map, isl_dim_param, 1);2075	space = isl_map_get_space(map);2076	diff = equate_parameter_to_length(space, param);2077	map = isl_map_intersect(map, diff);2078	map = isl_map_project_out(map, isl_dim_in, d, 1);2079	map = isl_map_project_out(map, isl_dim_out, d, 1);2080	map = isl_map_from_range(isl_map_wrap(map));2081	map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);2082 2083	map = isl_map_reset_space(map, target_space);2084 2085	return map;2086}2087 2088/* Compute a relation that maps each element in the range of the input2089 * relation to the lengths of all paths composed of edges in the input2090 * relation that end up in the given range element.2091 * The result may be an overapproximation, in which case *exact is set to 0.2092 * The resulting relation is very similar to the power relation.2093 * The difference are that the domain has been projected out, the2094 * range has become the domain and the exponent is the range instead2095 * of a parameter.2096 */2097__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,2098	isl_bool *exact)2099{2100	isl_space *space;2101	isl_map *diff;2102	isl_size d;2103	isl_size param;2104 2105	d = isl_map_dim(map, isl_dim_in);2106	param = isl_map_dim(map, isl_dim_param);2107	if (d < 0 || param < 0)2108		return isl_map_free(map);2109 2110	map = isl_map_compute_divs(map);2111	map = isl_map_coalesce(map);2112 2113	if (isl_map_plain_is_empty(map)) {2114		if (exact)2115			*exact = isl_bool_true;2116		map = isl_map_project_out(map, isl_dim_out, 0, d);2117		map = isl_map_add_dims(map, isl_dim_out, 1);2118		return map;2119	}2120 2121	map = map_power(map, exact, 0);2122 2123	map = isl_map_add_dims(map, isl_dim_param, 1);2124	space = isl_map_get_space(map);2125	diff = equate_parameter_to_length(space, param);2126	map = isl_map_intersect(map, diff);2127	map = isl_map_project_out(map, isl_dim_in, 0, d + 1);2128	map = isl_map_project_out(map, isl_dim_out, d, 1);2129	map = isl_map_reverse(map);2130	map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);2131 2132	return map;2133}2134 2135/* Given a map, compute the smallest superset of this map that is of the form2136 *2137 *	{ i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }2138 *2139 * (where p ranges over the (non-parametric) dimensions),2140 * compute the transitive closure of this map, i.e.,2141 *2142 *	{ i -> j : exists k > 0:2143 *		k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }2144 *2145 * and intersect domain and range of this transitive closure with2146 * the given domain and range.2147 *2148 * If with_id is set, then try to include as much of the identity mapping2149 * as possible, by computing2150 *2151 *	{ i -> j : exists k >= 0:2152 *		k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }2153 *2154 * instead (i.e., allow k = 0).2155 *2156 * In practice, we compute the difference set2157 *2158 *	delta  = { j - i | i -> j in map },2159 *2160 * look for stride constraint on the individual dimensions and compute2161 * (constant) lower and upper bounds for each individual dimension,2162 * adding a constraint for each bound not equal to infinity.2163 */2164static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,2165	__isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)2166{2167	int i;2168	int k;2169	unsigned d;2170	unsigned nparam;2171	unsigned total;2172	isl_space *space;2173	isl_set *delta;2174	isl_map *app = NULL;2175	isl_basic_set *aff = NULL;2176	isl_basic_map *bmap = NULL;2177	isl_vec *obj = NULL;2178	isl_int opt;2179 2180	isl_int_init(opt);2181 2182	delta = isl_map_deltas(isl_map_copy(map));2183 2184	aff = isl_set_affine_hull(isl_set_copy(delta));2185	if (!aff)2186		goto error;2187	space = isl_map_get_space(map);2188	d = isl_space_dim(space, isl_dim_in);2189	nparam = isl_space_dim(space, isl_dim_param);2190	total = isl_space_dim(space, isl_dim_all);2191	bmap = isl_basic_map_alloc_space(space,2192					aff->n_div + 1, aff->n_div, 2 * d + 1);2193	for (i = 0; i < aff->n_div + 1; ++i) {2194		k = isl_basic_map_alloc_div(bmap);2195		if (k < 0)2196			goto error;2197		isl_int_set_si(bmap->div[k][0], 0);2198	}2199	for (i = 0; i < aff->n_eq; ++i) {2200		if (!isl_basic_set_eq_is_stride(aff, i))2201			continue;2202		k = isl_basic_map_alloc_equality(bmap);2203		if (k < 0)2204			goto error;2205		isl_seq_clr(bmap->eq[k], 1 + nparam);2206		isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,2207				aff->eq[i] + 1 + nparam, d);2208		isl_seq_neg(bmap->eq[k] + 1 + nparam,2209				aff->eq[i] + 1 + nparam, d);2210		isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,2211				aff->eq[i] + 1 + nparam + d, aff->n_div);2212		isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);2213	}2214	obj = isl_vec_alloc(map->ctx, 1 + nparam + d);2215	if (!obj)2216		goto error;2217	isl_seq_clr(obj->el, 1 + nparam + d);2218	for (i = 0; i < d; ++ i) {2219		enum isl_lp_result res;2220 2221		isl_int_set_si(obj->el[1 + nparam + i], 1);2222 2223		res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,2224					NULL, NULL);2225		if (res == isl_lp_error)2226			goto error;2227		if (res == isl_lp_ok) {2228			k = isl_basic_map_alloc_inequality(bmap);2229			if (k < 0)2230				goto error;2231			isl_seq_clr(bmap->ineq[k],2232					1 + nparam + 2 * d + bmap->n_div);2233			isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);2234			isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);2235			isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);2236		}2237 2238		res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,2239					NULL, NULL);2240		if (res == isl_lp_error)2241			goto error;2242		if (res == isl_lp_ok) {2243			k = isl_basic_map_alloc_inequality(bmap);2244			if (k < 0)2245				goto error;2246			isl_seq_clr(bmap->ineq[k],2247					1 + nparam + 2 * d + bmap->n_div);2248			isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);2249			isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);2250			isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);2251		}2252 2253		isl_int_set_si(obj->el[1 + nparam + i], 0);2254	}2255	k = isl_basic_map_alloc_inequality(bmap);2256	if (k < 0)2257		goto error;2258	isl_seq_clr(bmap->ineq[k],2259			1 + nparam + 2 * d + bmap->n_div);2260	if (!with_id)2261		isl_int_set_si(bmap->ineq[k][0], -1);2262	isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);2263 2264	app = isl_map_from_domain_and_range(dom, ran);2265 2266	isl_vec_free(obj);2267	isl_basic_set_free(aff);2268	isl_map_free(map);2269	bmap = isl_basic_map_finalize(bmap);2270	isl_set_free(delta);2271	isl_int_clear(opt);2272 2273	map = isl_map_from_basic_map(bmap);2274	map = isl_map_intersect(map, app);2275 2276	return map;2277error:2278	isl_vec_free(obj);2279	isl_basic_map_free(bmap);2280	isl_basic_set_free(aff);2281	isl_set_free(dom);2282	isl_set_free(ran);2283	isl_map_free(map);2284	isl_set_free(delta);2285	isl_int_clear(opt);2286	return NULL;2287}2288 2289/* Given a map, compute the smallest superset of this map that is of the form2290 *2291 *	{ i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }2292 *2293 * (where p ranges over the (non-parametric) dimensions),2294 * compute the transitive closure of this map, i.e.,2295 *2296 *	{ i -> j : exists k > 0:2297 *		k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }2298 *2299 * and intersect domain and range of this transitive closure with2300 * domain and range of the original map.2301 */2302static __isl_give isl_map *box_closure(__isl_take isl_map *map)2303{2304	isl_set *domain;2305	isl_set *range;2306 2307	domain = isl_map_domain(isl_map_copy(map));2308	domain = isl_set_coalesce(domain);2309	range = isl_map_range(isl_map_copy(map));2310	range = isl_set_coalesce(range);2311 2312	return box_closure_on_domain(map, domain, range, 0);2313}2314 2315/* Given a map, compute the smallest superset of this map that is of the form2316 *2317 *	{ i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }2318 *2319 * (where p ranges over the (non-parametric) dimensions),2320 * compute the transitive and partially reflexive closure of this map, i.e.,2321 *2322 *	{ i -> j : exists k >= 0:2323 *		k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }2324 *2325 * and intersect domain and range of this transitive closure with2326 * the given domain.2327 */2328static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,2329	__isl_take isl_set *dom)2330{2331	return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);2332}2333 2334/* Check whether app is the transitive closure of map.2335 * In particular, check that app is acyclic and, if so,2336 * check that2337 *2338 *	app \subset (map \cup (map \circ app))2339 */2340static isl_bool check_exactness_omega(__isl_keep isl_map *map,2341	__isl_keep isl_map *app)2342{2343	isl_set *delta;2344	int i;2345	isl_bool is_empty, is_exact;2346	isl_size d;2347	isl_map *test;2348 2349	delta = isl_map_deltas(isl_map_copy(app));2350	d = isl_set_dim(delta, isl_dim_set);2351	if (d < 0)2352		delta = isl_set_free(delta);2353	for (i = 0; i < d; ++i)2354		delta = isl_set_fix_si(delta, isl_dim_set, i, 0);2355	is_empty = isl_set_is_empty(delta);2356	isl_set_free(delta);2357	if (is_empty < 0 || !is_empty)2358		return is_empty;2359 2360	test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));2361	test = isl_map_union(test, isl_map_copy(map));2362	is_exact = isl_map_is_subset(app, test);2363	isl_map_free(test);2364 2365	return is_exact;2366}2367 2368/* Check if basic map M_i can be combined with all the other2369 * basic maps such that2370 *2371 *	(\cup_j M_j)^+2372 *2373 * can be computed as2374 *2375 *	M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+2376 *2377 * In particular, check if we can compute a compact representation2378 * of2379 *2380 *		M_i^* \circ M_j \circ M_i^*2381 *2382 * for each j != i.2383 * Let M_i^? be an extension of M_i^+ that allows paths2384 * of length zero, i.e., the result of box_closure(., 1).2385 * The criterion, as proposed by Kelly et al., is that2386 * id = M_i^? - M_i^+ can be represented as a basic map2387 * and that2388 *2389 *	id \circ M_j \circ id = M_j2390 *2391 * for each j != i.2392 *2393 * If this function returns 1, then tc and qc are set to2394 * M_i^+ and M_i^?, respectively.2395 */2396static int can_be_split_off(__isl_keep isl_map *map, int i,2397	__isl_give isl_map **tc, __isl_give isl_map **qc)2398{2399	isl_map *map_i, *id = NULL;2400	int j = -1;2401	isl_set *C;2402 2403	*tc = NULL;2404	*qc = NULL;2405 2406	C = isl_set_union(isl_map_domain(isl_map_copy(map)),2407			  isl_map_range(isl_map_copy(map)));2408	C = isl_set_from_basic_set(isl_set_simple_hull(C));2409	if (!C)2410		goto error;2411 2412	map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));2413	*tc = box_closure(isl_map_copy(map_i));2414	*qc = box_closure_with_identity(map_i, C);2415	id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));2416 2417	if (!id || !*qc)2418		goto error;2419	if (id->n != 1 || (*qc)->n != 1)2420		goto done;2421 2422	for (j = 0; j < map->n; ++j) {2423		isl_map *map_j, *test;2424		int is_ok;2425 2426		if (i == j)2427			continue;2428		map_j = isl_map_from_basic_map(2429					isl_basic_map_copy(map->p[j]));2430		test = isl_map_apply_range(isl_map_copy(id),2431						isl_map_copy(map_j));2432		test = isl_map_apply_range(test, isl_map_copy(id));2433		is_ok = isl_map_is_equal(test, map_j);2434		isl_map_free(map_j);2435		isl_map_free(test);2436		if (is_ok < 0)2437			goto error;2438		if (!is_ok)2439			break;2440	}2441 2442done:2443	isl_map_free(id);2444	if (j == map->n)2445		return 1;2446 2447	isl_map_free(*qc);2448	isl_map_free(*tc);2449	*qc = NULL;2450	*tc = NULL;2451 2452	return 0;2453error:2454	isl_map_free(id);2455	isl_map_free(*qc);2456	isl_map_free(*tc);2457	*qc = NULL;2458	*tc = NULL;2459	return -1;2460}2461 2462static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,2463	isl_bool *exact)2464{2465	isl_map *app;2466 2467	app = box_closure(isl_map_copy(map));2468	if (exact) {2469		isl_bool is_exact = check_exactness_omega(map, app);2470 2471		if (is_exact < 0)2472			app = isl_map_free(app);2473		else2474			*exact = is_exact;2475	}2476 2477	isl_map_free(map);2478	return app;2479}2480 2481/* Compute an overapproximation of the transitive closure of "map"2482 * using a variation of the algorithm from2483 * "Transitive Closure of Infinite Graphs and its Applications"2484 * by Kelly et al.2485 *2486 * We first check whether we can can split of any basic map M_i and2487 * compute2488 *2489 *	(\cup_j M_j)^+2490 *2491 * as2492 *2493 *	M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+2494 *2495 * using a recursive call on the remaining map.2496 *2497 * If not, we simply call box_closure on the whole map.2498 */2499static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,2500	isl_bool *exact)2501{2502	int i, j;2503	isl_bool exact_i;2504	isl_map *app;2505 2506	if (!map)2507		return NULL;2508	if (map->n == 1)2509		return box_closure_with_check(map, exact);2510 2511	for (i = 0; i < map->n; ++i) {2512		int ok;2513		isl_map *qc, *tc;2514		ok = can_be_split_off(map, i, &tc, &qc);2515		if (ok < 0)2516			goto error;2517		if (!ok)2518			continue;2519 2520		app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);2521 2522		for (j = 0; j < map->n; ++j) {2523			if (j == i)2524				continue;2525			app = isl_map_add_basic_map(app,2526						isl_basic_map_copy(map->p[j]));2527		}2528 2529		app = isl_map_apply_range(isl_map_copy(qc), app);2530		app = isl_map_apply_range(app, qc);2531 2532		app = isl_map_union(tc, transitive_closure_omega(app, NULL));2533		exact_i = check_exactness_omega(map, app);2534		if (exact_i == isl_bool_true) {2535			if (exact)2536				*exact = exact_i;2537			isl_map_free(map);2538			return app;2539		}2540		isl_map_free(app);2541		if (exact_i < 0)2542			goto error;2543	}2544 2545	return box_closure_with_check(map, exact);2546error:2547	isl_map_free(map);2548	return NULL;2549}2550 2551/* Compute the transitive closure  of "map", or an overapproximation.2552 * If the result is exact, then *exact is set to 1.2553 * Simply use map_power to compute the powers of map, but tell2554 * it to project out the lengths of the paths instead of equating2555 * the length to a parameter.2556 */2557__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,2558	isl_bool *exact)2559{2560	isl_space *target_dim;2561	isl_bool closed;2562 2563	if (!map)2564		goto error;2565 2566	if (map->ctx->opt->closure == ISL_CLOSURE_BOX)2567		return transitive_closure_omega(map, exact);2568 2569	map = isl_map_compute_divs(map);2570	map = isl_map_coalesce(map);2571	closed = isl_map_is_transitively_closed(map);2572	if (closed < 0)2573		goto error;2574	if (closed) {2575		if (exact)2576			*exact = isl_bool_true;2577		return map;2578	}2579 2580	target_dim = isl_map_get_space(map);2581	map = map_power(map, exact, 1);2582	map = isl_map_reset_space(map, target_dim);2583 2584	return map;2585error:2586	isl_map_free(map);2587	return NULL;2588}2589 2590static isl_stat inc_count(__isl_take isl_map *map, void *user)2591{2592	int *n = user;2593 2594	*n += map->n;2595 2596	isl_map_free(map);2597 2598	return isl_stat_ok;2599}2600 2601static isl_stat collect_basic_map(__isl_take isl_map *map, void *user)2602{2603	int i;2604	isl_basic_map ***next = user;2605 2606	for (i = 0; i < map->n; ++i) {2607		**next = isl_basic_map_copy(map->p[i]);2608		if (!**next)2609			goto error;2610		(*next)++;2611	}2612 2613	isl_map_free(map);2614	return isl_stat_ok;2615error:2616	isl_map_free(map);2617	return isl_stat_error;2618}2619 2620/* Perform Floyd-Warshall on the given list of basic relations.2621 * The basic relations may live in different dimensions,2622 * but basic relations that get assigned to the diagonal of the2623 * grid have domains and ranges of the same dimension and so2624 * the standard algorithm can be used because the nested transitive2625 * closures are only applied to diagonal elements and because all2626 * compositions are performed on relations with compatible domains and ranges.2627 */2628static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,2629	__isl_keep isl_basic_map **list, int n, isl_bool *exact)2630{2631	int i, j, k;2632	int n_group;2633	int *group = NULL;2634	isl_set **set = NULL;2635	isl_map ***grid = NULL;2636	isl_union_map *app;2637 2638	group = setup_groups(ctx, list, n, &set, &n_group);2639	if (!group)2640		goto error;2641 2642	grid = isl_calloc_array(ctx, isl_map **, n_group);2643	if (!grid)2644		goto error;2645	for (i = 0; i < n_group; ++i) {2646		grid[i] = isl_calloc_array(ctx, isl_map *, n_group);2647		if (!grid[i])2648			goto error;2649		for (j = 0; j < n_group; ++j) {2650			isl_space *space1, *space2, *space;2651			space1 = isl_space_reverse(isl_set_get_space(set[i]));2652			space2 = isl_set_get_space(set[j]);2653			space = isl_space_join(space1, space2);2654			grid[i][j] = isl_map_empty(space);2655		}2656	}2657 2658	for (k = 0; k < n; ++k) {2659		i = group[2 * k];2660		j = group[2 * k + 1];2661		grid[i][j] = isl_map_union(grid[i][j],2662				isl_map_from_basic_map(2663					isl_basic_map_copy(list[k])));2664	}2665	2666	floyd_warshall_iterate(grid, n_group, exact);2667 2668	app = isl_union_map_empty(isl_map_get_space(grid[0][0]));2669 2670	for (i = 0; i < n_group; ++i) {2671		for (j = 0; j < n_group; ++j)2672			app = isl_union_map_add_map(app, grid[i][j]);2673		free(grid[i]);2674	}2675	free(grid);2676 2677	for (i = 0; i < 2 * n; ++i)2678		isl_set_free(set[i]);2679	free(set);2680 2681	free(group);2682	return app;2683error:2684	if (grid)2685		for (i = 0; i < n_group; ++i) {2686			if (!grid[i])2687				continue;2688			for (j = 0; j < n_group; ++j)2689				isl_map_free(grid[i][j]);2690			free(grid[i]);2691		}2692	free(grid);2693	if (set) {2694		for (i = 0; i < 2 * n; ++i)2695			isl_set_free(set[i]);2696		free(set);2697	}2698	free(group);2699	return NULL;2700}2701 2702/* Perform Floyd-Warshall on the given union relation.2703 * The implementation is very similar to that for non-unions.2704 * The main difference is that it is applied unconditionally.2705 * We first extract a list of basic maps from the union map2706 * and then perform the algorithm on this list.2707 */2708static __isl_give isl_union_map *union_floyd_warshall(2709	__isl_take isl_union_map *umap, isl_bool *exact)2710{2711	int i, n;2712	isl_ctx *ctx;2713	isl_basic_map **list = NULL;2714	isl_basic_map **next;2715	isl_union_map *res;2716 2717	n = 0;2718	if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)2719		goto error;2720 2721	ctx = isl_union_map_get_ctx(umap);2722	list = isl_calloc_array(ctx, isl_basic_map *, n);2723	if (!list)2724		goto error;2725 2726	next = list;2727	if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)2728		goto error;2729 2730	res = union_floyd_warshall_on_list(ctx, list, n, exact);2731 2732	if (list) {2733		for (i = 0; i < n; ++i)2734			isl_basic_map_free(list[i]);2735		free(list);2736	}2737 2738	isl_union_map_free(umap);2739	return res;2740error:2741	if (list) {2742		for (i = 0; i < n; ++i)2743			isl_basic_map_free(list[i]);2744		free(list);2745	}2746	isl_union_map_free(umap);2747	return NULL;2748}2749 2750/* Decompose the give union relation into strongly connected components.2751 * The implementation is essentially the same as that of2752 * construct_power_components with the major difference that all2753 * operations are performed on union maps.2754 */2755static __isl_give isl_union_map *union_components(2756	__isl_take isl_union_map *umap, isl_bool *exact)2757{2758	int i;2759	int n;2760	isl_ctx *ctx;2761	isl_basic_map **list = NULL;2762	isl_basic_map **next;2763	isl_union_map *path = NULL;2764	struct isl_tc_follows_data data;2765	struct isl_tarjan_graph *g = NULL;2766	int c, l;2767	int recheck = 0;2768 2769	n = 0;2770	if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)2771		goto error;2772 2773	if (n == 0)2774		return umap;2775	if (n <= 1)2776		return union_floyd_warshall(umap, exact);2777 2778	ctx = isl_union_map_get_ctx(umap);2779	list = isl_calloc_array(ctx, isl_basic_map *, n);2780	if (!list)2781		goto error;2782 2783	next = list;2784	if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)2785		goto error;2786 2787	data.list = list;2788	data.check_closed = 0;2789	g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);2790	if (!g)2791		goto error;2792 2793	c = 0;2794	i = 0;2795	l = n;2796	path = isl_union_map_empty(isl_union_map_get_space(umap));2797	while (l) {2798		isl_union_map *comp;2799		isl_union_map *path_comp, *path_comb;2800		comp = isl_union_map_empty(isl_union_map_get_space(umap));2801		while (g->order[i] != -1) {2802			comp = isl_union_map_add_map(comp,2803				    isl_map_from_basic_map(2804					isl_basic_map_copy(list[g->order[i]])));2805			--l;2806			++i;2807		}2808		path_comp = union_floyd_warshall(comp, exact);2809		path_comb = isl_union_map_apply_range(isl_union_map_copy(path),2810						isl_union_map_copy(path_comp));2811		path = isl_union_map_union(path, path_comp);2812		path = isl_union_map_union(path, path_comb);2813		++i;2814		++c;2815	}2816 2817	if (c > 1 && data.check_closed && !*exact) {2818		isl_bool closed;2819 2820		closed = isl_union_map_is_transitively_closed(path);2821		if (closed < 0)2822			goto error;2823		recheck = !closed;2824	}2825 2826	isl_tarjan_graph_free(g);2827 2828	for (i = 0; i < n; ++i)2829		isl_basic_map_free(list[i]);2830	free(list);2831 2832	if (recheck) {2833		isl_union_map_free(path);2834		return union_floyd_warshall(umap, exact);2835	}2836 2837	isl_union_map_free(umap);2838 2839	return path;2840error:2841	isl_tarjan_graph_free(g);2842	if (list) {2843		for (i = 0; i < n; ++i)2844			isl_basic_map_free(list[i]);2845		free(list);2846	}2847	isl_union_map_free(umap);2848	isl_union_map_free(path);2849	return NULL;2850}2851 2852/* Compute the transitive closure  of "umap", or an overapproximation.2853 * If the result is exact, then *exact is set to 1.2854 */2855__isl_give isl_union_map *isl_union_map_transitive_closure(2856	__isl_take isl_union_map *umap, isl_bool *exact)2857{2858	isl_bool closed;2859 2860	if (!umap)2861		return NULL;2862 2863	if (exact)2864		*exact = isl_bool_true;2865 2866	umap = isl_union_map_compute_divs(umap);2867	umap = isl_union_map_coalesce(umap);2868	closed = isl_union_map_is_transitively_closed(umap);2869	if (closed < 0)2870		goto error;2871	if (closed)2872		return umap;2873	umap = union_components(umap, exact);2874	return umap;2875error:2876	isl_union_map_free(umap);2877	return NULL;2878}2879 2880struct isl_union_power {2881	isl_union_map *pow;2882	isl_bool *exact;2883};2884 2885static isl_stat power(__isl_take isl_map *map, void *user)2886{2887	struct isl_union_power *up = user;2888 2889	map = isl_map_power(map, up->exact);2890	up->pow = isl_union_map_from_map(map);2891 2892	return isl_stat_error;2893}2894 2895/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "space".2896 */2897static __isl_give isl_union_map *deltas_map(__isl_take isl_space *space)2898{2899	isl_basic_map *bmap;2900 2901	space = isl_space_add_dims(space, isl_dim_in, 1);2902	space = isl_space_add_dims(space, isl_dim_out, 1);2903	bmap = isl_basic_map_universe(space);2904	bmap = isl_basic_map_deltas_map(bmap);2905 2906	return isl_union_map_from_map(isl_map_from_basic_map(bmap));2907}2908 2909/* Compute the positive powers of "map", or an overapproximation.2910 * The result maps the exponent to a nested copy of the corresponding power.2911 * If the result is exact, then *exact is set to 1.2912 */2913__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,2914	isl_bool *exact)2915{2916	isl_size n;2917	isl_union_map *inc;2918	isl_union_map *dm;2919 2920	n = isl_union_map_n_map(umap);2921	if (n < 0)2922		return isl_union_map_free(umap);2923	if (n == 0)2924		return umap;2925	if (n == 1) {2926		struct isl_union_power up = { NULL, exact };2927		isl_union_map_foreach_map(umap, &power, &up);2928		isl_union_map_free(umap);2929		return up.pow;2930	}2931	inc = isl_union_map_from_map(increment(isl_union_map_get_space(umap)));2932	umap = isl_union_map_product(inc, umap);2933	umap = isl_union_map_transitive_closure(umap, exact);2934	umap = isl_union_map_zip(umap);2935	dm = deltas_map(isl_union_map_get_space(umap));2936	umap = isl_union_map_apply_domain(umap, dm);2937	2938	return umap;2939}2940 2941#undef TYPE2942#define TYPE isl_map2943#include "isl_power_templ.c"2944 2945#undef TYPE2946#define TYPE isl_union_map2947#include "isl_power_templ.c"2948