1/*
2 * Copyright 2008-2009 Katholieke Universiteit Leuven
3 *
4 * Use of this software is governed by the MIT license
5 *
6 * Written by Sven Verdoolaege, K.U.Leuven, Departement
7 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
8 */
9
10#include <isl_map_private.h>
11#include <isl/seq.h>
12#include <isl/set.h>
13#include <isl/map.h>
14#include "isl_tab.h"
15#include <isl_point_private.h>
16
17/* Expand the constraint "c" into "v".  The initial "dim" dimensions
18 * are the same, but "v" may have more divs than "c" and the divs of "c"
19 * may appear in different positions in "v".
20 * The number of divs in "c" is given by "n_div" and the mapping
21 * of divs in "c" to divs in "v" is given by "div_map".
22 *
23 * Although it shouldn't happen in practice, it is theoretically
24 * possible that two or more divs in "c" are mapped to the same div in "v".
25 * These divs are then necessarily the same, so we simply add their
26 * coefficients.
27 */
28static void expand_constraint(isl_vec *v, unsigned dim,
29	isl_int *c, int *div_map, unsigned n_div)
30{
31	int i;
32
33	isl_seq_cpy(v->el, c, 1 + dim);
34	isl_seq_clr(v->el + 1 + dim, v->size - (1 + dim));
35
36	for (i = 0; i < n_div; ++i) {
37		int pos = 1 + dim + div_map[i];
38		isl_int_add(v->el[pos], v->el[pos], c[1 + dim + i]);
39	}
40}
41
42/* Add all constraints of bmap to tab.  The equalities of bmap
43 * are added as a pair of inequalities.
44 */
45static int tab_add_constraints(struct isl_tab *tab,
46	__isl_keep isl_basic_map *bmap, int *div_map)
47{
48	int i;
49	unsigned dim;
50	unsigned tab_total;
51	unsigned bmap_total;
52	isl_vec *v;
53
54	if (!tab || !bmap)
55		return -1;
56
57	tab_total = isl_basic_map_total_dim(tab->bmap);
58	bmap_total = isl_basic_map_total_dim(bmap);
59	dim = isl_space_dim(tab->bmap->dim, isl_dim_all);
60
61	if (isl_tab_extend_cons(tab, 2 * bmap->n_eq + bmap->n_ineq) < 0)
62		return -1;
63
64	v = isl_vec_alloc(bmap->ctx, 1 + tab_total);
65	if (!v)
66		return -1;
67
68	for (i = 0; i < bmap->n_eq; ++i) {
69		expand_constraint(v, dim, bmap->eq[i], div_map, bmap->n_div);
70		if (isl_tab_add_ineq(tab, v->el) < 0)
71			goto error;
72		isl_seq_neg(bmap->eq[i], bmap->eq[i], 1 + bmap_total);
73		expand_constraint(v, dim, bmap->eq[i], div_map, bmap->n_div);
74		if (isl_tab_add_ineq(tab, v->el) < 0)
75			goto error;
76		isl_seq_neg(bmap->eq[i], bmap->eq[i], 1 + bmap_total);
77		if (tab->empty)
78			break;
79	}
80
81	for (i = 0; i < bmap->n_ineq; ++i) {
82		expand_constraint(v, dim, bmap->ineq[i], div_map, bmap->n_div);
83		if (isl_tab_add_ineq(tab, v->el) < 0)
84			goto error;
85		if (tab->empty)
86			break;
87	}
88
89	isl_vec_free(v);
90	return 0;
91error:
92	isl_vec_free(v);
93	return -1;
94}
95
96/* Add a specific constraint of bmap (or its opposite) to tab.
97 * The position of the constraint is specified by "c", where
98 * the equalities of bmap are counted twice, once for the inequality
99 * that is equal to the equality, and once for its negation.
100 */
101static int tab_add_constraint(struct isl_tab *tab,
102	__isl_keep isl_basic_map *bmap, int *div_map, int c, int oppose)
103{
104	unsigned dim;
105	unsigned tab_total;
106	unsigned bmap_total;
107	isl_vec *v;
108	int r;
109
110	if (!tab || !bmap)
111		return -1;
112
113	tab_total = isl_basic_map_total_dim(tab->bmap);
114	bmap_total = isl_basic_map_total_dim(bmap);
115	dim = isl_space_dim(tab->bmap->dim, isl_dim_all);
116
117	v = isl_vec_alloc(bmap->ctx, 1 + tab_total);
118	if (!v)
119		return -1;
120
121	if (c < 2 * bmap->n_eq) {
122		if ((c % 2) != oppose)
123			isl_seq_neg(bmap->eq[c/2], bmap->eq[c/2],
124					1 + bmap_total);
125		if (oppose)
126			isl_int_sub_ui(bmap->eq[c/2][0], bmap->eq[c/2][0], 1);
127		expand_constraint(v, dim, bmap->eq[c/2], div_map, bmap->n_div);
128		r = isl_tab_add_ineq(tab, v->el);
129		if (oppose)
130			isl_int_add_ui(bmap->eq[c/2][0], bmap->eq[c/2][0], 1);
131		if ((c % 2) != oppose)
132			isl_seq_neg(bmap->eq[c/2], bmap->eq[c/2],
133					1 + bmap_total);
134	} else {
135		c -= 2 * bmap->n_eq;
136		if (oppose) {
137			isl_seq_neg(bmap->ineq[c], bmap->ineq[c],
138					1 + bmap_total);
139			isl_int_sub_ui(bmap->ineq[c][0], bmap->ineq[c][0], 1);
140		}
141		expand_constraint(v, dim, bmap->ineq[c], div_map, bmap->n_div);
142		r = isl_tab_add_ineq(tab, v->el);
143		if (oppose) {
144			isl_int_add_ui(bmap->ineq[c][0], bmap->ineq[c][0], 1);
145			isl_seq_neg(bmap->ineq[c], bmap->ineq[c],
146					1 + bmap_total);
147		}
148	}
149
150	isl_vec_free(v);
151	return r;
152}
153
154static int tab_add_divs(struct isl_tab *tab, __isl_keep isl_basic_map *bmap,
155	int **div_map)
156{
157	int i, j;
158	struct isl_vec *vec;
159	unsigned total;
160	unsigned dim;
161
162	if (!bmap)
163		return -1;
164	if (!bmap->n_div)
165		return 0;
166
167	if (!*div_map)
168		*div_map = isl_alloc_array(bmap->ctx, int, bmap->n_div);
169	if (!*div_map)
170		return -1;
171
172	total = isl_basic_map_total_dim(tab->bmap);
173	dim = total - tab->bmap->n_div;
174	vec = isl_vec_alloc(bmap->ctx, 2 + total + bmap->n_div);
175	if (!vec)
176		return -1;
177
178	for (i = 0; i < bmap->n_div; ++i) {
179		isl_seq_cpy(vec->el, bmap->div[i], 2 + dim);
180		isl_seq_clr(vec->el + 2 + dim, tab->bmap->n_div);
181		for (j = 0; j < i; ++j)
182			isl_int_set(vec->el[2 + dim + (*div_map)[j]],
183					bmap->div[i][2 + dim + j]);
184		for (j = 0; j < tab->bmap->n_div; ++j)
185			if (isl_seq_eq(tab->bmap->div[j],
186					vec->el, 2 + dim + tab->bmap->n_div))
187				break;
188		(*div_map)[i] = j;
189		if (j == tab->bmap->n_div) {
190			vec->size = 2 + dim + tab->bmap->n_div;
191			if (isl_tab_add_div(tab, vec, NULL, NULL) < 0)
192				goto error;
193		}
194	}
195
196	isl_vec_free(vec);
197
198	return 0;
199error:
200	isl_vec_free(vec);
201
202	return -1;
203}
204
205/* Freeze all constraints of tableau tab.
206 */
207static int tab_freeze_constraints(struct isl_tab *tab)
208{
209	int i;
210
211	for (i = 0; i < tab->n_con; ++i)
212		if (isl_tab_freeze_constraint(tab, i) < 0)
213			return -1;
214
215	return 0;
216}
217
218/* Check for redundant constraints starting at offset.
219 * Put the indices of the redundant constraints in index
220 * and return the number of redundant constraints.
221 */
222static int n_non_redundant(isl_ctx *ctx, struct isl_tab *tab,
223	int offset, int **index)
224{
225	int i, n;
226	int n_test = tab->n_con - offset;
227
228	if (isl_tab_detect_redundant(tab) < 0)
229		return -1;
230
231	if (n_test == 0)
232		return 0;
233	if (!*index)
234		*index = isl_alloc_array(ctx, int, n_test);
235	if (!*index)
236		return -1;
237
238	for (n = 0, i = 0; i < n_test; ++i) {
239		int r;
240		r = isl_tab_is_redundant(tab, offset + i);
241		if (r < 0)
242			return -1;
243		if (r)
244			continue;
245		(*index)[n++] = i;
246	}
247
248	return n;
249}
250
251/* basic_map_collect_diff calls add on each of the pieces of
252 * the set difference between bmap and map until the add method
253 * return a negative value.
254 */
255struct isl_diff_collector {
256	int (*add)(struct isl_diff_collector *dc,
257		    __isl_take isl_basic_map *bmap);
258};
259
260/* Compute the set difference between bmap and map and call
261 * dc->add on each of the piece until this function returns
262 * a negative value.
263 * Return 0 on success and -1 on error.  dc->add returning
264 * a negative value is treated as an error, but the calling
265 * function can interpret the results based on the state of dc.
266 *
267 * Assumes that map has known divs.
268 *
269 * The difference is computed by a backtracking algorithm.
270 * Each level corresponds to a basic map in "map".
271 * When a node in entered for the first time, we check
272 * if the corresonding basic map intersects the current piece
273 * of "bmap".  If not, we move to the next level.
274 * Otherwise, we split the current piece into as many
275 * pieces as there are non-redundant constraints of the current
276 * basic map in the intersection.  Each of these pieces is
277 * handled by a child of the current node.
278 * In particular, if there are n non-redundant constraints,
279 * then for each 0 <= i < n, a piece is cut off by adding
280 * constraints 0 <= j < i and adding the opposite of constraint i.
281 * If there are no non-redundant constraints, meaning that the current
282 * piece is a subset of the current basic map, then we simply backtrack.
283 *
284 * In the leaves, we check if the remaining piece has any integer points
285 * and if so, pass it along to dc->add.  As a special case, if nothing
286 * has been removed when we end up in a leaf, we simply pass along
287 * the original basic map.
288 */
289static int basic_map_collect_diff(__isl_take isl_basic_map *bmap,
290	__isl_take isl_map *map, struct isl_diff_collector *dc)
291{
292	int i;
293	int modified;
294	int level;
295	int init;
296	int empty;
297	isl_ctx *ctx;
298	struct isl_tab *tab = NULL;
299	struct isl_tab_undo **snap = NULL;
300	int *k = NULL;
301	int *n = NULL;
302	int **index = NULL;
303	int **div_map = NULL;
304
305	empty = isl_basic_map_is_empty(bmap);
306	if (empty) {
307		isl_basic_map_free(bmap);
308		isl_map_free(map);
309		return empty < 0 ? -1 : 0;
310	}
311
312	bmap = isl_basic_map_cow(bmap);
313	map = isl_map_cow(map);
314
315	if (!bmap || !map)
316		goto error;
317
318	ctx = map->ctx;
319	snap = isl_alloc_array(map->ctx, struct isl_tab_undo *, map->n);
320	k = isl_alloc_array(map->ctx, int, map->n);
321	n = isl_alloc_array(map->ctx, int, map->n);
322	index = isl_calloc_array(map->ctx, int *, map->n);
323	div_map = isl_calloc_array(map->ctx, int *, map->n);
324	if (!snap || !k || !n || !index || !div_map)
325		goto error;
326
327	bmap = isl_basic_map_order_divs(bmap);
328	map = isl_map_order_divs(map);
329
330	tab = isl_tab_from_basic_map(bmap, 1);
331	if (!tab)
332		goto error;
333
334	modified = 0;
335	level = 0;
336	init = 1;
337
338	while (level >= 0) {
339		if (level >= map->n) {
340			int empty;
341			struct isl_basic_map *bm;
342			if (!modified) {
343				if (dc->add(dc, isl_basic_map_copy(bmap)) < 0)
344					goto error;
345				break;
346			}
347			bm = isl_basic_map_copy(tab->bmap);
348			bm = isl_basic_map_cow(bm);
349			bm = isl_basic_map_update_from_tab(bm, tab);
350			bm = isl_basic_map_simplify(bm);
351			bm = isl_basic_map_finalize(bm);
352			empty = isl_basic_map_is_empty(bm);
353			if (empty)
354				isl_basic_map_free(bm);
355			else if (dc->add(dc, bm) < 0)
356				goto error;
357			if (empty < 0)
358				goto error;
359			level--;
360			init = 0;
361			continue;
362		}
363		if (init) {
364			int offset;
365			struct isl_tab_undo *snap2;
366			snap2 = isl_tab_snap(tab);
367			if (tab_add_divs(tab, map->p[level],
368					 &div_map[level]) < 0)
369				goto error;
370			offset = tab->n_con;
371			snap[level] = isl_tab_snap(tab);
372			if (tab_freeze_constraints(tab) < 0)
373				goto error;
374			if (tab_add_constraints(tab, map->p[level],
375						div_map[level]) < 0)
376				goto error;
377			k[level] = 0;
378			n[level] = 0;
379			if (tab->empty) {
380				if (isl_tab_rollback(tab, snap2) < 0)
381					goto error;
382				level++;
383				continue;
384			}
385			modified = 1;
386			n[level] = n_non_redundant(ctx, tab, offset,
387						    &index[level]);
388			if (n[level] < 0)
389				goto error;
390			if (n[level] == 0) {
391				level--;
392				init = 0;
393				continue;
394			}
395			if (isl_tab_rollback(tab, snap[level]) < 0)
396				goto error;
397			if (tab_add_constraint(tab, map->p[level],
398					div_map[level], index[level][0], 1) < 0)
399				goto error;
400			level++;
401			continue;
402		} else {
403			if (k[level] + 1 >= n[level]) {
404				level--;
405				continue;
406			}
407			if (isl_tab_rollback(tab, snap[level]) < 0)
408				goto error;
409			if (tab_add_constraint(tab, map->p[level],
410						div_map[level],
411						index[level][k[level]], 0) < 0)
412				goto error;
413			snap[level] = isl_tab_snap(tab);
414			k[level]++;
415			if (tab_add_constraint(tab, map->p[level],
416						div_map[level],
417						index[level][k[level]], 1) < 0)
418				goto error;
419			level++;
420			init = 1;
421			continue;
422		}
423	}
424
425	isl_tab_free(tab);
426	free(snap);
427	free(n);
428	free(k);
429	for (i = 0; index && i < map->n; ++i)
430		free(index[i]);
431	free(index);
432	for (i = 0; div_map && i < map->n; ++i)
433		free(div_map[i]);
434	free(div_map);
435
436	isl_basic_map_free(bmap);
437	isl_map_free(map);
438
439	return 0;
440error:
441	isl_tab_free(tab);
442	free(snap);
443	free(n);
444	free(k);
445	for (i = 0; index && i < map->n; ++i)
446		free(index[i]);
447	free(index);
448	for (i = 0; div_map && i < map->n; ++i)
449		free(div_map[i]);
450	free(div_map);
451	isl_basic_map_free(bmap);
452	isl_map_free(map);
453	return -1;
454}
455
456/* A diff collector that actually collects all parts of the
457 * set difference in the field diff.
458 */
459struct isl_subtract_diff_collector {
460	struct isl_diff_collector dc;
461	struct isl_map *diff;
462};
463
464/* isl_subtract_diff_collector callback.
465 */
466static int basic_map_subtract_add(struct isl_diff_collector *dc,
467			    __isl_take isl_basic_map *bmap)
468{
469	struct isl_subtract_diff_collector *sdc;
470	sdc = (struct isl_subtract_diff_collector *)dc;
471
472	sdc->diff = isl_map_union_disjoint(sdc->diff,
473			isl_map_from_basic_map(bmap));
474
475	return sdc->diff ? 0 : -1;
476}
477
478/* Return the set difference between bmap and map.
479 */
480static __isl_give isl_map *basic_map_subtract(__isl_take isl_basic_map *bmap,
481	__isl_take isl_map *map)
482{
483	struct isl_subtract_diff_collector sdc;
484	sdc.dc.add = &basic_map_subtract_add;
485	sdc.diff = isl_map_empty_like_basic_map(bmap);
486	if (basic_map_collect_diff(bmap, map, &sdc.dc) < 0) {
487		isl_map_free(sdc.diff);
488		sdc.diff = NULL;
489	}
490	return sdc.diff;
491}
492
493/* Return the set difference between map1 and map2.
494 * (U_i A_i) \ (U_j B_j) is computed as U_i (A_i \ (U_j B_j))
495 */
496static __isl_give isl_map *map_subtract( __isl_take isl_map *map1,
497	__isl_take isl_map *map2)
498{
499	int i;
500	struct isl_map *diff;
501
502	if (!map1 || !map2)
503		goto error;
504
505	isl_assert(map1->ctx, isl_space_is_equal(map1->dim, map2->dim), goto error);
506
507	if (isl_map_is_empty(map2)) {
508		isl_map_free(map2);
509		return map1;
510	}
511
512	map1 = isl_map_compute_divs(map1);
513	map2 = isl_map_compute_divs(map2);
514	if (!map1 || !map2)
515		goto error;
516
517	map1 = isl_map_remove_empty_parts(map1);
518	map2 = isl_map_remove_empty_parts(map2);
519
520	diff = isl_map_empty_like(map1);
521	for (i = 0; i < map1->n; ++i) {
522		struct isl_map *d;
523		d = basic_map_subtract(isl_basic_map_copy(map1->p[i]),
524				       isl_map_copy(map2));
525		if (ISL_F_ISSET(map1, ISL_MAP_DISJOINT))
526			diff = isl_map_union_disjoint(diff, d);
527		else
528			diff = isl_map_union(diff, d);
529	}
530
531	isl_map_free(map1);
532	isl_map_free(map2);
533
534	return diff;
535error:
536	isl_map_free(map1);
537	isl_map_free(map2);
538	return NULL;
539}
540
541__isl_give isl_map *isl_map_subtract( __isl_take isl_map *map1,
542	__isl_take isl_map *map2)
543{
544	return isl_map_align_params_map_map_and(map1, map2, &map_subtract);
545}
546
547struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
548{
549	return (struct isl_set *)
550		isl_map_subtract(
551			(struct isl_map *)set1, (struct isl_map *)set2);
552}
553
554/* Remove the elements of "dom" from the domain of "map".
555 */
556static __isl_give isl_map *map_subtract_domain(__isl_take isl_map *map,
557	__isl_take isl_set *dom)
558{
559	isl_map *ext_dom;
560
561	if (!isl_map_compatible_domain(map, dom))
562		isl_die(isl_set_get_ctx(dom), isl_error_invalid,
563			"incompatible spaces", goto error);
564
565	ext_dom = isl_map_universe(isl_map_get_space(map));
566	ext_dom = isl_map_intersect_domain(ext_dom, dom);
567	return isl_map_subtract(map, ext_dom);
568error:
569	isl_map_free(map);
570	isl_set_free(dom);
571	return NULL;
572}
573
574__isl_give isl_map *isl_map_subtract_domain(__isl_take isl_map *map,
575	__isl_take isl_set *dom)
576{
577	return isl_map_align_params_map_map_and(map, dom, &map_subtract_domain);
578}
579
580/* Remove the elements of "dom" from the range of "map".
581 */
582static __isl_give isl_map *map_subtract_range(__isl_take isl_map *map,
583	__isl_take isl_set *dom)
584{
585	isl_map *ext_dom;
586
587	if (!isl_map_compatible_range(map, dom))
588		isl_die(isl_set_get_ctx(dom), isl_error_invalid,
589			"incompatible spaces", goto error);
590
591	ext_dom = isl_map_universe(isl_map_get_space(map));
592	ext_dom = isl_map_intersect_range(ext_dom, dom);
593	return isl_map_subtract(map, ext_dom);
594error:
595	isl_map_free(map);
596	isl_set_free(dom);
597	return NULL;
598}
599
600__isl_give isl_map *isl_map_subtract_range(__isl_take isl_map *map,
601	__isl_take isl_set *dom)
602{
603	return isl_map_align_params_map_map_and(map, dom, &map_subtract_range);
604}
605
606/* A diff collector that aborts as soon as its add function is called,
607 * setting empty to 0.
608 */
609struct isl_is_empty_diff_collector {
610	struct isl_diff_collector dc;
611	int empty;
612};
613
614/* isl_is_empty_diff_collector callback.
615 */
616static int basic_map_is_empty_add(struct isl_diff_collector *dc,
617			    __isl_take isl_basic_map *bmap)
618{
619	struct isl_is_empty_diff_collector *edc;
620	edc = (struct isl_is_empty_diff_collector *)dc;
621
622	edc->empty = 0;
623
624	isl_basic_map_free(bmap);
625	return -1;
626}
627
628/* Check if bmap \ map is empty by computing this set difference
629 * and breaking off as soon as the difference is known to be non-empty.
630 */
631static int basic_map_diff_is_empty(__isl_keep isl_basic_map *bmap,
632	__isl_keep isl_map *map)
633{
634	int r;
635	struct isl_is_empty_diff_collector edc;
636
637	r = isl_basic_map_plain_is_empty(bmap);
638	if (r)
639		return r;
640
641	edc.dc.add = &basic_map_is_empty_add;
642	edc.empty = 1;
643	r = basic_map_collect_diff(isl_basic_map_copy(bmap),
644				   isl_map_copy(map), &edc.dc);
645	if (!edc.empty)
646		return 0;
647
648	return r < 0 ? -1 : 1;
649}
650
651/* Check if map1 \ map2 is empty by checking if the set difference is empty
652 * for each of the basic maps in map1.
653 */
654static int map_diff_is_empty(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
655{
656	int i;
657	int is_empty = 1;
658
659	if (!map1 || !map2)
660		return -1;
661
662	for (i = 0; i < map1->n; ++i) {
663		is_empty = basic_map_diff_is_empty(map1->p[i], map2);
664		if (is_empty < 0 || !is_empty)
665			 break;
666	}
667
668	return is_empty;
669}
670
671/* Return 1 if "bmap" contains a single element.
672 */
673int isl_basic_map_plain_is_singleton(__isl_keep isl_basic_map *bmap)
674{
675	if (!bmap)
676		return -1;
677	if (bmap->n_div)
678		return 0;
679	if (bmap->n_ineq)
680		return 0;
681	return bmap->n_eq == isl_basic_map_total_dim(bmap);
682}
683
684/* Return 1 if "map" contains a single element.
685 */
686int isl_map_plain_is_singleton(__isl_keep isl_map *map)
687{
688	if (!map)
689		return -1;
690	if (map->n != 1)
691		return 0;
692
693	return isl_basic_map_plain_is_singleton(map->p[0]);
694}
695
696/* Given a singleton basic map, extract the single element
697 * as an isl_point.
698 */
699static __isl_give isl_point *singleton_extract_point(
700	__isl_keep isl_basic_map *bmap)
701{
702	int j;
703	unsigned dim;
704	struct isl_vec *point;
705	isl_int m;
706
707	if (!bmap)
708		return NULL;
709
710	dim = isl_basic_map_total_dim(bmap);
711	isl_assert(bmap->ctx, bmap->n_eq == dim, return NULL);
712	point = isl_vec_alloc(bmap->ctx, 1 + dim);
713	if (!point)
714		return NULL;
715
716	isl_int_init(m);
717
718	isl_int_set_si(point->el[0], 1);
719	for (j = 0; j < bmap->n_eq; ++j) {
720		int i = dim - 1 - j;
721		isl_assert(bmap->ctx,
722		    isl_seq_first_non_zero(bmap->eq[j] + 1, i) == -1,
723		    goto error);
724		isl_assert(bmap->ctx,
725		    isl_int_is_one(bmap->eq[j][1 + i]) ||
726		    isl_int_is_negone(bmap->eq[j][1 + i]),
727		    goto error);
728		isl_assert(bmap->ctx,
729		    isl_seq_first_non_zero(bmap->eq[j]+1+i+1, dim-i-1) == -1,
730		    goto error);
731
732		isl_int_gcd(m, point->el[0], bmap->eq[j][1 + i]);
733		isl_int_divexact(m, bmap->eq[j][1 + i], m);
734		isl_int_abs(m, m);
735		isl_seq_scale(point->el, point->el, m, 1 + i);
736		isl_int_divexact(m, point->el[0], bmap->eq[j][1 + i]);
737		isl_int_neg(m, m);
738		isl_int_mul(point->el[1 + i], m, bmap->eq[j][0]);
739	}
740
741	isl_int_clear(m);
742	return isl_point_alloc(isl_basic_map_get_space(bmap), point);
743error:
744	isl_int_clear(m);
745	isl_vec_free(point);
746	return NULL;
747}
748
749/* Return 1 is the singleton map "map1" is a subset of "map2",
750 * i.e., if the single element of "map1" is also an element of "map2".
751 * Assumes "map2" has known divs.
752 */
753static int map_is_singleton_subset(__isl_keep isl_map *map1,
754	__isl_keep isl_map *map2)
755{
756	int i;
757	int is_subset = 0;
758	struct isl_point *point;
759
760	if (!map1 || !map2)
761		return -1;
762	if (map1->n != 1)
763		return -1;
764
765	point = singleton_extract_point(map1->p[0]);
766	if (!point)
767		return -1;
768
769	for (i = 0; i < map2->n; ++i) {
770		is_subset = isl_basic_map_contains_point(map2->p[i], point);
771		if (is_subset)
772			break;
773	}
774
775	isl_point_free(point);
776	return is_subset;
777}
778
779static int map_is_subset(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
780{
781	int is_subset = 0;
782	int empty;
783	int rat1, rat2;
784
785	if (!map1 || !map2)
786		return -1;
787
788	if (!isl_map_has_equal_space(map1, map2))
789		return 0;
790
791	empty = isl_map_is_empty(map1);
792	if (empty < 0)
793		return -1;
794	if (empty)
795		return 1;
796
797	empty = isl_map_is_empty(map2);
798	if (empty < 0)
799		return -1;
800	if (empty)
801		return 0;
802
803	rat1 = isl_map_has_rational(map1);
804	rat2 = isl_map_has_rational(map2);
805	if (rat1 < 0 || rat2 < 0)
806		return -1;
807	if (rat1 && !rat2)
808		return 0;
809
810	if (isl_map_plain_is_universe(map2))
811		return 1;
812
813	map2 = isl_map_compute_divs(isl_map_copy(map2));
814	if (isl_map_plain_is_singleton(map1)) {
815		is_subset = map_is_singleton_subset(map1, map2);
816		isl_map_free(map2);
817		return is_subset;
818	}
819	is_subset = map_diff_is_empty(map1, map2);
820	isl_map_free(map2);
821
822	return is_subset;
823}
824
825int isl_map_is_subset(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
826{
827	return isl_map_align_params_map_map_and_test(map1, map2,
828							&map_is_subset);
829}
830
831int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
832{
833	return isl_map_is_subset(
834			(struct isl_map *)set1, (struct isl_map *)set2);
835}
836
837__isl_give isl_map *isl_map_make_disjoint(__isl_take isl_map *map)
838{
839	int i;
840	struct isl_subtract_diff_collector sdc;
841	sdc.dc.add = &basic_map_subtract_add;
842
843	if (!map)
844		return NULL;
845	if (ISL_F_ISSET(map, ISL_MAP_DISJOINT))
846		return map;
847	if (map->n <= 1)
848		return map;
849
850	map = isl_map_compute_divs(map);
851	map = isl_map_remove_empty_parts(map);
852
853	if (!map || map->n <= 1)
854		return map;
855
856	sdc.diff = isl_map_from_basic_map(isl_basic_map_copy(map->p[0]));
857
858	for (i = 1; i < map->n; ++i) {
859		struct isl_basic_map *bmap = isl_basic_map_copy(map->p[i]);
860		struct isl_map *copy = isl_map_copy(sdc.diff);
861		if (basic_map_collect_diff(bmap, copy, &sdc.dc) < 0) {
862			isl_map_free(sdc.diff);
863			sdc.diff = NULL;
864			break;
865		}
866	}
867
868	isl_map_free(map);
869
870	return sdc.diff;
871}
872
873__isl_give isl_set *isl_set_make_disjoint(__isl_take isl_set *set)
874{
875	return (struct isl_set *)isl_map_make_disjoint((struct isl_map *)set);
876}
877
878__isl_give isl_map *isl_map_complement(__isl_take isl_map *map)
879{
880	isl_map *universe;
881
882	if (!map)
883		return NULL;
884
885	universe = isl_map_universe(isl_map_get_space(map));
886
887	return isl_map_subtract(universe, map);
888}
889
890__isl_give isl_set *isl_set_complement(__isl_take isl_set *set)
891{
892	return isl_map_complement(set);
893}
894