1	/*
2	* Copyright 2010 INRIA Saclay
3	*
4	* Use of this software is governed by the MIT license
5	*
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
23	isl_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(map1: isl_map_copy(map), map2: isl_map_copy(map));
29	closed = isl_map_is_subset(map1: map2, map2: map);
30	isl_map_free(map: map2);
31
32	return closed;
33	}
34
35	isl_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(umap1: isl_union_map_copy(umap),
41	umap2: isl_union_map_copy(umap));
42	closed = isl_union_map_is_subset(umap1: umap2, umap2: umap);
43	isl_union_map_free(umap: umap2);
44
45	return closed;
46	}
47
48	/* Given a map that represents a path with the length of the path
49	* encoded as the difference between the last output coordindate
50	* and the last input coordinate, set this length to either
51	* exactly "length" (if "exactly" is set) or at least "length"
52	* (if "exactly" is not set).
53	*/
54	static __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, type: isl_dim_in);
70	nparam = isl_space_dim(space, type: isl_dim_param);
71	total = isl_space_dim(space, type: 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, extra: 0, n_eq: 1, n_ineq: 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(p: c, len: 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(map1: map, map2: isl_map_from_basic_map(bmap));
93
94	return map;
95	error:
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 exactly
102	* the power of "map". Let R be "map" and A_k the overapproximation.
103	* The approximation is exact if
104	*
105	* A_1 = R
106	* A_k = A_{k-1} \circ R k >= 2
107	*
108	* Since A_k is known to be an overapproximation, we only need to check
109	*
110	* A_1 \subset R
111	* A_k \subset A_{k-1} \circ R k >= 2
112	*
113	* In practice, "app" has an extra input and output coordinate
114	* to encode the length of the path. So, we first need to add
115	* this coordinate to "map" and set the length of the path to
116	* one.
117	*/
118	static 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, type: isl_dim_in, n: 1);
126	map = isl_map_add_dims(map, type: isl_dim_out, n: 1);
127	map = set_path_length(map, exactly: 1, length: 1);
128
129	app_1 = set_path_length(map: isl_map_copy(map: app), exactly: 1, length: 1);
130
131	exact = isl_map_is_subset(map1: app_1, map2: map);
132	isl_map_free(map: app_1);
133
134	if (!exact || exact < 0) {
135	isl_map_free(map: app);
136	isl_map_free(map);
137	return exact;
138	}
139
140	app_1 = set_path_length(map: isl_map_copy(map: app), exactly: 0, length: 1);
141	app_2 = set_path_length(map: app, exactly: 0, length: 2);
142	app_1 = isl_map_apply_range(map1: map, map2: app_1);
143
144	exact = isl_map_is_subset(map1: app_2, map2: app_1);
145
146	isl_map_free(map: app_1);
147	isl_map_free(map: app_2);
148
149	return exact;
150	}
151
152	/* Check whether the overapproximation of the power of "map" is exactly
153	* 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 check
158	*
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 check
164	*
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 encode
170	* the length of the part. If we are only interested in the transitive
171	* closure, then we can simply project out these coordinates first.
172	*/
173	static 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, type: isl_dim_in);
184	if (d < 0)
185	app = isl_map_free(map: app);
186	app = set_path_length(map: app, exactly: 0, length: 1);
187	app = isl_map_project_out(map: app, type: isl_dim_in, first: d, n: 1);
188	app = isl_map_project_out(map: app, type: isl_dim_out, first: d, n: 1);
189
190	app = isl_map_reset_space(map: app, space: isl_map_get_space(map));
191
192	test = isl_map_apply_range(map1: isl_map_copy(map), map2: isl_map_copy(map: app));
193	test = isl_map_union(map1: test, map2: isl_map_copy(map));
194
195	exact = isl_map_is_subset(map1: app, map2: test);
196
197	isl_map_free(map: app);
198	isl_map_free(map: test);
199
200	isl_map_free(map);
201
202	return exact;
203	}
204
205	/*
206	* The transitive closure implementation is based on the paper
207	* "Computing the Transitive Closure of a Union of Affine Integer
208	* Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
209	* Albert Cohen.
210	*/
211
212	/* Given a set of n offsets v_i (the rows of "steps"), construct a relation
213	* of the given dimension specification (Z^{n+1} -> Z^{n+1})
214	* that maps an element x to any element that can be reached
215	* by taking a non-negative number of steps along any of
216	* the extended offsets v'_i = [v_i 1].
217	* That is, construct
218	*
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 taken
222	* is equal to the difference in the final coordinates.
223	*/
224	static __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, type: isl_dim_in);
235	nparam = isl_space_dim(space, type: 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(space: isl_space_copy(space), extra: n, n_eq: d, n_ineq: n);
242
243	for (i = 0; i < n; ++i) {
244	k = isl_basic_map_alloc_div(bmap: 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(bmap: path, type: isl_dim_all);
252	if (total < 0)
253	goto error;
254	for (i = 0; i < d; ++i) {
255	k = isl_basic_map_alloc_equality(bmap: path);
256	if (k < 0)
257	goto error;
258	isl_seq_clr(p: path->eq[k], len: 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	else
265	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(bmap: path);
272	if (k < 0)
273	goto error;
274	isl_seq_clr(p: path->ineq[k], len: 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(bmap: path);
281	path = isl_basic_map_finalize(bmap: path);
282	return isl_map_from_basic_map(bmap: path);
283	error:
284	isl_space_free(space);
285	isl_basic_map_free(bmap: path);
286	return NULL;
287	}
288
289	#define IMPURE 0
290	#define PURE_PARAM 1
291	#define PURE_VAR 2
292	#define MIXED 3
293
294	/* Check whether the parametric constant term of constraint c is never
295	* positive in "bset".
296	*/
297	static 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, type: isl_dim_div);
309	d = isl_basic_set_dim(bset, type: isl_dim_set);
310	nparam = isl_basic_set_dim(bset, type: isl_dim_param);
311	total = isl_basic_set_dim(bset, type: 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(base: bset, n_eq: 0, n_ineq: 1);
318	k = isl_basic_set_alloc_inequality(bset);
319	if (k < 0)
320	goto error;
321	isl_seq_clr(p: bset->ineq[k], len: 1 + total);
322	isl_seq_cpy(dst: bset->ineq[k], src: c, len: 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;
334	error:
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 set
342	* variables are non-zero and if moreover the parametric constant
343	* can never attain positive values.
344	* Return IMPURE otherwise.
345	*/
346	static 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, type: isl_dim_div);
357	d = isl_basic_set_dim(bset, type: isl_dim_set);
358	nparam = isl_basic_set_dim(bset, type: 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(p: c + 1, len: nparam) == -1)
372	return PURE_VAR;
373	if (!v && isl_seq_first_non_zero(p: c + 1 + nparam, len: 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(dst: c, src: c, len: 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	*/
392	static __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, type: isl_dim_div);
401	d = isl_basic_set_dim(bset, type: isl_dim_set);
402	nparam = isl_basic_set_dim(bset, type: 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(p: bset->div[i] + 2, len: nparam) != -1)
417	p = 1;
418	if (isl_seq_first_non_zero(p: bset->div[i] + 2 + nparam, len: 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 identity
437	* mapping.
438	*/
439	static 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(bmap: path);
447	test = isl_basic_map_fix_si(bmap: test, type: isl_dim_div, pos, value: 0);
448	id = isl_basic_map_identity(space: isl_basic_map_get_space(bmap: path));
449	is_id = isl_basic_map_is_equal(bmap1: test, bmap2: id);
450	isl_basic_map_free(bmap: test);
451	isl_basic_map_free(bmap: id);
452	return is_id;
453	}
454
455	/* If any of the constraints is found to be impure then this function
456	* sets *impurity to 1.
457	*
458	* If impurity is NULL then we are dealing with a non-parametric set
459	* and so the constraints are obviously PURE_VAR.
460	*/
461	static __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(bset: delta, type: isl_dim_div);
472	total = isl_basic_map_dim(bmap: path, type: isl_dim_all);
473	if (n_div < 0 || total < 0)
474	return isl_basic_map_free(bmap: 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(bset: delta, c: 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(bmap: path);
489	if (k < 0)
490	goto error;
491	path_c = path->eq[k];
492	} else {
493	k = isl_basic_map_alloc_inequality(bmap: path);
494	if (k < 0)
495	goto error;
496	path_c = path->ineq[k];
497	}
498	isl_seq_clr(p: path_c, len: 1 + total);
499	if (p == PURE_VAR) {
500	isl_seq_cpy(dst: path_c + off,
501	src: delta_c[i] + 1 + nparam, len: d);
502	isl_int_set(path_c[off + d], delta_c[i][0]);
503	} else if (p == PURE_PARAM) {
504	isl_seq_cpy(dst: path_c, src: delta_c[i], len: 1 + nparam);
505	} else {
506	isl_seq_cpy(dst: path_c + off,
507	src: delta_c[i] + 1 + nparam, len: d);
508	isl_seq_cpy(dst: path_c, src: delta_c[i], len: 1 + nparam);
509	}
510	isl_seq_cpy(dst: path_c + off - n_div,
511	src: delta_c[i] + 1 + nparam + d, len: n_div);
512	}
513
514	return path;
515	error:
516	isl_basic_map_free(bmap: path);
517	return NULL;
518	}
519
520	/* Given a set of offsets "delta", construct a relation of the
521	* given dimension specification (Z^{n+1} -> Z^{n+1}) that
522	* is an overapproximation of the relations that
523	* maps an element x to any element that can be reached
524	* by taking a non-negative number of steps along any of
525	* the elements in "delta".
526	* That is, construct an approximation of
527	*
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 taken
532	* is equal to the difference in the final coordinates.
533	*
534	* In particular, let delta be defined as
535	*
536	* \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and
537	* C x + C'p + c >= 0 and
538	* D x + D'p + d >= 0 }
539	*
540	* where the constraints C x + C'p + c >= 0 are such that the parametric
541	* 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 as
543	*
544	* { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
545	* A f + k a >= 0 and B p + b >= 0 and
546	* 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 identity
550	* mapping, then we return
551	*
552	* { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
553	* A f + k a >= 0 and B p + b >= 0 and
554	* C f + C'p + c >= 0 and k >= 0 }
555	*
556	* instead.
557	*
558	* Existentially quantified variables in \delta are handled by
559	* classifying them as independent of the parameters, purely
560	* parameter dependent and others. Constraints containing
561	* 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 eliminate
565	* the parameters from \delta, resulting in a set
566	*
567	* \delta' = { [x] : E x + e >= 0 }
568	*
569	* and add the constraints
570	*
571	* E f + k e >= 0
572	*
573	* to the constructed relation.
574	*/
575	static __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(bset: delta, type: isl_dim_div);
590	d = isl_basic_set_dim(bset: delta, type: isl_dim_set);
591	nparam = isl_basic_set_dim(bset: delta, type: isl_dim_param);
592	if (n_div < 0 || d < 0 || nparam < 0)
593	goto error;
594	path = isl_basic_map_alloc_space(space: isl_space_copy(space), extra: n_div + d + 1,
595	n_eq: d + 1 + delta->n_eq, n_ineq: 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(bmap: 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(bmap: path, type: 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(bmap: path);
610	if (k < 0)
611	goto error;
612	isl_seq_clr(p: path->eq[k], len: 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(bset: delta);
619	if (n_div && !div_purity)
620	goto error;
621
622	path = add_delta_constraints(path, delta, off, nparam, d,
623	div_purity, eq: 1, impurity: &impurity);
624	path = add_delta_constraints(path, delta, off, nparam, d,
625	div_purity, eq: 0, impurity: &impurity);
626	if (impurity) {
627	isl_space *space = isl_basic_set_get_space(bset: delta);
628	delta = isl_basic_set_project_out(bset: delta,
629	type: isl_dim_param, first: 0, n: nparam);
630	delta = isl_basic_set_add_dims(bset: delta, type: isl_dim_param, n: nparam);
631	delta = isl_basic_set_reset_space(bset: delta, space);
632	if (!delta)
633	goto error;
634	path = isl_basic_map_extend_constraints(base: path, n_eq: delta->n_eq,
635	n_ineq: delta->n_ineq + 1);
636	path = add_delta_constraints(path, delta, off, nparam, d,
637	NULL, eq: 1, NULL);
638	path = add_delta_constraints(path, delta, off, nparam, d,
639	NULL, eq: 0, NULL);
640	path = isl_basic_map_gauss(bmap: path, NULL);
641	}
642
643	is_id = empty_path_is_identity(path, pos: n_div + d);
644	if (is_id < 0)
645	goto error;
646
647	k = isl_basic_map_alloc_inequality(bmap: path);
648	if (k < 0)
649	goto error;
650	isl_seq_clr(p: path->ineq[k], len: 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(ptr: div_purity);
656	isl_basic_set_free(bset: delta);
657	path = isl_basic_map_finalize(bmap: path);
658	if (is_id) {
659	isl_space_free(space);
660	return isl_map_from_basic_map(bmap: path);
661	}
662	return isl_basic_map_union(bmap1: path, bmap2: isl_basic_map_identity(space));
663	error:
664	free(ptr: div_purity);
665	isl_space_free(space);
666	isl_basic_set_free(bset: delta);
667	isl_basic_map_free(bmap: 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 difference
673	* in the final coordinates and imposes that this difference is positive.
674	* That is, construct
675	*
676	* { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
677	*/
678	static __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, type: isl_dim_in);
688	nparam = isl_space_dim(space, type: isl_dim_param);
689	total = isl_space_dim(space, type: 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, extra: 0, n_eq: 1, n_ineq: 1);
693	k = isl_basic_map_alloc_equality(bmap);
694	if (k < 0)
695	goto error;
696	isl_seq_clr(p: bmap->eq[k], len: 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(p: bmap->ineq[k], len: 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);
710	error:
711	isl_basic_map_free(bmap);
712	return NULL;
713	}
714
715	/* Check whether "path" is acyclic, where the last coordinates of domain
716	* and range of path encode the number of steps taken.
717	* That is, check whether
718	*
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	*/
724	static 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(map: path);
732	dim = isl_set_dim(set: delta, type: isl_dim_set);
733	if (dim < 0)
734	delta = isl_set_free(set: delta);
735	for (i = 0; i < dim; ++i) {
736	if (i == dim -1)
737	delta = isl_set_lower_bound_si(set: delta, type: isl_dim_set, pos: i, value: 1);
738	else
739	delta = isl_set_fix_si(set: delta, type: isl_dim_set, pos: i, value: 0);
740	}
741
742	acyclic = isl_set_is_empty(set: delta);
743	isl_set_free(set: delta);
744
745	return acyclic;
746	}
747
748	/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
749	* and a dimension specification (Z^{n+1} -> Z^{n+1}),
750	* construct a map that is an overapproximation of the map
751	* that takes an element from the space D \times Z to another
752	* element from the same space, such that the first n coordinates of the
753	* difference between them is a sum of differences between images
754	* and pre-images in one of the R_i and such that the last coordinate
755	* is equal to the number of steps taken.
756	* That is, let
757	*
758	* \Delta_i = { y - x | (x, y) in R_i }
759	*
760	* then the constructed map is an overapproximation of
761	*
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 the
766	* 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 overapproximation
769	* of the paths along elements of \Delta_i.
770	* Since each of these paths performs an addition, composition is
771	* symmetric and we can simply compose all resulting paths in any order.
772	*/
773	static __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, type: isl_dim_in);
782	if (d < 0)
783	goto error;
784
785	path = isl_map_identity(space: isl_space_copy(space));
786
787	steps = isl_mat_alloc(ctx: map->ctx, n_row: map->n, n_col: 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(bmap: isl_basic_map_copy(bmap: 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(bset: delta, dim: j,
801	val: &steps->row[n][j]);
802	if (fixed < 0) {
803	isl_basic_set_free(bset: delta);
804	goto error;
805	}
806	if (!fixed)
807	break;
808	}
809
810
811	if (j < d) {
812	path = isl_map_apply_range(map1: path,
813	map2: path_along_delta(space: isl_space_copy(space), delta));
814	path = isl_map_coalesce(map: path);
815	} else {
816	isl_basic_set_free(bset: delta);
817	++n;
818	}
819	}
820
821	if (n > 0) {
822	steps->n_row = n;
823	path = isl_map_apply_range(map1: path,
824	map2: path_along_steps(space: isl_space_copy(space), steps));
825	}
826
827	if (project && *project) {
828	*project = is_acyclic(path: isl_map_copy(map: path));
829	if (*project < 0)
830	goto error;
831	}
832
833	isl_space_free(space);
834	isl_mat_free(mat: steps);
835	return path;
836	error:
837	isl_space_free(space);
838	isl_mat_free(mat: steps);
839	isl_map_free(map: path);
840	return NULL;
841	}
842
843	static isl_bool isl_set_overlaps(__isl_keep isl_set *set1,
844	__isl_keep isl_set *set2)
845	{
846	return isl_bool_not(b: isl_set_is_disjoint(set1, set2));
847	}
848
849	/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
850	* and a dimension specification (Z^{n+1} -> Z^{n+1}),
851	* construct a map that is an overapproximation of the map
852	* that takes an element from the dom R \times Z to an
853	* element from ran R \times Z, such that the first n coordinates of the
854	* difference between them is a sum of differences between images
855	* and pre-images in one of the R_i and such that the last coordinate
856	* is equal to the number of steps taken.
857	* That is, let
858	*
859	* \Delta_i = { y - x | (x, y) in R_i }
860	*
861	* then the constructed map is an overapproximation of
862	*
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) and
865	* x in dom R and x + d in ran R and
866	* \sum_i k_i >= 1 }
867	*/
868	static __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(bmap: isl_map_copy(map));
879	domain = isl_set_coalesce(set: domain);
880	range = isl_map_range(map: isl_map_copy(map));
881	range = isl_set_coalesce(set: range);
882	overlaps = isl_set_overlaps(set1: domain, set2: range);
883	if (overlaps < 0 || !overlaps) {
884	isl_set_free(set: domain);
885	isl_set_free(set: 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, type: isl_dim_in, n: 1);
892	map = isl_map_add_dims(map, type: isl_dim_out, n: 1);
893	map = set_path_length(map, exactly: 1, length: 1);
894	return map;
895	}
896	app = isl_map_from_domain_and_range(domain, range);
897	app = isl_map_add_dims(map: app, type: isl_dim_in, n: 1);
898	app = isl_map_add_dims(map: app, type: isl_dim_out, n: 1);
899
900	check = exact && *exact == isl_bool_true;
901	path = construct_extended_path(space: isl_space_copy(space), map,
902	project: check ? &project : NULL);
903	app = isl_map_intersect(map1: app, map2: path);
904
905	if (check &&
906	(*exact = check_exactness(map: isl_map_copy(map), app: isl_map_copy(map: app),
907	project)) < 0)
908	goto error;
909
910	isl_space_free(space);
911	app = set_path_length(map: app, exactly: 0, length: 1);
912	return app;
913	error:
914	isl_space_free(space);
915	isl_map_free(map: app);
916	return NULL;
917	}
918
919	/* Call construct_component and, if "project" is set, project out
920	* the final coordinates.
921	*/
922	static __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, type: isl_dim_in);
932
933	app = construct_component(space, map, exact, project);
934	if (project) {
935	app = isl_map_project_out(map: app, type: isl_dim_in, first: d - 1, n: 1);
936	app = isl_map_project_out(map: app, type: isl_dim_out, first: d - 1, n: 1);
937	}
938	return app;
939	}
940
941	/* Compute an extended version, i.e., with path lengths, of
942	* an overapproximation of the transitive closure of "bmap"
943	* with path lengths greater than or equal to zero and with
944	* domain and range equal to "dom".
945	*/
946	static __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(set: dom, type: isl_dim_set, n: 1);
956	app = isl_map_from_domain_and_range(domain: dom, range: isl_set_copy(set: dom));
957	map = isl_map_from_basic_map(bmap: isl_basic_map_copy(bmap));
958	path = construct_extended_path(space, map, project: &project);
959	app = isl_map_intersect(map1: app, map2: path);
960
961	if ((*exact = check_exactness(map, app: isl_map_copy(map: app), project)) < 0)
962	goto error;
963
964	return app;
965	error:
966	isl_map_free(map: app);
967	return NULL;
968	}
969
970	/* Check whether qc has any elements of length at least one
971	* with domain and/or range outside of dom and ran.
972	*/
973	static 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(map: qc, type: isl_dim_in);
981	if (d < 0 || !dom || !ran)
982	return isl_bool_error;
983
984	qc = isl_map_copy(map: qc);
985	qc = set_path_length(map: qc, exactly: 0, length: 1);
986	qc = isl_map_project_out(map: qc, type: isl_dim_in, first: d - 1, n: 1);
987	qc = isl_map_project_out(map: qc, type: isl_dim_out, first: d - 1, n: 1);
988
989	s = isl_map_domain(bmap: isl_map_copy(map: qc));
990	subset = isl_set_is_subset(set1: s, set2: dom);
991	isl_set_free(set: s);
992	if (subset < 0)
993	goto error;
994	if (!subset) {
995	isl_map_free(map: qc);
996	return isl_bool_true;
997	}
998
999	s = isl_map_range(map: qc);
1000	subset = isl_set_is_subset(set1: s, set2: ran);
1001	isl_set_free(set: s);
1002
1003	return isl_bool_not(b: subset);
1004	error:
1005	isl_map_free(map: qc);
1006	return isl_bool_error;
1007	}
1008
1009	#define LEFT 2
1010	#define RIGHT 1
1011
1012	/* For each basic map in "map", except i, check whether it combines
1013	* with the transitive closure that is reflexive on C combines
1014	* to the left and to the right.
1015	*
1016	* In particular, if
1017	*
1018	* dom map_j \subseteq C
1019	*
1020	* then right[j] is set to 1. Otherwise, if
1021	*
1022	* ran map_i \cap dom map_j = \emptyset
1023	*
1024	* then right[j] is set to 0. Otherwise, composing to the right
1025	* is impossible.
1026	*
1027	* Similar, for composing to the left, we have if
1028	*
1029	* ran map_j \subseteq C
1030	*
1031	* then left[j] is set to 1. Otherwise, if
1032	*
1033	* dom map_i \cap ran map_j = \emptyset
1034	*
1035	* then left[j] is set to 0. Otherwise, composing to the left
1036	* is impossible.
1037	*
1038	* The return value is or'd with LEFT if composing to the left
1039	* is possible and with RIGHT if composing to the right is possible.
1040	*/
1041	static 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	bset: isl_basic_map_domain(
1058	bmap: isl_basic_map_copy(bmap: map->p[j])));
1059	if (!dom[j])
1060	return -1;
1061	overlaps = isl_set_overlaps(set1: ran[i], set2: dom[j]);
1062	if (overlaps < 0)
1063	return -1;
1064	if (!overlaps)
1065	right[j] = 0;
1066	else {
1067	subset = isl_set_is_subset(set1: dom[j], set2: C);
1068	if (subset < 0)
1069	return -1;
1070	if (subset)
1071	right[j] = 1;
1072	else
1073	ok &= ~RIGHT;
1074	}
1075	}
1076
1077	if (ok & LEFT) {
1078	if (!ran[j])
1079	ran[j] = isl_set_from_basic_set(
1080	bset: isl_basic_map_range(
1081	bmap: isl_basic_map_copy(bmap: map->p[j])));
1082	if (!ran[j])
1083	return -1;
1084	overlaps = isl_set_overlaps(set1: dom[i], set2: ran[j]);
1085	if (overlaps < 0)
1086	return -1;
1087	if (!overlaps)
1088	left[j] = 0;
1089	else {
1090	subset = isl_set_is_subset(set1: ran[j], set2: C);
1091	if (subset < 0)
1092	return -1;
1093	if (subset)
1094	left[j] = 1;
1095	else
1096	ok &= ~LEFT;
1097	}
1098	}
1099	}
1100
1101	return ok;
1102	}
1103
1104	static __isl_give isl_map *anonymize(__isl_take isl_map *map)
1105	{
1106	map = isl_map_reset(map, type: isl_dim_in);
1107	map = isl_map_reset(map, type: 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	*/
1115	static __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(space: 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(bmap: isl_basic_map_copy(bmap: map->p[j]));
1129	map_j = anonymize(map: map_j);
1130	if (left && left[j])
1131	map_j = isl_map_apply_range(map1: map_j, map2: isl_map_copy(map: qc));
1132	if (right && right[j])
1133	map_j = isl_map_apply_range(map1: isl_map_copy(map: qc), map2: map_j);
1134	comp = isl_map_union(map1: comp, map2: map_j);
1135	}
1136
1137	comp = isl_map_compute_divs(map: comp);
1138	comp = isl_map_coalesce(map: comp);
1139
1140	isl_map_free(map: qc);
1141
1142	return comp;
1143	}
1144
1145	/* Compute the transitive closure of "map" incrementally by
1146	* computing
1147	*
1148	* map_i^+ \cup qc^+
1149	*
1150	* or
1151	*
1152	* map_i^+ \cup ((id \cup map_i^) \circ qc^+)
1153	*
1154	* or
1155	*
1156	* map_i^+ \cup (qc^+ \circ (id \cup map_i^))
1157	*
1158	* depending on whether left or right are NULL.
1159	*/
1160	static __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(bmap: isl_basic_map_copy(bmap: map->p[i]));
1173	tc = construct_projected_component(space: isl_space_copy(space), map: map_i,
1174	exact, project: 1);
1175	isl_map_free(map: map_i);
1176
1177	if (*exact)
1178	qc = isl_map_transitive_closure(map: qc, exact);
1179
1180	if (!*exact) {
1181	isl_space_free(space);
1182	isl_map_free(map: tc);
1183	isl_map_free(map: qc);
1184	return isl_map_universe(space: isl_map_get_space(map));
1185	}
1186
1187	if (!left || !right)
1188	rtc = isl_map_union(map1: isl_map_copy(map: tc),
1189	map2: isl_map_identity(space: isl_map_get_space(map: tc)));
1190	if (!right)
1191	qc = isl_map_apply_range(map1: rtc, map2: qc);
1192	if (!left)
1193	qc = isl_map_apply_range(map1: qc, map2: rtc);
1194	qc = isl_map_union(map1: tc, map2: qc);
1195
1196	isl_space_free(space);
1197
1198	return qc;
1199	error:
1200	isl_space_free(space);
1201	isl_map_free(map: qc);
1202	return NULL;
1203	}
1204
1205	/* Given a map "map", try to find a basic map such that
1206	* map^+ can be computed as
1207	*
1208	* map^+ = map_i^+ \cup
1209	* \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 zero
1213	* 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 transitive
1216	* closures are exact and if the number of basic maps in the union,
1217	* after computing the integer divisions, is smaller than the number
1218	* of basic maps in the input map.
1219	*/
1220	static 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, type: isl_dim_in);
1232	if (d < 0)
1233	return isl_bool_error;
1234
1235	C = isl_set_union(set1: isl_map_domain(bmap: isl_map_copy(map)),
1236	set2: isl_map_range(map: isl_map_copy(map)));
1237	C = isl_set_from_basic_set(bset: isl_set_simple_hull(set: C));
1238	if (!C)
1239	return isl_bool_error;
1240	if (C->n != 1) {
1241	isl_set_free(set: 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(bset: isl_basic_map_domain(
1251	bmap: isl_basic_map_copy(bmap: map->p[i])));
1252	ran[i] = isl_set_from_basic_set(bset: isl_basic_map_range(
1253	bmap: isl_basic_map_copy(bmap: map->p[i])));
1254	qc = q_closure(space: isl_space_copy(space), dom: isl_set_copy(set: C),
1255	bmap: map->p[i], exact: &exact_i);
1256	if (!qc)
1257	goto error;
1258	if (!exact_i) {
1259	isl_map_free(map: qc);
1260	continue;
1261	}
1262	spurious = has_spurious_elements(qc, dom: dom[i], ran: ran[i]);
1263	if (spurious) {
1264	isl_map_free(map: qc);
1265	if (spurious < 0)
1266	goto error;
1267	continue;
1268	}
1269	qc = isl_map_project_out(map: qc, type: isl_dim_in, first: d, n: 1);
1270	qc = isl_map_project_out(map: qc, type: isl_dim_out, first: d, n: 1);
1271	qc = isl_map_compute_divs(map: 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(map: qc);
1279	continue;
1280	}
1281	*res = compute_incremental(space: isl_space_copy(space), map, i, qc,
1282	left, right, exact: &exact_i);
1283	if (!*res)
1284	goto error;
1285	if (exact_i)
1286	break;
1287	isl_map_free(map: *res);
1288	*res = NULL;
1289	}
1290
1291	isl_set_free(set: C);
1292
1293	return isl_bool_ok(b: *res != NULL);
1294	error:
1295	isl_set_free(set: C);
1296	return isl_bool_error;
1297	}
1298
1299	/* Try and compute the transitive closure of "map" as
1300	*
1301	* map^+ = map_i^+ \cup
1302	* \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 entire
1305	* map or the simple hull of domain and range of map_i.
1306	*/
1307	static __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, type: 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: &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	bset: isl_basic_map_domain(
1351	bmap: isl_basic_map_copy(bmap: map->p[i])));
1352	if (!dom[i])
1353	goto error;
1354	if (!ran[i])
1355	ran[i] = isl_set_from_basic_set(
1356	bset: isl_basic_map_range(
1357	bmap: isl_basic_map_copy(bmap: map->p[i])));
1358	if (!ran[i])
1359	goto error;
1360	C = isl_set_union(set1: isl_set_copy(set: dom[i]),
1361	set2: isl_set_copy(set: ran[i]));
1362	C = isl_set_from_basic_set(bset: isl_set_simple_hull(set: C));
1363	if (!C)
1364	goto error;
1365	if (C->n != 1) {
1366	isl_set_free(set: C);
1367	continue;
1368	}
1369	comp = composability(C, i, dom, ran, left, right, map);
1370	if (!comp || comp < 0) {
1371	isl_set_free(set: C);
1372	if (comp < 0)
1373	goto error;
1374	continue;
1375	}
1376	qc = q_closure(space: isl_space_copy(space), dom: C, bmap: map->p[i], exact: &exact_i);
1377	if (!qc)
1378	goto error;
1379	if (!exact_i) {
1380	isl_map_free(map: qc);
1381	continue;
1382	}
1383	spurious = has_spurious_elements(qc, dom: dom[i], ran: ran[i]);
1384	if (spurious) {
1385	isl_map_free(map: qc);
1386	if (spurious < 0)
1387	goto error;
1388	continue;
1389	}
1390	qc = isl_map_project_out(map: qc, type: isl_dim_in, first: d, n: 1);
1391	qc = isl_map_project_out(map: qc, type: isl_dim_out, first: d, n: 1);
1392	qc = isl_map_compute_divs(map: qc);
1393	qc = compose(map, i, qc, left: (comp & LEFT) ? left : NULL,
1394	right: (comp & RIGHT) ? right : NULL);
1395	if (!qc)
1396	goto error;
1397	if (qc->n >= map->n) {
1398	isl_map_free(map: qc);
1399	continue;
1400	}
1401	res = compute_incremental(space: isl_space_copy(space), map, i, qc,
1402	left: (comp & LEFT) ? left : NULL,
1403	right: (comp & RIGHT) ? right : NULL, exact: &exact_i);
1404	if (!res)
1405	goto error;
1406	if (exact_i)
1407	break;
1408	isl_map_free(map: res);
1409	res = NULL;
1410	}
1411
1412	for (i = 0; i < map->n; ++i) {
1413	isl_set_free(set: dom[i]);
1414	isl_set_free(set: ran[i]);
1415	}
1416	free(ptr: dom);
1417	free(ptr: ran);
1418	free(ptr: left);
1419	free(ptr: right);
1420
1421	if (res) {
1422	isl_space_free(space);
1423	return res;
1424	}
1425
1426	return construct_projected_component(space, map, exact, project);
1427	error:
1428	if (dom)
1429	for (i = 0; i < map->n; ++i)
1430	isl_set_free(set: dom[i]);
1431	free(ptr: dom);
1432	if (ran)
1433	for (i = 0; i < map->n; ++i)
1434	isl_set_free(set: ran[i]);
1435	free(ptr: ran);
1436	free(ptr: left);
1437	free(ptr: 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", into
1445	* 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 with
1448	* group[i] = i, as the other sets have already been combined
1449	* with one of the group leaders.
1450	*/
1451	static 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(set: 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(set1: set[i], set2: dom);
1465	if (o < 0)
1466	goto error;
1467	if (!o)
1468	continue;
1469
1470	set[i] = isl_set_union(set1: set[i], set2: 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(set: dom);
1479	return 0;
1480	error:
1481	isl_set_free(set: dom);
1482	return -1;
1483	}
1484
1485	/* Construct a map [x] -> [x+1], with parameters prescribed by "space".
1486	*/
1487	static __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, type: isl_dim_set, n: 1);
1495	space = isl_space_map_from_set(space);
1496	bmap = isl_basic_map_alloc_space(space, extra: 0, n_eq: 1, n_ineq: 0);
1497	total = isl_basic_map_dim(bmap, type: isl_dim_all);
1498	k = isl_basic_map_alloc_equality(bmap);
1499	if (total < 0 || k < 0)
1500	goto error;
1501	isl_seq_clr(p: bmap->eq[k], len: 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);
1506	error:
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 product
1512	* with the relation { [i] -> [i + 1] }.
1513	*/
1514	static 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(space: 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(map1: grid[i][j],
1529	map2: isl_map_copy(map: step));
1530
1531	isl_map_free(map: 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 whole
1540	* 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 diagonal
1543	* element corresponding to the current vertex is replaced by its
1544	* transitive closure to account for all indirect paths that stay
1545	* in the current vertex.
1546	*/
1547	static 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(map: grid[r][r],
1555	exact: 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	map1: isl_map_copy(map: grid[p][r]),
1566	map2: isl_map_copy(map: grid[r][q]));
1567	grid[p][q] = isl_map_union(map1: grid[p][q], map2: loop);
1568	loop = isl_map_apply_range(
1569	map1: isl_map_copy(map: grid[p][r]),
1570	map2: isl_map_apply_range(
1571	map1: isl_map_copy(map: grid[r][r]),
1572	map2: isl_map_copy(map: grid[r][q])));
1573	grid[p][q] = isl_map_union(map1: grid[p][q], map2: loop);
1574	grid[p][q] = isl_map_coalesce(map: 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 partition
1581	* as vertices.
1582	*
1583	* In particular, there are "n" elements in the partition and "group" is
1584	* 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 the
1588	* 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 the
1592	* path lengths encoded as the difference between an extra pair of
1593	* coordinates. We therefore apply the nested transitive closures
1594	* to relations that include these lengths. In particular, we replace
1595	* the input relation by the cross product with the unit length relation
1596	* { [i] -> [i + 1] }.
1597	*/
1598	static __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(ptr: 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(space: 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(map1: grid[i][j],
1629	map2: isl_map_from_basic_map(
1630	bmap: isl_basic_map_copy(bmap: 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(space: isl_map_get_space(map: grid[0][0]));
1639
1640	for (i = 0; i < n; ++i) {
1641	for (j = 0; j < n; ++j)
1642	app = isl_map_union(map1: app, map2: grid[i][j]);
1643	free(ptr: grid[i]);
1644	}
1645	free(ptr: grid);
1646
1647	free(ptr: group);
1648	isl_space_free(space);
1649
1650	return app;
1651	error:
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(map: grid[i][j]);
1658	free(ptr: grid[i]);
1659	}
1660	free(ptr: grid);
1661	free(ptr: group);
1662	isl_space_free(space);
1663	return NULL;
1664	}
1665
1666	/* Partition the domains and ranges of the n basic relations in list
1667	* into disjoint cells.
1668	*
1669	* To find the partition, we simply consider all of the domains
1670	* and ranges in turn and combine those that overlap.
1671	* "set" contains the partition elements and "group" indicates
1672	* 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 and
1677	* ranges in the corresponding group, or is equal to some l < k,
1678	* with l another domain or range in the same group.
1679	*/
1680	static 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(bset: isl_basic_map_domain(
1696	bmap: isl_basic_map_copy(bmap: list[i])));
1697	if (merge(set: *set, group, dom, pos: 2 * i) < 0)
1698	goto error;
1699	dom = isl_set_from_basic_set(bset: isl_basic_map_range(
1700	bmap: isl_basic_map_copy(bmap: list[i])));
1701	if (merge(set: *set, group, dom, pos: 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	} else
1714	group[i] = group[group[i]];
1715
1716	*n_group = g;
1717
1718	return group;
1719	error:
1720	if (*set) {
1721	for (i = 0; i < 2 * n; ++i)
1722	isl_set_free(set: (*set)[i]);
1723	free(ptr: *set);
1724	*set = NULL;
1725	}
1726	free(ptr: group);
1727	return NULL;
1728	}
1729
1730	/* Check if the domains and ranges of the basic maps in "map" can
1731	* be partitioned, and if so, apply Floyd-Warshall on the elements
1732	* of the partition. Note that we also apply this algorithm
1733	* 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 recursive
1735	* calls inside the Floyd-Warshall algorithm typically result in
1736	* non-linear path lengths quite quickly.
1737	*/
1738	static __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(ctx: map->ctx, list: map->p, n: map->n, set: &set, n_group: &n);
1752	if (!group)
1753	goto error;
1754
1755	for (i = 0; i < 2 * map->n; ++i)
1756	isl_set_free(set: set[i]);
1757
1758	free(ptr: set);
1759
1760	return floyd_warshall_with_groups(space, map, exact, project, group, n);
1761	error:
1762	isl_space_free(space);
1763	return NULL;
1764	}
1765
1766	/* Structure for representing the nodes of the graph of which
1767	* strongly connected components are being computed.
1768	*
1769	* list contains the actual nodes
1770	* check_closed is set if we may have used the fact that
1771	* a pair of basic maps can be interchanged
1772	*/
1773	struct isl_tc_follows_data {
1774	isl_basic_map **list;
1775	int check_closed;
1776	};
1777
1778	/* Check whether in the computation of the transitive closure
1779	* "list[i]" (R_1) should follow (or be part of the same component as)
1780	* "list[j]" (R_2).
1781	*
1782	* That is check whether
1783	*
1784	* R_1 \circ R_2
1785	*
1786	* is a subset of
1787	*
1788	* R_2 \circ R_1
1789	*
1790	* If so, then there is no reason for R_1 to immediately follow R_2
1791	* in any path.
1792	*
1793	* *check_closed is set if the subset relation holds while
1794	* R_1 \circ R_2 is not empty.
1795	*/
1796	static 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(bmap1: data->list[j], bmap2: 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	bmap: isl_basic_map_apply_range(
1811	bmap1: isl_basic_map_copy(bmap: data->list[j]),
1812	bmap2: isl_basic_map_copy(bmap: data->list[i])));
1813	subset = isl_map_is_empty(map: map21);
1814	if (subset < 0)
1815	goto error;
1816	if (subset) {
1817	isl_map_free(map: map21);
1818	return isl_bool_false;
1819	}
1820
1821	if (!isl_basic_map_is_transformation(bmap: data->list[i]) ||
1822	!isl_basic_map_is_transformation(bmap: data->list[j])) {
1823	isl_map_free(map: map21);
1824	return isl_bool_true;
1825	}
1826
1827	map12 = isl_map_from_basic_map(
1828	bmap: isl_basic_map_apply_range(
1829	bmap1: isl_basic_map_copy(bmap: data->list[i]),
1830	bmap2: isl_basic_map_copy(bmap: data->list[j])));
1831
1832	subset = isl_map_is_subset(map1: map21, map2: map12);
1833
1834	isl_map_free(map: map12);
1835	isl_map_free(map: map21);
1836
1837	if (subset)
1838	data->check_closed = 1;
1839
1840	return isl_bool_not(b: subset);
1841	error:
1842	isl_map_free(map: map21);
1843	return isl_bool_error;
1844	}
1845
1846	/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
1847	* and a dimension specification (Z^{n+1} -> Z^{n+1}),
1848	* construct a map that is an overapproximation of the map
1849	* that takes an element from the dom R \times Z to an
1850	* element from ran R \times Z, such that the first n coordinates of the
1851	* difference between them is a sum of differences between images
1852	* and pre-images in one of the R_i and such that the last coordinate
1853	* 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, let
1857	*
1858	* \Delta_i = { y - x | (x, y) in R_i }
1859	*
1860	* then the constructed map is an overapproximation of
1861	*
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) and
1864	* x in dom R and x + d in ran R }
1865	*
1866	* or
1867	*
1868	* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1869	* d = (\sum_i k_i \delta_i) and
1870	* 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, perform
1875	* the above on each component and then join the results in the correct
1876	* order, at each join also taking in the union of both arguments
1877	* to allow for paths that do not go through one of the two arguments.
1878	*/
1879	static __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(ctx: map->ctx, len: map->n, follows: &basic_map_follows, user: &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(space: isl_map_get_space(map));
1910	else
1911	path = isl_map_empty(space: isl_space_copy(space));
1912	path = anonymize(map: path);
1913	while (n) {
1914	struct isl_map *comp;
1915	isl_map *path_comp, *path_comb;
1916	comp = isl_map_alloc_space(space: isl_map_get_space(map), n, flags: 0);
1917	while (g->order[i] != -1) {
1918	comp = isl_map_add_basic_map(map: comp,
1919	bmap: isl_basic_map_copy(bmap: map->p[g->order[i]]));
1920	--n;
1921	++i;
1922	}
1923	path_comp = floyd_warshall(space: isl_space_copy(space),
1924	map: comp, exact, project);
1925	path_comp = anonymize(map: path_comp);
1926	path_comb = isl_map_apply_range(map1: isl_map_copy(map: path),
1927	map2: isl_map_copy(map: path_comp));
1928	path = isl_map_union(map1: path, map2: path_comp);
1929	path = isl_map_union(map1: path, map2: path_comb);
1930	isl_map_free(map: 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(map: path);
1939	if (closed < 0)
1940	goto error;
1941	if (!closed) {
1942	isl_tarjan_graph_free(g);
1943	isl_map_free(map: path);
1944	return floyd_warshall(space, map, exact: orig_exact, project);
1945	}
1946	}
1947
1948	isl_tarjan_graph_free(g);
1949	isl_space_free(space);
1950
1951	return path;
1952	error:
1953	isl_tarjan_graph_free(g);
1954	isl_space_free(space);
1955	isl_map_free(map: 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 map
1961	* that takes an element from the space D to another
1962	* element from the same space, such that the difference between
1963	* them is a strictly positive sum of differences between images
1964	* 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, let
1967	*
1968	* \Delta_i = { y - x | (x, y) in R_i }
1969	*
1970	* then the constructed map is an overapproximation of
1971	*
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	* or
1975	*
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, then
1982	* we construct an extended mapping with an extra coordinate
1983	* that indicates the number of steps taken. In particular,
1984	* the difference in the last coordinate is equal to the number
1985	* of steps taken to move from a domain element to the corresponding
1986	* image element(s).
1987	*/
1988	static __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, type: isl_dim_in, n: 1);
2000	space = isl_space_add_dims(space, type: isl_dim_out, n: 1);
2001
2002	app = construct_power_components(space: 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 transitive
2014	* closure, so we can use a more relaxed exactness check.
2015	* The lengths of the paths are also projected out instead of being
2016	* encoded as the difference between an extra pair of final coordinates.
2017	*/
2018	static __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 lengths
2039	* encoded as the difference between the final coordinates.
2040	* In the final step, this difference is equated to an extra parameter
2041	* and made positive. The extra coordinates are subsequently projected out
2042	* 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, type: isl_dim_in);
2053	param = isl_map_dim(map, type: 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(set: isl_map_wrap(map));
2062	map = isl_map_add_dims(map, type: isl_dim_in, n: 1);
2063	map = isl_map_set_dim_name(map, type: isl_dim_in, pos: 0, s: "k");
2064	return map;
2065	}
2066
2067	target_space = isl_map_get_space(map);
2068	target_space = isl_space_from_range(space: isl_space_wrap(space: target_space));
2069	target_space = isl_space_add_dims(space: target_space, type: isl_dim_in, n: 1);
2070	target_space = isl_space_set_dim_name(space: target_space, type: isl_dim_in, pos: 0, name: "k");
2071
2072	map = map_power(map, exact, project: 0);
2073
2074	map = isl_map_add_dims(map, type: isl_dim_param, n: 1);
2075	space = isl_map_get_space(map);
2076	diff = equate_parameter_to_length(space, param);
2077	map = isl_map_intersect(map1: map, map2: diff);
2078	map = isl_map_project_out(map, type: isl_dim_in, first: d, n: 1);
2079	map = isl_map_project_out(map, type: isl_dim_out, first: d, n: 1);
2080	map = isl_map_from_range(set: isl_map_wrap(map));
2081	map = isl_map_move_dims(map, dst_type: isl_dim_in, dst_pos: 0, src_type: isl_dim_param, src_pos: param, n: 1);
2082
2083	map = isl_map_reset_space(map, space: target_space);
2084
2085	return map;
2086	}
2087
2088	/* Compute a relation that maps each element in the range of the input
2089	* relation to the lengths of all paths composed of edges in the input
2090	* 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, the
2094	* range has become the domain and the exponent is the range instead
2095	* 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, type: isl_dim_in);
2106	param = isl_map_dim(map, type: 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, type: isl_dim_out, first: 0, n: d);
2117	map = isl_map_add_dims(map, type: isl_dim_out, n: 1);
2118	return map;
2119	}
2120
2121	map = map_power(map, exact, project: 0);
2122
2123	map = isl_map_add_dims(map, type: isl_dim_param, n: 1);
2124	space = isl_map_get_space(map);
2125	diff = equate_parameter_to_length(space, param);
2126	map = isl_map_intersect(map1: map, map2: diff);
2127	map = isl_map_project_out(map, type: isl_dim_in, first: 0, n: d + 1);
2128	map = isl_map_project_out(map, type: isl_dim_out, first: d, n: 1);
2129	map = isl_map_reverse(map);
2130	map = isl_map_move_dims(map, dst_type: isl_dim_out, dst_pos: 0, src_type: isl_dim_param, src_pos: param, n: 1);
2131
2132	return map;
2133	}
2134
2135	/* Given a map, compute the smallest superset of this map that is of the form
2136	*
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 with
2146	* the given domain and range.
2147	*
2148	* If with_id is set, then try to include as much of the identity mapping
2149	* as possible, by computing
2150	*
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 set
2157	*
2158	* delta = { j - i | i -> j in map },
2159	*
2160	* look for stride constraint on the individual dimensions and compute
2161	* (constant) lower and upper bounds for each individual dimension,
2162	* adding a constraint for each bound not equal to infinity.
2163	*/
2164	static __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(map: isl_map_copy(map));
2183
2184	aff = isl_set_affine_hull(set: isl_set_copy(set: delta));
2185	if (!aff)
2186	goto error;
2187	space = isl_map_get_space(map);
2188	d = isl_space_dim(space, type: isl_dim_in);
2189	nparam = isl_space_dim(space, type: isl_dim_param);
2190	total = isl_space_dim(space, type: isl_dim_all);
2191	bmap = isl_basic_map_alloc_space(space,
2192	extra: aff->n_div + 1, n_eq: aff->n_div, n_ineq: 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(bset: aff, i))
2201	continue;
2202	k = isl_basic_map_alloc_equality(bmap);
2203	if (k < 0)
2204	goto error;
2205	isl_seq_clr(p: bmap->eq[k], len: 1 + nparam);
2206	isl_seq_cpy(dst: bmap->eq[k] + 1 + nparam + d,
2207	src: aff->eq[i] + 1 + nparam, len: d);
2208	isl_seq_neg(dst: bmap->eq[k] + 1 + nparam,
2209	src: aff->eq[i] + 1 + nparam, len: d);
2210	isl_seq_cpy(dst: bmap->eq[k] + 1 + nparam + 2 * d,
2211	src: aff->eq[i] + 1 + nparam + d, len: aff->n_div);
2212	isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);
2213	}
2214	obj = isl_vec_alloc(ctx: map->ctx, size: 1 + nparam + d);
2215	if (!obj)
2216	goto error;
2217	isl_seq_clr(p: obj->el, len: 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(set: delta, max: 0, f: obj->el, denom: map->ctx->one, opt: &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(p: bmap->ineq[k],
2232	len: 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(set: delta, max: 1, f: obj->el, denom: map->ctx->one, opt: &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(p: bmap->ineq[k],
2247	len: 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(p: bmap->ineq[k],
2259	len: 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(domain: dom, range: ran);
2265
2266	isl_vec_free(vec: obj);
2267	isl_basic_set_free(bset: aff);
2268	isl_map_free(map);
2269	bmap = isl_basic_map_finalize(bmap);
2270	isl_set_free(set: delta);
2271	isl_int_clear(opt);
2272
2273	map = isl_map_from_basic_map(bmap);
2274	map = isl_map_intersect(map1: map, map2: app);
2275
2276	return map;
2277	error:
2278	isl_vec_free(vec: obj);
2279	isl_basic_map_free(bmap);
2280	isl_basic_set_free(bset: aff);
2281	isl_set_free(set: dom);
2282	isl_set_free(set: ran);
2283	isl_map_free(map);
2284	isl_set_free(set: 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 form
2290	*
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 with
2300	* domain and range of the original map.
2301	*/
2302	static __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(bmap: isl_map_copy(map));
2308	domain = isl_set_coalesce(set: domain);
2309	range = isl_map_range(map: isl_map_copy(map));
2310	range = isl_set_coalesce(set: range);
2311
2312	return box_closure_on_domain(map, dom: domain, ran: range, with_id: 0);
2313	}
2314
2315	/* Given a map, compute the smallest superset of this map that is of the form
2316	*
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 with
2326	* the given domain.
2327	*/
2328	static __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, ran: isl_set_copy(set: dom), with_id: 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 that
2337	*
2338	* app \subset (map \cup (map \circ app))
2339	*/
2340	static 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(map: isl_map_copy(map: app));
2350	d = isl_set_dim(set: delta, type: isl_dim_set);
2351	if (d < 0)
2352	delta = isl_set_free(set: delta);
2353	for (i = 0; i < d; ++i)
2354	delta = isl_set_fix_si(set: delta, type: isl_dim_set, pos: i, value: 0);
2355	is_empty = isl_set_is_empty(set: delta);
2356	isl_set_free(set: delta);
2357	if (is_empty < 0 || !is_empty)
2358	return is_empty;
2359
2360	test = isl_map_apply_range(map1: isl_map_copy(map: app), map2: isl_map_copy(map));
2361	test = isl_map_union(map1: test, map2: isl_map_copy(map));
2362	is_exact = isl_map_is_subset(map1: app, map2: test);
2363	isl_map_free(map: test);
2364
2365	return is_exact;
2366	}
2367
2368	/* Check if basic map M_i can be combined with all the other
2369	* basic maps such that
2370	*
2371	* (\cup_j M_j)^+
2372	*
2373	* can be computed as
2374	*
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 representation
2378	* of
2379	*
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 paths
2384	* of length zero, i.e., the result of box_closure(., 1).
2385	* The criterion, as proposed by Kelly et al., is that
2386	* id = M_i^? - M_i^+ can be represented as a basic map
2387	* and that
2388	*
2389	* id \circ M_j \circ id = M_j
2390	*
2391	* for each j != i.
2392	*
2393	* If this function returns 1, then tc and qc are set to
2394	* M_i^+ and M_i^?, respectively.
2395	*/
2396	static 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(set1: isl_map_domain(bmap: isl_map_copy(map)),
2407	set2: isl_map_range(map: isl_map_copy(map)));
2408	C = isl_set_from_basic_set(bset: isl_set_simple_hull(set: C));
2409	if (!C)
2410	goto error;
2411
2412	map_i = isl_map_from_basic_map(bmap: isl_basic_map_copy(bmap: map->p[i]));
2413	*tc = box_closure(map: isl_map_copy(map: map_i));
2414	*qc = box_closure_with_identity(map: map_i, dom: C);
2415	id = isl_map_subtract(map1: isl_map_copy(map: *qc), map2: isl_map_copy(map: *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	bmap: isl_basic_map_copy(bmap: map->p[j]));
2430	test = isl_map_apply_range(map1: isl_map_copy(map: id),
2431	map2: isl_map_copy(map: map_j));
2432	test = isl_map_apply_range(map1: test, map2: isl_map_copy(map: id));
2433	is_ok = isl_map_is_equal(map1: test, map2: map_j);
2434	isl_map_free(map: map_j);
2435	isl_map_free(map: test);
2436	if (is_ok < 0)
2437	goto error;
2438	if (!is_ok)
2439	break;
2440	}
2441
2442	done:
2443	isl_map_free(map: id);
2444	if (j == map->n)
2445	return 1;
2446
2447	isl_map_free(map: *qc);
2448	isl_map_free(map: *tc);
2449	*qc = NULL;
2450	*tc = NULL;
2451
2452	return 0;
2453	error:
2454	isl_map_free(map: id);
2455	isl_map_free(map: *qc);
2456	isl_map_free(map: *tc);
2457	*qc = NULL;
2458	*tc = NULL;
2459	return -1;
2460	}
2461
2462	static __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(map: 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(map: app);
2473	else
2474	*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 from
2483	* "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 and
2487	* compute
2488	*
2489	* (\cup_j M_j)^+
2490	*
2491	* as
2492	*
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	*/
2499	static __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: &tc, qc: &qc);
2515	if (ok < 0)
2516	goto error;
2517	if (!ok)
2518	continue;
2519
2520	app = isl_map_alloc_space(space: isl_map_get_space(map), n: map->n - 1, flags: 0);
2521
2522	for (j = 0; j < map->n; ++j) {
2523	if (j == i)
2524	continue;
2525	app = isl_map_add_basic_map(map: app,
2526	bmap: isl_basic_map_copy(bmap: map->p[j]));
2527	}
2528
2529	app = isl_map_apply_range(map1: isl_map_copy(map: qc), map2: app);
2530	app = isl_map_apply_range(map1: app, map2: qc);
2531
2532	app = isl_map_union(map1: tc, map2: transitive_closure_omega(map: 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(map: app);
2541	if (exact_i < 0)
2542	goto error;
2543	}
2544
2545	return box_closure_with_check(map, exact);
2546	error:
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 tell
2554	* it to project out the lengths of the paths instead of equating
2555	* 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, project: 1);
2582	map = isl_map_reset_space(map, space: target_dim);
2583
2584	return map;
2585	error:
2586	isl_map_free(map);
2587	return NULL;
2588	}
2589
2590	static 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
2601	static 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(bmap: map->p[i]);
2608	if (!**next)
2609	goto error;
2610	(*next)++;
2611	}
2612
2613	isl_map_free(map);
2614	return isl_stat_ok;
2615	error:
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 the
2623	* grid have domains and ranges of the same dimension and so
2624	* the standard algorithm can be used because the nested transitive
2625	* closures are only applied to diagonal elements and because all
2626	* compositions are performed on relations with compatible domains and ranges.
2627	*/
2628	static __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: &set, n_group: &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(space: isl_set_get_space(set: set[i]));
2652	space2 = isl_set_get_space(set: set[j]);
2653	space = isl_space_join(left: space1, right: 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(map1: grid[i][j],
2662	map2: isl_map_from_basic_map(
2663	bmap: isl_basic_map_copy(bmap: list[k])));
2664	}
2665
2666	floyd_warshall_iterate(grid, n: n_group, exact);
2667
2668	app = isl_union_map_empty(space: isl_map_get_space(map: 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(umap: app, map: grid[i][j]);
2673	free(ptr: grid[i]);
2674	}
2675	free(ptr: grid);
2676
2677	for (i = 0; i < 2 * n; ++i)
2678	isl_set_free(set: set[i]);
2679	free(ptr: set);
2680
2681	free(ptr: group);
2682	return app;
2683	error:
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(map: grid[i][j]);
2690	free(ptr: grid[i]);
2691	}
2692	free(ptr: grid);
2693	if (set) {
2694	for (i = 0; i < 2 * n; ++i)
2695	isl_set_free(set: set[i]);
2696	free(ptr: set);
2697	}
2698	free(ptr: 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 map
2706	* and then perform the algorithm on this list.
2707	*/
2708	static __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, fn: inc_count, user: &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, fn: collect_basic_map, user: &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(bmap: list[i]);
2735	free(ptr: list);
2736	}
2737
2738	isl_union_map_free(umap);
2739	return res;
2740	error:
2741	if (list) {
2742	for (i = 0; i < n; ++i)
2743	isl_basic_map_free(bmap: list[i]);
2744	free(ptr: 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 of
2752	* construct_power_components with the major difference that all
2753	* operations are performed on union maps.
2754	*/
2755	static __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, fn: inc_count, user: &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, fn: collect_basic_map, user: &next) < 0)
2785	goto error;
2786
2787	data.list = list;
2788	data.check_closed = 0;
2789	g = isl_tarjan_graph_init(ctx, len: n, follows: &basic_map_follows, user: &data);
2790	if (!g)
2791	goto error;
2792
2793	c = 0;
2794	i = 0;
2795	l = n;
2796	path = isl_union_map_empty(space: 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(space: isl_union_map_get_space(umap));
2801	while (g->order[i] != -1) {
2802	comp = isl_union_map_add_map(umap: comp,
2803	map: isl_map_from_basic_map(
2804	bmap: isl_basic_map_copy(bmap: list[g->order[i]])));
2805	--l;
2806	++i;
2807	}
2808	path_comp = union_floyd_warshall(umap: comp, exact);
2809	path_comb = isl_union_map_apply_range(umap1: isl_union_map_copy(umap: path),
2810	umap2: isl_union_map_copy(umap: path_comp));
2811	path = isl_union_map_union(umap1: path, umap2: path_comp);
2812	path = isl_union_map_union(umap1: path, umap2: 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(umap: 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(bmap: list[i]);
2830	free(ptr: list);
2831
2832	if (recheck) {
2833	isl_union_map_free(umap: path);
2834	return union_floyd_warshall(umap, exact);
2835	}
2836
2837	isl_union_map_free(umap);
2838
2839	return path;
2840	error:
2841	isl_tarjan_graph_free(g);
2842	if (list) {
2843	for (i = 0; i < n; ++i)
2844	isl_basic_map_free(bmap: list[i]);
2845	free(ptr: list);
2846	}
2847	isl_union_map_free(umap);
2848	isl_union_map_free(umap: 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;
2875	error:
2876	isl_union_map_free(umap);
2877	return NULL;
2878	}
2879
2880	struct isl_union_power {
2881	isl_union_map *pow;
2882	isl_bool *exact;
2883	};
2884
2885	static 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, exact: 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	*/
2897	static __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, type: isl_dim_in, n: 1);
2902	space = isl_space_add_dims(space, type: isl_dim_out, n: 1);
2903	bmap = isl_basic_map_universe(space);
2904	bmap = isl_basic_map_deltas_map(bmap);
2905
2906	return isl_union_map_from_map(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, fn: &power, user: &up);
2928	isl_union_map_free(umap);
2929	return up.pow;
2930	}
2931	inc = isl_union_map_from_map(map: increment(space: isl_union_map_get_space(umap)));
2932	umap = isl_union_map_product(umap1: inc, umap2: umap);
2933	umap = isl_union_map_transitive_closure(umap, exact);
2934	umap = isl_union_map_zip(umap);
2935	dm = deltas_map(space: isl_union_map_get_space(umap));
2936	umap = isl_union_map_apply_domain(umap1: umap, umap2: dm);
2937
2938	return umap;
2939	}
2940
2941	#undef TYPE
2942	#define TYPE isl_map
2943	#include "isl_power_templ.c"
2944
2945	#undef TYPE
2946	#define TYPE isl_union_map
2947	#include "isl_power_templ.c"
2948