1	/*
2	* Copyright 2008-2009 Katholieke Universiteit Leuven
3	* Copyright 2012-2013 Ecole Normale Superieure
4	* Copyright 2014-2015 INRIA Rocquencourt
5	* Copyright 2016 Sven Verdoolaege
6	*
7	* Use of this software is governed by the MIT license
8	*
9	* Written by Sven Verdoolaege, K.U.Leuven, Departement
10	* Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11	* and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
12	* and Inria Paris - Rocquencourt, Domaine de Voluceau - Rocquencourt,
13	* B.P. 105 - 78153 Le Chesnay, France
14	*/
15
16	#include <isl_ctx_private.h>
17	#include <isl_map_private.h>
18	#include "isl_equalities.h"
19	#include <isl/map.h>
20	#include <isl_seq.h>
21	#include "isl_tab.h"
22	#include <isl_space_private.h>
23	#include <isl_mat_private.h>
24	#include <isl_vec_private.h>
25
26	#include <bset_to_bmap.c>
27	#include <bset_from_bmap.c>
28	#include <set_to_map.c>
29	#include <set_from_map.c>
30
31	static void swap_equality(__isl_keep isl_basic_map *bmap, int a, int b)
32	{
33	isl_int *t = bmap->eq[a];
34	bmap->eq[a] = bmap->eq[b];
35	bmap->eq[b] = t;
36	}
37
38	static void swap_inequality(__isl_keep isl_basic_map *bmap, int a, int b)
39	{
40	if (a != b) {
41	isl_int *t = bmap->ineq[a];
42	bmap->ineq[a] = bmap->ineq[b];
43	bmap->ineq[b] = t;
44	}
45	}
46
47	__isl_give isl_basic_map *isl_basic_map_normalize_constraints(
48	__isl_take isl_basic_map *bmap)
49	{
50	int i;
51	isl_int gcd;
52	isl_size total = isl_basic_map_dim(bmap, type: isl_dim_all);
53
54	if (total < 0)
55	return isl_basic_map_free(bmap);
56
57	isl_int_init(gcd);
58	for (i = bmap->n_eq - 1; i >= 0; --i) {
59	isl_seq_gcd(p: bmap->eq[i]+1, len: total, gcd: &gcd);
60	if (isl_int_is_zero(gcd)) {
61	if (!isl_int_is_zero(bmap->eq[i][0])) {
62	bmap = isl_basic_map_set_to_empty(bmap);
63	break;
64	}
65	if (isl_basic_map_drop_equality(bmap, pos: i) < 0)
66	goto error;
67	continue;
68	}
69	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
70	isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
71	if (isl_int_is_one(gcd))
72	continue;
73	if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
74	bmap = isl_basic_map_set_to_empty(bmap);
75	break;
76	}
77	isl_seq_scale_down(dst: bmap->eq[i], src: bmap->eq[i], f: gcd, len: 1+total);
78	}
79
80	for (i = bmap->n_ineq - 1; i >= 0; --i) {
81	isl_seq_gcd(p: bmap->ineq[i]+1, len: total, gcd: &gcd);
82	if (isl_int_is_zero(gcd)) {
83	if (isl_int_is_neg(bmap->ineq[i][0])) {
84	bmap = isl_basic_map_set_to_empty(bmap);
85	break;
86	}
87	if (isl_basic_map_drop_inequality(bmap, pos: i) < 0)
88	goto error;
89	continue;
90	}
91	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
92	isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
93	if (isl_int_is_one(gcd))
94	continue;
95	isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
96	isl_seq_scale_down(dst: bmap->ineq[i]+1, src: bmap->ineq[i]+1, f: gcd, len: total);
97	}
98	isl_int_clear(gcd);
99
100	return bmap;
101	error:
102	isl_int_clear(gcd);
103	isl_basic_map_free(bmap);
104	return NULL;
105	}
106
107	__isl_give isl_basic_set *isl_basic_set_normalize_constraints(
108	__isl_take isl_basic_set *bset)
109	{
110	isl_basic_map *bmap = bset_to_bmap(bset);
111	return bset_from_bmap(bmap: isl_basic_map_normalize_constraints(bmap));
112	}
113
114	/* Reduce the coefficient of the variable at position "pos"
115	* in integer division "div", such that it lies in the half-open
116	* interval (1/2,1/2], extracting any excess value from this integer division.
117	* "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
118	* corresponds to the constant term.
119	*
120	* That is, the integer division is of the form
121	*
122	* floor((... + (c * d + r) * x_pos + ...)/d)
123	*
124	* with -d < 2 * r <= d.
125	* Replace it by
126	*
127	* floor((... + r * x_pos + ...)/d) + c * x_pos
128	*
129	* If 2 * ((c * d + r) % d) <= d, then c = floor((c * d + r)/d).
130	* Otherwise, c = floor((c * d + r)/d) + 1.
131	*
132	* This is the same normalization that is performed by isl_aff_floor.
133	*/
134	static __isl_give isl_basic_map *reduce_coefficient_in_div(
135	__isl_take isl_basic_map *bmap, int div, int pos)
136	{
137	isl_int shift;
138	int add_one;
139
140	isl_int_init(shift);
141	isl_int_fdiv_r(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
142	isl_int_mul_ui(shift, shift, 2);
143	add_one = isl_int_gt(shift, bmap->div[div][0]);
144	isl_int_fdiv_q(shift, bmap->div[div][1 + pos], bmap->div[div][0]);
145	if (add_one)
146	isl_int_add_ui(shift, shift, 1);
147	isl_int_neg(shift, shift);
148	bmap = isl_basic_map_shift_div(bmap, div, pos, shift);
149	isl_int_clear(shift);
150
151	return bmap;
152	}
153
154	/* Does the coefficient of the variable at position "pos"
155	* in integer division "div" need to be reduced?
156	* That is, does it lie outside the half-open interval (1/2,1/2]?
157	* The coefficient c/d lies outside this interval if abs(2 * c) >= d and
158	* 2 * c != d.
159	*/
160	static isl_bool needs_reduction(__isl_keep isl_basic_map *bmap, int div,
161	int pos)
162	{
163	isl_bool r;
164
165	if (isl_int_is_zero(bmap->div[div][1 + pos]))
166	return isl_bool_false;
167
168	isl_int_mul_ui(bmap->div[div][1 + pos], bmap->div[div][1 + pos], 2);
169	r = isl_int_abs_ge(bmap->div[div][1 + pos], bmap->div[div][0]) &&
170	!isl_int_eq(bmap->div[div][1 + pos], bmap->div[div][0]);
171	isl_int_divexact_ui(bmap->div[div][1 + pos],
172	bmap->div[div][1 + pos], 2);
173
174	return r;
175	}
176
177	/* Reduce the coefficients (including the constant term) of
178	* integer division "div", if needed.
179	* In particular, make sure all coefficients lie in
180	* the half-open interval (1/2,1/2].
181	*/
182	static __isl_give isl_basic_map *reduce_div_coefficients_of_div(
183	__isl_take isl_basic_map *bmap, int div)
184	{
185	int i;
186	isl_size total;
187
188	total = isl_basic_map_dim(bmap, type: isl_dim_all);
189	if (total < 0)
190	return isl_basic_map_free(bmap);
191	for (i = 0; i < 1 + total; ++i) {
192	isl_bool reduce;
193
194	reduce = needs_reduction(bmap, div, pos: i);
195	if (reduce < 0)
196	return isl_basic_map_free(bmap);
197	if (!reduce)
198	continue;
199	bmap = reduce_coefficient_in_div(bmap, div, pos: i);
200	if (!bmap)
201	break;
202	}
203
204	return bmap;
205	}
206
207	/* Reduce the coefficients (including the constant term) of
208	* the known integer divisions, if needed
209	* In particular, make sure all coefficients lie in
210	* the half-open interval (1/2,1/2].
211	*/
212	static __isl_give isl_basic_map *reduce_div_coefficients(
213	__isl_take isl_basic_map *bmap)
214	{
215	int i;
216
217	if (!bmap)
218	return NULL;
219	if (bmap->n_div == 0)
220	return bmap;
221
222	for (i = 0; i < bmap->n_div; ++i) {
223	if (isl_int_is_zero(bmap->div[i][0]))
224	continue;
225	bmap = reduce_div_coefficients_of_div(bmap, div: i);
226	if (!bmap)
227	break;
228	}
229
230	return bmap;
231	}
232
233	/* Remove any common factor in numerator and denominator of the div expression,
234	* not taking into account the constant term.
235	* That is, if the div is of the form
236	*
237	* floor((a + m f(x))/(m d))
238	*
239	* then replace it by
240	*
241	* floor((floor(a/m) + f(x))/d)
242	*
243	* The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
244	* and can therefore not influence the result of the floor.
245	*/
246	static __isl_give isl_basic_map *normalize_div_expression(
247	__isl_take isl_basic_map *bmap, int div)
248	{
249	isl_size total = isl_basic_map_dim(bmap, type: isl_dim_all);
250	isl_ctx *ctx = bmap->ctx;
251
252	if (total < 0)
253	return isl_basic_map_free(bmap);
254	if (isl_int_is_zero(bmap->div[div][0]))
255	return bmap;
256	isl_seq_gcd(p: bmap->div[div] + 2, len: total, gcd: &ctx->normalize_gcd);
257	isl_int_gcd(ctx->normalize_gcd, ctx->normalize_gcd, bmap->div[div][0]);
258	if (isl_int_is_one(ctx->normalize_gcd))
259	return bmap;
260	isl_int_fdiv_q(bmap->div[div][1], bmap->div[div][1],
261	ctx->normalize_gcd);
262	isl_int_divexact(bmap->div[div][0], bmap->div[div][0],
263	ctx->normalize_gcd);
264	isl_seq_scale_down(dst: bmap->div[div] + 2, src: bmap->div[div] + 2,
265	f: ctx->normalize_gcd, len: total);
266
267	return bmap;
268	}
269
270	/* Remove any common factor in numerator and denominator of a div expression,
271	* not taking into account the constant term.
272	* That is, look for any div of the form
273	*
274	* floor((a + m f(x))/(m d))
275	*
276	* and replace it by
277	*
278	* floor((floor(a/m) + f(x))/d)
279	*
280	* The difference {a/m}/d in the argument satisfies 0 <= {a/m}/d < 1/d
281	* and can therefore not influence the result of the floor.
282	*/
283	static __isl_give isl_basic_map *normalize_div_expressions(
284	__isl_take isl_basic_map *bmap)
285	{
286	int i;
287
288	if (!bmap)
289	return NULL;
290	if (bmap->n_div == 0)
291	return bmap;
292
293	for (i = 0; i < bmap->n_div; ++i)
294	bmap = normalize_div_expression(bmap, div: i);
295
296	return bmap;
297	}
298
299	/* Assumes divs have been ordered if keep_divs is set.
300	*/
301	static __isl_give isl_basic_map *eliminate_var_using_equality(
302	__isl_take isl_basic_map *bmap,
303	unsigned pos, isl_int *eq, int keep_divs, int *progress)
304	{
305	isl_size total;
306	isl_size v_div;
307	int k;
308	int last_div;
309
310	total = isl_basic_map_dim(bmap, type: isl_dim_all);
311	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
312	if (total < 0 || v_div < 0)
313	return isl_basic_map_free(bmap);
314	last_div = isl_seq_last_non_zero(p: eq + 1 + v_div, len: bmap->n_div);
315	for (k = 0; k < bmap->n_eq; ++k) {
316	if (bmap->eq[k] == eq)
317	continue;
318	if (isl_int_is_zero(bmap->eq[k][1+pos]))
319	continue;
320	if (progress)
321	*progress = 1;
322	isl_seq_elim(dst: bmap->eq[k], src: eq, pos: 1+pos, len: 1+total, NULL);
323	isl_seq_normalize(ctx: bmap->ctx, p: bmap->eq[k], len: 1 + total);
324	}
325
326	for (k = 0; k < bmap->n_ineq; ++k) {
327	if (isl_int_is_zero(bmap->ineq[k][1+pos]))
328	continue;
329	if (progress)
330	*progress = 1;
331	isl_seq_elim(dst: bmap->ineq[k], src: eq, pos: 1+pos, len: 1+total, NULL);
332	isl_seq_normalize(ctx: bmap->ctx, p: bmap->ineq[k], len: 1 + total);
333	ISL_F_CLR(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
334	ISL_F_CLR(bmap, ISL_BASIC_MAP_SORTED);
335	}
336
337	for (k = 0; k < bmap->n_div; ++k) {
338	if (isl_int_is_zero(bmap->div[k][0]))
339	continue;
340	if (isl_int_is_zero(bmap->div[k][1+1+pos]))
341	continue;
342	if (progress)
343	*progress = 1;
344	/* We need to be careful about circular definitions,
345	* so for now we just remove the definition of div k
346	* if the equality contains any divs.
347	* If keep_divs is set, then the divs have been ordered
348	* and we can keep the definition as long as the result
349	* is still ordered.
350	*/
351	if (last_div == -1 || (keep_divs && last_div < k)) {
352	isl_seq_elim(dst: bmap->div[k]+1, src: eq,
353	pos: 1+pos, len: 1+total, m: &bmap->div[k][0]);
354	bmap = normalize_div_expression(bmap, div: k);
355	if (!bmap)
356	return NULL;
357	} else
358	isl_seq_clr(p: bmap->div[k], len: 1 + total);
359	}
360
361	return bmap;
362	}
363
364	/* Assumes divs have been ordered if keep_divs is set.
365	*/
366	static __isl_give isl_basic_map *eliminate_div(__isl_take isl_basic_map *bmap,
367	isl_int *eq, unsigned div, int keep_divs)
368	{
369	isl_size v_div;
370	unsigned pos;
371
372	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
373	if (v_div < 0)
374	return isl_basic_map_free(bmap);
375	pos = v_div + div;
376	bmap = eliminate_var_using_equality(bmap, pos, eq, keep_divs, NULL);
377
378	bmap = isl_basic_map_drop_div(bmap, div);
379
380	return bmap;
381	}
382
383	/* Check if elimination of div "div" using equality "eq" would not
384	* result in a div depending on a later div.
385	*/
386	static isl_bool ok_to_eliminate_div(__isl_keep isl_basic_map *bmap, isl_int *eq,
387	unsigned div)
388	{
389	int k;
390	int last_div;
391	isl_size v_div;
392	unsigned pos;
393
394	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
395	if (v_div < 0)
396	return isl_bool_error;
397	pos = v_div + div;
398
399	last_div = isl_seq_last_non_zero(p: eq + 1 + v_div, len: bmap->n_div);
400	if (last_div < 0 || last_div <= div)
401	return isl_bool_true;
402
403	for (k = 0; k <= last_div; ++k) {
404	if (isl_int_is_zero(bmap->div[k][0]))
405	continue;
406	if (!isl_int_is_zero(bmap->div[k][1 + 1 + pos]))
407	return isl_bool_false;
408	}
409
410	return isl_bool_true;
411	}
412
413	/* Eliminate divs based on equalities
414	*/
415	static __isl_give isl_basic_map *eliminate_divs_eq(
416	__isl_take isl_basic_map *bmap, int *progress)
417	{
418	int d;
419	int i;
420	int modified = 0;
421	unsigned off;
422
423	bmap = isl_basic_map_order_divs(bmap);
424
425	if (!bmap)
426	return NULL;
427
428	off = isl_basic_map_offset(bmap, type: isl_dim_div);
429
430	for (d = bmap->n_div - 1; d >= 0 ; --d) {
431	for (i = 0; i < bmap->n_eq; ++i) {
432	isl_bool ok;
433
434	if (!isl_int_is_one(bmap->eq[i][off + d]) &&
435	!isl_int_is_negone(bmap->eq[i][off + d]))
436	continue;
437	ok = ok_to_eliminate_div(bmap, eq: bmap->eq[i], div: d);
438	if (ok < 0)
439	return isl_basic_map_free(bmap);
440	if (!ok)
441	continue;
442	modified = 1;
443	*progress = 1;
444	bmap = eliminate_div(bmap, eq: bmap->eq[i], div: d, keep_divs: 1);
445	if (isl_basic_map_drop_equality(bmap, pos: i) < 0)
446	return isl_basic_map_free(bmap);
447	break;
448	}
449	}
450	if (modified)
451	return eliminate_divs_eq(bmap, progress);
452	return bmap;
453	}
454
455	/* Eliminate divs based on inequalities
456	*/
457	static __isl_give isl_basic_map *eliminate_divs_ineq(
458	__isl_take isl_basic_map *bmap, int *progress)
459	{
460	int d;
461	int i;
462	unsigned off;
463	struct isl_ctx *ctx;
464
465	if (!bmap)
466	return NULL;
467
468	ctx = bmap->ctx;
469	off = isl_basic_map_offset(bmap, type: isl_dim_div);
470
471	for (d = bmap->n_div - 1; d >= 0 ; --d) {
472	for (i = 0; i < bmap->n_eq; ++i)
473	if (!isl_int_is_zero(bmap->eq[i][off + d]))
474	break;
475	if (i < bmap->n_eq)
476	continue;
477	for (i = 0; i < bmap->n_ineq; ++i)
478	if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
479	break;
480	if (i < bmap->n_ineq)
481	continue;
482	*progress = 1;
483	bmap = isl_basic_map_eliminate_vars(bmap, pos: (off-1)+d, n: 1);
484	if (!bmap || ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
485	break;
486	bmap = isl_basic_map_drop_div(bmap, div: d);
487	if (!bmap)
488	break;
489	}
490	return bmap;
491	}
492
493	/* Does the equality constraint at position "eq" in "bmap" involve
494	* any local variables in the range [first, first + n)
495	* that are not marked as having an explicit representation?
496	*/
497	static isl_bool bmap_eq_involves_unknown_divs(__isl_keep isl_basic_map *bmap,
498	int eq, unsigned first, unsigned n)
499	{
500	unsigned o_div;
501	int i;
502
503	if (!bmap)
504	return isl_bool_error;
505
506	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
507	for (i = 0; i < n; ++i) {
508	isl_bool unknown;
509
510	if (isl_int_is_zero(bmap->eq[eq][o_div + first + i]))
511	continue;
512	unknown = isl_basic_map_div_is_marked_unknown(bmap, div: first + i);
513	if (unknown < 0)
514	return isl_bool_error;
515	if (unknown)
516	return isl_bool_true;
517	}
518
519	return isl_bool_false;
520	}
521
522	/* The last local variable involved in the equality constraint
523	* at position "eq" in "bmap" is the local variable at position "div".
524	* It can therefore be used to extract an explicit representation
525	* for that variable.
526	* Do so unless the local variable already has an explicit representation or
527	* the explicit representation would involve any other local variables
528	* that in turn do not have an explicit representation.
529	* An equality constraint involving local variables without an explicit
530	* representation can be used in isl_basic_map_drop_redundant_divs
531	* to separate out an independent local variable. Introducing
532	* an explicit representation here would block this transformation,
533	* while the partial explicit representation in itself is not very useful.
534	* Set *progress if anything is changed.
535	*
536	* The equality constraint is of the form
537	*
538	* f(x) + n e >= 0
539	*
540	* with n a positive number. The explicit representation derived from
541	* this constraint is
542	*
543	* floor((-f(x))/n)
544	*/
545	static __isl_give isl_basic_map *set_div_from_eq(__isl_take isl_basic_map *bmap,
546	int div, int eq, int *progress)
547	{
548	isl_size total;
549	unsigned o_div;
550	isl_bool involves;
551
552	if (!bmap)
553	return NULL;
554
555	if (!isl_int_is_zero(bmap->div[div][0]))
556	return bmap;
557
558	involves = bmap_eq_involves_unknown_divs(bmap, eq, first: 0, n: div);
559	if (involves < 0)
560	return isl_basic_map_free(bmap);
561	if (involves)
562	return bmap;
563
564	total = isl_basic_map_dim(bmap, type: isl_dim_all);
565	if (total < 0)
566	return isl_basic_map_free(bmap);
567	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
568	isl_seq_neg(dst: bmap->div[div] + 1, src: bmap->eq[eq], len: 1 + total);
569	isl_int_set_si(bmap->div[div][1 + o_div + div], 0);
570	isl_int_set(bmap->div[div][0], bmap->eq[eq][o_div + div]);
571	if (progress)
572	*progress = 1;
573
574	return bmap;
575	}
576
577	/* Perform fangcheng (Gaussian elimination) on the equality
578	* constraints of "bmap".
579	* That is, put them into row-echelon form, starting from the last column
580	* backward and use them to eliminate the corresponding coefficients
581	* from all constraints.
582	*
583	* If "progress" is not NULL, then it gets set if the elimination
584	* results in any changes.
585	* The elimination process may result in some equality constraints
586	* getting interchanged or removed.
587	* If "swap" or "drop" are not NULL, then they get called when
588	* two equality constraints get interchanged or
589	* when a number of final equality constraints get removed.
590	* As a special case, if the input turns out to be empty,
591	* then drop gets called with the number of removed equality
592	* constraints set to the total number of equality constraints.
593	* If "swap" or "drop" are not NULL, then the local variables (if any)
594	* are assumed to be in a valid order.
595	*/
596	__isl_give isl_basic_map *isl_basic_map_gauss5(__isl_take isl_basic_map *bmap,
597	int *progress,
598	isl_stat (*swap)(unsigned a, unsigned b, void *user),
599	isl_stat (*drop)(unsigned n, void *user), void *user)
600	{
601	int k;
602	int done;
603	int last_var;
604	unsigned total_var;
605	isl_size total;
606	unsigned n_drop;
607
608	if (!swap && !drop)
609	bmap = isl_basic_map_order_divs(bmap);
610
611	total = isl_basic_map_dim(bmap, type: isl_dim_all);
612	if (total < 0)
613	return isl_basic_map_free(bmap);
614
615	total_var = total - bmap->n_div;
616
617	last_var = total - 1;
618	for (done = 0; done < bmap->n_eq; ++done) {
619	for (; last_var >= 0; --last_var) {
620	for (k = done; k < bmap->n_eq; ++k)
621	if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
622	break;
623	if (k < bmap->n_eq)
624	break;
625	}
626	if (last_var < 0)
627	break;
628	if (k != done) {
629	swap_equality(bmap, a: k, b: done);
630	if (swap && swap(k, done, user) < 0)
631	return isl_basic_map_free(bmap);
632	}
633	if (isl_int_is_neg(bmap->eq[done][1+last_var]))
634	isl_seq_neg(dst: bmap->eq[done], src: bmap->eq[done], len: 1+total);
635
636	bmap = eliminate_var_using_equality(bmap, pos: last_var,
637	eq: bmap->eq[done], keep_divs: 1, progress);
638
639	if (last_var >= total_var)
640	bmap = set_div_from_eq(bmap, div: last_var - total_var,
641	eq: done, progress);
642	if (!bmap)
643	return NULL;
644	}
645	if (done == bmap->n_eq)
646	return bmap;
647	for (k = done; k < bmap->n_eq; ++k) {
648	if (isl_int_is_zero(bmap->eq[k][0]))
649	continue;
650	if (drop && drop(bmap->n_eq, user) < 0)
651	return isl_basic_map_free(bmap);
652	return isl_basic_map_set_to_empty(bmap);
653	}
654	n_drop = bmap->n_eq - done;
655	bmap = isl_basic_map_free_equality(bmap, n: n_drop);
656	if (drop && drop(n_drop, user) < 0)
657	return isl_basic_map_free(bmap);
658	return bmap;
659	}
660
661	__isl_give isl_basic_map *isl_basic_map_gauss(__isl_take isl_basic_map *bmap,
662	int *progress)
663	{
664	return isl_basic_map_gauss5(bmap, progress, NULL, NULL, NULL);
665	}
666
667	__isl_give isl_basic_set *isl_basic_set_gauss(
668	__isl_take isl_basic_set *bset, int *progress)
669	{
670	return bset_from_bmap(bmap: isl_basic_map_gauss(bmap: bset_to_bmap(bset),
671	progress));
672	}
673
674
675	static unsigned int round_up(unsigned int v)
676	{
677	int old_v = v;
678
679	while (v) {
680	old_v = v;
681	v ^= v & -v;
682	}
683	return old_v << 1;
684	}
685
686	/* Hash table of inequalities in a basic map.
687	* "index" is an array of addresses of inequalities in the basic map, some
688	* of which are NULL. The inequalities are hashed on the coefficients
689	* except the constant term.
690	* "size" is the number of elements in the array and is always a power of two
691	* "bits" is the number of bits need to represent an index into the array.
692	* "total" is the total dimension of the basic map.
693	*/
694	struct isl_constraint_index {
695	unsigned int size;
696	int bits;
697	isl_int ***index;
698	isl_size total;
699	};
700
701	/* Fill in the "ci" data structure for holding the inequalities of "bmap".
702	*/
703	static isl_stat create_constraint_index(struct isl_constraint_index *ci,
704	__isl_keep isl_basic_map *bmap)
705	{
706	isl_ctx *ctx;
707
708	ci->index = NULL;
709	if (!bmap)
710	return isl_stat_error;
711	ci->total = isl_basic_map_dim(bmap, type: isl_dim_all);
712	if (ci->total < 0)
713	return isl_stat_error;
714	if (bmap->n_ineq == 0)
715	return isl_stat_ok;
716	ci->size = round_up(v: 4 * (bmap->n_ineq + 1) / 3 - 1);
717	ci->bits = ffs(i: ci->size) - 1;
718	ctx = isl_basic_map_get_ctx(bmap);
719	ci->index = isl_calloc_array(ctx, isl_int **, ci->size);
720	if (!ci->index)
721	return isl_stat_error;
722
723	return isl_stat_ok;
724	}
725
726	/* Free the memory allocated by create_constraint_index.
727	*/
728	static void constraint_index_free(struct isl_constraint_index *ci)
729	{
730	free(ptr: ci->index);
731	}
732
733	/* Return the position in ci->index that contains the address of
734	* an inequality that is equal to *ineq up to the constant term,
735	* provided this address is not identical to "ineq".
736	* If there is no such inequality, then return the position where
737	* such an inequality should be inserted.
738	*/
739	static int hash_index_ineq(struct isl_constraint_index *ci, isl_int **ineq)
740	{
741	int h;
742	uint32_t hash = isl_seq_get_hash_bits(p: (*ineq) + 1, len: ci->total, bits: ci->bits);
743	for (h = hash; ci->index[h]; h = (h+1) % ci->size)
744	if (ineq != ci->index[h] &&
745	isl_seq_eq(p1: (*ineq) + 1, p2: ci->index[h][0]+1, len: ci->total))
746	break;
747	return h;
748	}
749
750	/* Return the position in ci->index that contains the address of
751	* an inequality that is equal to the k'th inequality of "bmap"
752	* up to the constant term, provided it does not point to the very
753	* same inequality.
754	* If there is no such inequality, then return the position where
755	* such an inequality should be inserted.
756	*/
757	static int hash_index(struct isl_constraint_index *ci,
758	__isl_keep isl_basic_map *bmap, int k)
759	{
760	return hash_index_ineq(ci, ineq: &bmap->ineq[k]);
761	}
762
763	static int set_hash_index(struct isl_constraint_index *ci,
764	__isl_keep isl_basic_set *bset, int k)
765	{
766	return hash_index(ci, bmap: bset, k);
767	}
768
769	/* Fill in the "ci" data structure with the inequalities of "bset".
770	*/
771	static isl_stat setup_constraint_index(struct isl_constraint_index *ci,
772	__isl_keep isl_basic_set *bset)
773	{
774	int k, h;
775
776	if (create_constraint_index(ci, bmap: bset) < 0)
777	return isl_stat_error;
778
779	for (k = 0; k < bset->n_ineq; ++k) {
780	h = set_hash_index(ci, bset, k);
781	ci->index[h] = &bset->ineq[k];
782	}
783
784	return isl_stat_ok;
785	}
786
787	/* Is the inequality ineq (obviously) redundant with respect
788	* to the constraints in "ci"?
789	*
790	* Look for an inequality in "ci" with the same coefficients and then
791	* check if the contant term of "ineq" is greater than or equal
792	* to the constant term of that inequality. If so, "ineq" is clearly
793	* redundant.
794	*
795	* Note that hash_index_ineq ignores a stored constraint if it has
796	* the same address as the passed inequality. It is ok to pass
797	* the address of a local variable here since it will never be
798	* the same as the address of a constraint in "ci".
799	*/
800	static isl_bool constraint_index_is_redundant(struct isl_constraint_index *ci,
801	isl_int *ineq)
802	{
803	int h;
804
805	h = hash_index_ineq(ci, ineq: &ineq);
806	if (!ci->index[h])
807	return isl_bool_false;
808	return isl_int_ge(ineq[0], (*ci->index[h])[0]);
809	}
810
811	/* If we can eliminate more than one div, then we need to make
812	* sure we do it from last div to first div, in order not to
813	* change the position of the other divs that still need to
814	* be removed.
815	*/
816	static __isl_give isl_basic_map *remove_duplicate_divs(
817	__isl_take isl_basic_map *bmap, int *progress)
818	{
819	unsigned int size;
820	int *index;
821	int *elim_for;
822	int k, l, h;
823	int bits;
824	struct isl_blk eq;
825	isl_size v_div;
826	unsigned total;
827	struct isl_ctx *ctx;
828
829	bmap = isl_basic_map_order_divs(bmap);
830	if (!bmap || bmap->n_div <= 1)
831	return bmap;
832
833	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
834	if (v_div < 0)
835	return isl_basic_map_free(bmap);
836	total = v_div + bmap->n_div;
837
838	ctx = bmap->ctx;
839	for (k = bmap->n_div - 1; k >= 0; --k)
840	if (!isl_int_is_zero(bmap->div[k][0]))
841	break;
842	if (k <= 0)
843	return bmap;
844
845	size = round_up(v: 4 * bmap->n_div / 3 - 1);
846	if (size == 0)
847	return bmap;
848	elim_for = isl_calloc_array(ctx, int, bmap->n_div);
849	bits = ffs(i: size) - 1;
850	index = isl_calloc_array(ctx, int, size);
851	if (!elim_for || !index)
852	goto out;
853	eq = isl_blk_alloc(ctx, n: 1+total);
854	if (isl_blk_is_error(block: eq))
855	goto out;
856
857	isl_seq_clr(p: eq.data, len: 1+total);
858	index[isl_seq_get_hash_bits(p: bmap->div[k], len: 2+total, bits)] = k + 1;
859	for (--k; k >= 0; --k) {
860	uint32_t hash;
861
862	if (isl_int_is_zero(bmap->div[k][0]))
863	continue;
864
865	hash = isl_seq_get_hash_bits(p: bmap->div[k], len: 2+total, bits);
866	for (h = hash; index[h]; h = (h+1) % size)
867	if (isl_seq_eq(p1: bmap->div[k],
868	p2: bmap->div[index[h]-1], len: 2+total))
869	break;
870	if (index[h]) {
871	*progress = 1;
872	l = index[h] - 1;
873	elim_for[l] = k + 1;
874	}
875	index[h] = k+1;
876	}
877	for (l = bmap->n_div - 1; l >= 0; --l) {
878	if (!elim_for[l])
879	continue;
880	k = elim_for[l] - 1;
881	isl_int_set_si(eq.data[1 + v_div + k], -1);
882	isl_int_set_si(eq.data[1 + v_div + l], 1);
883	bmap = eliminate_div(bmap, eq: eq.data, div: l, keep_divs: 1);
884	if (!bmap)
885	break;
886	isl_int_set_si(eq.data[1 + v_div + k], 0);
887	isl_int_set_si(eq.data[1 + v_div + l], 0);
888	}
889
890	isl_blk_free(ctx, block: eq);
891	out:
892	free(ptr: index);
893	free(ptr: elim_for);
894	return bmap;
895	}
896
897	static int n_pure_div_eq(__isl_keep isl_basic_map *bmap)
898	{
899	int i, j;
900	isl_size v_div;
901
902	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
903	if (v_div < 0)
904	return -1;
905	for (i = 0, j = bmap->n_div-1; i < bmap->n_eq; ++i) {
906	while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + v_div + j]))
907	--j;
908	if (j < 0)
909	break;
910	if (isl_seq_first_non_zero(p: bmap->eq[i] + 1 + v_div, len: j) != -1)
911	return 0;
912	}
913	return i;
914	}
915
916	/* Normalize divs that appear in equalities.
917	*
918	* In particular, we assume that bmap contains some equalities
919	* of the form
920	*
921	* a x = m * e_i
922	*
923	* and we want to replace the set of e_i by a minimal set and
924	* such that the new e_i have a canonical representation in terms
925	* of the vector x.
926	* If any of the equalities involves more than one divs, then
927	* we currently simply bail out.
928	*
929	* Let us first additionally assume that all equalities involve
930	* a div. The equalities then express modulo constraints on the
931	* remaining variables and we can use "parameter compression"
932	* to find a minimal set of constraints. The result is a transformation
933	*
934	* x = T(x') = x_0 + G x'
935	*
936	* with G a lower-triangular matrix with all elements below the diagonal
937	* non-negative and smaller than the diagonal element on the same row.
938	* We first normalize x_0 by making the same property hold in the affine
939	* T matrix.
940	* The rows i of G with a 1 on the diagonal do not impose any modulo
941	* constraint and simply express x_i = x'_i.
942	* For each of the remaining rows i, we introduce a div and a corresponding
943	* equality. In particular
944	*
945	* g_ii e_j = x_i - g_i(x')
946	*
947	* where each x'_k is replaced either by x_k (if g_kk = 1) or the
948	* corresponding div (if g_kk != 1).
949	*
950	* If there are any equalities not involving any div, then we
951	* first apply a variable compression on the variables x:
952	*
953	* x = C x'' x'' = C_2 x
954	*
955	* and perform the above parameter compression on A C instead of on A.
956	* The resulting compression is then of the form
957	*
958	* x'' = T(x') = x_0 + G x'
959	*
960	* and in constructing the new divs and the corresponding equalities,
961	* we have to replace each x'', i.e., the x'_k with (g_kk = 1),
962	* by the corresponding row from C_2.
963	*/
964	static __isl_give isl_basic_map *normalize_divs(__isl_take isl_basic_map *bmap,
965	int *progress)
966	{
967	int i, j, k;
968	isl_size v_div;
969	int div_eq;
970	struct isl_mat *B;
971	struct isl_vec *d;
972	struct isl_mat *T = NULL;
973	struct isl_mat *C = NULL;
974	struct isl_mat *C2 = NULL;
975	isl_int v;
976	int *pos = NULL;
977	int dropped, needed;
978
979	if (!bmap)
980	return NULL;
981
982	if (bmap->n_div == 0)
983	return bmap;
984
985	if (bmap->n_eq == 0)
986	return bmap;
987
988	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS))
989	return bmap;
990
991	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
992	div_eq = n_pure_div_eq(bmap);
993	if (v_div < 0 || div_eq < 0)
994	return isl_basic_map_free(bmap);
995	if (div_eq == 0)
996	return bmap;
997
998	if (div_eq < bmap->n_eq) {
999	B = isl_mat_sub_alloc6(ctx: bmap->ctx, row: bmap->eq, first_row: div_eq,
1000	n_row: bmap->n_eq - div_eq, first_col: 0, n_col: 1 + v_div);
1001	C = isl_mat_variable_compression(B, T2: &C2);
1002	if (!C || !C2)
1003	goto error;
1004	if (C->n_col == 0) {
1005	bmap = isl_basic_map_set_to_empty(bmap);
1006	isl_mat_free(mat: C);
1007	isl_mat_free(mat: C2);
1008	goto done;
1009	}
1010	}
1011
1012	d = isl_vec_alloc(ctx: bmap->ctx, size: div_eq);
1013	if (!d)
1014	goto error;
1015	for (i = 0, j = bmap->n_div-1; i < div_eq; ++i) {
1016	while (j >= 0 && isl_int_is_zero(bmap->eq[i][1 + v_div + j]))
1017	--j;
1018	isl_int_set(d->block.data[i], bmap->eq[i][1 + v_div + j]);
1019	}
1020	B = isl_mat_sub_alloc6(ctx: bmap->ctx, row: bmap->eq, first_row: 0, n_row: div_eq, first_col: 0, n_col: 1 + v_div);
1021
1022	if (C) {
1023	B = isl_mat_product(left: B, right: C);
1024	C = NULL;
1025	}
1026
1027	T = isl_mat_parameter_compression(B, d);
1028	if (!T)
1029	goto error;
1030	if (T->n_col == 0) {
1031	bmap = isl_basic_map_set_to_empty(bmap);
1032	isl_mat_free(mat: C2);
1033	isl_mat_free(mat: T);
1034	goto done;
1035	}
1036	isl_int_init(v);
1037	for (i = 0; i < T->n_row - 1; ++i) {
1038	isl_int_fdiv_q(v, T->row[1 + i][0], T->row[1 + i][1 + i]);
1039	if (isl_int_is_zero(v))
1040	continue;
1041	isl_mat_col_submul(mat: T, dst_col: 0, f: v, src_col: 1 + i);
1042	}
1043	isl_int_clear(v);
1044	pos = isl_alloc_array(bmap->ctx, int, T->n_row);
1045	if (!pos)
1046	goto error;
1047	/* We have to be careful because dropping equalities may reorder them */
1048	dropped = 0;
1049	for (j = bmap->n_div - 1; j >= 0; --j) {
1050	for (i = 0; i < bmap->n_eq; ++i)
1051	if (!isl_int_is_zero(bmap->eq[i][1 + v_div + j]))
1052	break;
1053	if (i < bmap->n_eq) {
1054	bmap = isl_basic_map_drop_div(bmap, div: j);
1055	if (isl_basic_map_drop_equality(bmap, pos: i) < 0)
1056	goto error;
1057	++dropped;
1058	}
1059	}
1060	pos[0] = 0;
1061	needed = 0;
1062	for (i = 1; i < T->n_row; ++i) {
1063	if (isl_int_is_one(T->row[i][i]))
1064	pos[i] = i;
1065	else
1066	needed++;
1067	}
1068	if (needed > dropped) {
1069	bmap = isl_basic_map_extend(base: bmap, extra: needed, n_eq: needed, n_ineq: 0);
1070	if (!bmap)
1071	goto error;
1072	}
1073	for (i = 1; i < T->n_row; ++i) {
1074	if (isl_int_is_one(T->row[i][i]))
1075	continue;
1076	k = isl_basic_map_alloc_div(bmap);
1077	pos[i] = 1 + v_div + k;
1078	isl_seq_clr(p: bmap->div[k] + 1, len: 1 + v_div + bmap->n_div);
1079	isl_int_set(bmap->div[k][0], T->row[i][i]);
1080	if (C2)
1081	isl_seq_cpy(dst: bmap->div[k] + 1, src: C2->row[i], len: 1 + v_div);
1082	else
1083	isl_int_set_si(bmap->div[k][1 + i], 1);
1084	for (j = 0; j < i; ++j) {
1085	if (isl_int_is_zero(T->row[i][j]))
1086	continue;
1087	if (pos[j] < T->n_row && C2)
1088	isl_seq_submul(dst: bmap->div[k] + 1, f: T->row[i][j],
1089	src: C2->row[pos[j]], len: 1 + v_div);
1090	else
1091	isl_int_neg(bmap->div[k][1 + pos[j]],
1092	T->row[i][j]);
1093	}
1094	j = isl_basic_map_alloc_equality(bmap);
1095	isl_seq_neg(dst: bmap->eq[j], src: bmap->div[k]+1, len: 1+v_div+bmap->n_div);
1096	isl_int_set(bmap->eq[j][pos[i]], bmap->div[k][0]);
1097	}
1098	free(ptr: pos);
1099	isl_mat_free(mat: C2);
1100	isl_mat_free(mat: T);
1101
1102	if (progress)
1103	*progress = 1;
1104	done:
1105	ISL_F_SET(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
1106
1107	return bmap;
1108	error:
1109	free(ptr: pos);
1110	isl_mat_free(mat: C);
1111	isl_mat_free(mat: C2);
1112	isl_mat_free(mat: T);
1113	isl_basic_map_free(bmap);
1114	return NULL;
1115	}
1116
1117	static __isl_give isl_basic_map *set_div_from_lower_bound(
1118	__isl_take isl_basic_map *bmap, int div, int ineq)
1119	{
1120	unsigned total = isl_basic_map_offset(bmap, type: isl_dim_div);
1121
1122	isl_seq_neg(dst: bmap->div[div] + 1, src: bmap->ineq[ineq], len: total + bmap->n_div);
1123	isl_int_set(bmap->div[div][0], bmap->ineq[ineq][total + div]);
1124	isl_int_add(bmap->div[div][1], bmap->div[div][1], bmap->div[div][0]);
1125	isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1126	isl_int_set_si(bmap->div[div][1 + total + div], 0);
1127
1128	return bmap;
1129	}
1130
1131	/* Check whether it is ok to define a div based on an inequality.
1132	* To avoid the introduction of circular definitions of divs, we
1133	* do not allow such a definition if the resulting expression would refer to
1134	* any other undefined divs or if any known div is defined in
1135	* terms of the unknown div.
1136	*/
1137	static isl_bool ok_to_set_div_from_bound(__isl_keep isl_basic_map *bmap,
1138	int div, int ineq)
1139	{
1140	int j;
1141	unsigned total = isl_basic_map_offset(bmap, type: isl_dim_div);
1142
1143	/* Not defined in terms of unknown divs */
1144	for (j = 0; j < bmap->n_div; ++j) {
1145	if (div == j)
1146	continue;
1147	if (isl_int_is_zero(bmap->ineq[ineq][total + j]))
1148	continue;
1149	if (isl_int_is_zero(bmap->div[j][0]))
1150	return isl_bool_false;
1151	}
1152
1153	/* No other div defined in terms of this one => avoid loops */
1154	for (j = 0; j < bmap->n_div; ++j) {
1155	if (div == j)
1156	continue;
1157	if (isl_int_is_zero(bmap->div[j][0]))
1158	continue;
1159	if (!isl_int_is_zero(bmap->div[j][1 + total + div]))
1160	return isl_bool_false;
1161	}
1162
1163	return isl_bool_true;
1164	}
1165
1166	/* Would an expression for div "div" based on inequality "ineq" of "bmap"
1167	* be a better expression than the current one?
1168	*
1169	* If we do not have any expression yet, then any expression would be better.
1170	* Otherwise we check if the last variable involved in the inequality
1171	* (disregarding the div that it would define) is in an earlier position
1172	* than the last variable involved in the current div expression.
1173	*/
1174	static isl_bool better_div_constraint(__isl_keep isl_basic_map *bmap,
1175	int div, int ineq)
1176	{
1177	unsigned total = isl_basic_map_offset(bmap, type: isl_dim_div);
1178	int last_div;
1179	int last_ineq;
1180
1181	if (isl_int_is_zero(bmap->div[div][0]))
1182	return isl_bool_true;
1183
1184	if (isl_seq_last_non_zero(p: bmap->ineq[ineq] + total + div + 1,
1185	len: bmap->n_div - (div + 1)) >= 0)
1186	return isl_bool_false;
1187
1188	last_ineq = isl_seq_last_non_zero(p: bmap->ineq[ineq], len: total + div);
1189	last_div = isl_seq_last_non_zero(p: bmap->div[div] + 1,
1190	len: total + bmap->n_div);
1191
1192	return last_ineq < last_div;
1193	}
1194
1195	/* Given two constraints "k" and "l" that are opposite to each other,
1196	* except for the constant term, check if we can use them
1197	* to obtain an expression for one of the hitherto unknown divs or
1198	* a "better" expression for a div for which we already have an expression.
1199	* "sum" is the sum of the constant terms of the constraints.
1200	* If this sum is strictly smaller than the coefficient of one
1201	* of the divs, then this pair can be used to define the div.
1202	* To avoid the introduction of circular definitions of divs, we
1203	* do not use the pair if the resulting expression would refer to
1204	* any other undefined divs or if any known div is defined in
1205	* terms of the unknown div.
1206	*/
1207	static __isl_give isl_basic_map *check_for_div_constraints(
1208	__isl_take isl_basic_map *bmap, int k, int l, isl_int sum,
1209	int *progress)
1210	{
1211	int i;
1212	unsigned total = isl_basic_map_offset(bmap, type: isl_dim_div);
1213
1214	for (i = 0; i < bmap->n_div; ++i) {
1215	isl_bool set_div;
1216
1217	if (isl_int_is_zero(bmap->ineq[k][total + i]))
1218	continue;
1219	if (isl_int_abs_ge(sum, bmap->ineq[k][total + i]))
1220	continue;
1221	set_div = better_div_constraint(bmap, div: i, ineq: k);
1222	if (set_div >= 0 && set_div)
1223	set_div = ok_to_set_div_from_bound(bmap, div: i, ineq: k);
1224	if (set_div < 0)
1225	return isl_basic_map_free(bmap);
1226	if (!set_div)
1227	break;
1228	if (isl_int_is_pos(bmap->ineq[k][total + i]))
1229	bmap = set_div_from_lower_bound(bmap, div: i, ineq: k);
1230	else
1231	bmap = set_div_from_lower_bound(bmap, div: i, ineq: l);
1232	if (progress)
1233	*progress = 1;
1234	break;
1235	}
1236	return bmap;
1237	}
1238
1239	__isl_give isl_basic_map *isl_basic_map_remove_duplicate_constraints(
1240	__isl_take isl_basic_map *bmap, int *progress, int detect_divs)
1241	{
1242	struct isl_constraint_index ci;
1243	int k, l, h;
1244	isl_size total = isl_basic_map_dim(bmap, type: isl_dim_all);
1245	isl_int sum;
1246
1247	if (total < 0 || bmap->n_ineq <= 1)
1248	return bmap;
1249
1250	if (create_constraint_index(ci: &ci, bmap) < 0)
1251	return bmap;
1252
1253	h = isl_seq_get_hash_bits(p: bmap->ineq[0] + 1, len: total, bits: ci.bits);
1254	ci.index[h] = &bmap->ineq[0];
1255	for (k = 1; k < bmap->n_ineq; ++k) {
1256	h = hash_index(ci: &ci, bmap, k);
1257	if (!ci.index[h]) {
1258	ci.index[h] = &bmap->ineq[k];
1259	continue;
1260	}
1261	if (progress)
1262	*progress = 1;
1263	l = ci.index[h] - &bmap->ineq[0];
1264	if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
1265	swap_inequality(bmap, a: k, b: l);
1266	isl_basic_map_drop_inequality(bmap, pos: k);
1267	--k;
1268	}
1269	isl_int_init(sum);
1270	for (k = 0; bmap && k < bmap->n_ineq-1; ++k) {
1271	isl_seq_neg(dst: bmap->ineq[k]+1, src: bmap->ineq[k]+1, len: total);
1272	h = hash_index(ci: &ci, bmap, k);
1273	isl_seq_neg(dst: bmap->ineq[k]+1, src: bmap->ineq[k]+1, len: total);
1274	if (!ci.index[h])
1275	continue;
1276	l = ci.index[h] - &bmap->ineq[0];
1277	isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
1278	if (isl_int_is_pos(sum)) {
1279	if (detect_divs)
1280	bmap = check_for_div_constraints(bmap, k, l,
1281	sum, progress);
1282	continue;
1283	}
1284	if (isl_int_is_zero(sum)) {
1285	/* We need to break out of the loop after these
1286	* changes since the contents of the hash
1287	* will no longer be valid.
1288	* Plus, we probably we want to regauss first.
1289	*/
1290	if (progress)
1291	*progress = 1;
1292	isl_basic_map_drop_inequality(bmap, pos: l);
1293	isl_basic_map_inequality_to_equality(bmap, pos: k);
1294	} else
1295	bmap = isl_basic_map_set_to_empty(bmap);
1296	break;
1297	}
1298	isl_int_clear(sum);
1299
1300	constraint_index_free(ci: &ci);
1301	return bmap;
1302	}
1303
1304	/* Detect all pairs of inequalities that form an equality.
1305	*
1306	* isl_basic_map_remove_duplicate_constraints detects at most one such pair.
1307	* Call it repeatedly while it is making progress.
1308	*/
1309	__isl_give isl_basic_map *isl_basic_map_detect_inequality_pairs(
1310	__isl_take isl_basic_map *bmap, int *progress)
1311	{
1312	int duplicate;
1313
1314	do {
1315	duplicate = 0;
1316	bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1317	progress: &duplicate, detect_divs: 0);
1318	if (progress && duplicate)
1319	*progress = 1;
1320	} while (duplicate);
1321
1322	return bmap;
1323	}
1324
1325	/* Given a known integer division "div" that is not integral
1326	* (with denominator 1), eliminate it from the constraints in "bmap"
1327	* where it appears with a (positive or negative) unit coefficient.
1328	* If "progress" is not NULL, then it gets set if the elimination
1329	* results in any changes.
1330	*
1331	* That is, replace
1332	*
1333	* floor(e/m) + f >= 0
1334	*
1335	* by
1336	*
1337	* e + m f >= 0
1338	*
1339	* and
1340	*
1341	* -floor(e/m) + f >= 0
1342	*
1343	* by
1344	*
1345	* -e + m f + m - 1 >= 0
1346	*
1347	* The first conversion is valid because floor(e/m) >= -f is equivalent
1348	* to e/m >= -f because -f is an integral expression.
1349	* The second conversion follows from the fact that
1350	*
1351	* -floor(e/m) = ceil(-e/m) = floor((-e + m - 1)/m)
1352	*
1353	*
1354	* Note that one of the div constraints may have been eliminated
1355	* due to being redundant with respect to the constraint that is
1356	* being modified by this function. The modified constraint may
1357	* no longer imply this div constraint, so we add it back to make
1358	* sure we do not lose any information.
1359	*/
1360	static __isl_give isl_basic_map *eliminate_unit_div(
1361	__isl_take isl_basic_map *bmap, int div, int *progress)
1362	{
1363	int j;
1364	isl_size v_div, dim;
1365	isl_ctx *ctx;
1366
1367	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
1368	dim = isl_basic_map_dim(bmap, type: isl_dim_all);
1369	if (v_div < 0 || dim < 0)
1370	return isl_basic_map_free(bmap);
1371
1372	ctx = isl_basic_map_get_ctx(bmap);
1373
1374	for (j = 0; j < bmap->n_ineq; ++j) {
1375	int s;
1376
1377	if (!isl_int_is_one(bmap->ineq[j][1 + v_div + div]) &&
1378	!isl_int_is_negone(bmap->ineq[j][1 + v_div + div]))
1379	continue;
1380
1381	if (progress)
1382	*progress = 1;
1383
1384	s = isl_int_sgn(bmap->ineq[j][1 + v_div + div]);
1385	isl_int_set_si(bmap->ineq[j][1 + v_div + div], 0);
1386	if (s < 0)
1387	isl_seq_combine(dst: bmap->ineq[j],
1388	m1: ctx->negone, src1: bmap->div[div] + 1,
1389	m2: bmap->div[div][0], src2: bmap->ineq[j], len: 1 + dim);
1390	else
1391	isl_seq_combine(dst: bmap->ineq[j],
1392	m1: ctx->one, src1: bmap->div[div] + 1,
1393	m2: bmap->div[div][0], src2: bmap->ineq[j], len: 1 + dim);
1394	if (s < 0) {
1395	isl_int_add(bmap->ineq[j][0],
1396	bmap->ineq[j][0], bmap->div[div][0]);
1397	isl_int_sub_ui(bmap->ineq[j][0],
1398	bmap->ineq[j][0], 1);
1399	}
1400
1401	bmap = isl_basic_map_extend_constraints(base: bmap, n_eq: 0, n_ineq: 1);
1402	bmap = isl_basic_map_add_div_constraint(bmap, div, sign: s);
1403	if (!bmap)
1404	return NULL;
1405	}
1406
1407	return bmap;
1408	}
1409
1410	/* Eliminate selected known divs from constraints where they appear with
1411	* a (positive or negative) unit coefficient.
1412	* In particular, only handle those for which "select" returns isl_bool_true.
1413	* If "progress" is not NULL, then it gets set if the elimination
1414	* results in any changes.
1415	*
1416	* We skip integral divs, i.e., those with denominator 1, as we would
1417	* risk eliminating the div from the div constraints. We do not need
1418	* to handle those divs here anyway since the div constraints will turn
1419	* out to form an equality and this equality can then be used to eliminate
1420	* the div from all constraints.
1421	*/
1422	static __isl_give isl_basic_map *eliminate_selected_unit_divs(
1423	__isl_take isl_basic_map *bmap,
1424	isl_bool (*select)(__isl_keep isl_basic_map *bmap, int div),
1425	int *progress)
1426	{
1427	int i;
1428
1429	if (!bmap)
1430	return NULL;
1431
1432	for (i = 0; i < bmap->n_div; ++i) {
1433	isl_bool selected;
1434
1435	if (isl_int_is_zero(bmap->div[i][0]))
1436	continue;
1437	if (isl_int_is_one(bmap->div[i][0]))
1438	continue;
1439	selected = select(bmap, i);
1440	if (selected < 0)
1441	return isl_basic_map_free(bmap);
1442	if (!selected)
1443	continue;
1444	bmap = eliminate_unit_div(bmap, div: i, progress);
1445	if (!bmap)
1446	return NULL;
1447	}
1448
1449	return bmap;
1450	}
1451
1452	/* eliminate_selected_unit_divs callback that selects every
1453	* integer division.
1454	*/
1455	static isl_bool is_any_div(__isl_keep isl_basic_map *bmap, int div)
1456	{
1457	return isl_bool_true;
1458	}
1459
1460	/* Eliminate known divs from constraints where they appear with
1461	* a (positive or negative) unit coefficient.
1462	* If "progress" is not NULL, then it gets set if the elimination
1463	* results in any changes.
1464	*/
1465	static __isl_give isl_basic_map *eliminate_unit_divs(
1466	__isl_take isl_basic_map *bmap, int *progress)
1467	{
1468	return eliminate_selected_unit_divs(bmap, select: &is_any_div, progress);
1469	}
1470
1471	/* eliminate_selected_unit_divs callback that selects
1472	* integer divisions that only appear with
1473	* a (positive or negative) unit coefficient
1474	* (outside their div constraints).
1475	*/
1476	static isl_bool is_pure_unit_div(__isl_keep isl_basic_map *bmap, int div)
1477	{
1478	int i;
1479	isl_size v_div, n_ineq;
1480
1481	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
1482	n_ineq = isl_basic_map_n_inequality(bmap);
1483	if (v_div < 0 || n_ineq < 0)
1484	return isl_bool_error;
1485
1486	for (i = 0; i < n_ineq; ++i) {
1487	isl_bool skip;
1488
1489	if (isl_int_is_zero(bmap->ineq[i][1 + v_div + div]))
1490	continue;
1491	skip = isl_basic_map_is_div_constraint(bmap,
1492	constraint: bmap->ineq[i], div);
1493	if (skip < 0)
1494	return isl_bool_error;
1495	if (skip)
1496	continue;
1497	if (!isl_int_is_one(bmap->ineq[i][1 + v_div + div]) &&
1498	!isl_int_is_negone(bmap->ineq[i][1 + v_div + div]))
1499	return isl_bool_false;
1500	}
1501
1502	return isl_bool_true;
1503	}
1504
1505	/* Eliminate known divs from constraints where they appear with
1506	* a (positive or negative) unit coefficient,
1507	* but only if they do not appear in any other constraints
1508	* (other than the div constraints).
1509	*/
1510	__isl_give isl_basic_map *isl_basic_map_eliminate_pure_unit_divs(
1511	__isl_take isl_basic_map *bmap)
1512	{
1513	return eliminate_selected_unit_divs(bmap, select: &is_pure_unit_div, NULL);
1514	}
1515
1516	__isl_give isl_basic_map *isl_basic_map_simplify(__isl_take isl_basic_map *bmap)
1517	{
1518	int progress = 1;
1519	if (!bmap)
1520	return NULL;
1521	while (progress) {
1522	isl_bool empty;
1523
1524	progress = 0;
1525	empty = isl_basic_map_plain_is_empty(bmap);
1526	if (empty < 0)
1527	return isl_basic_map_free(bmap);
1528	if (empty)
1529	break;
1530	bmap = isl_basic_map_normalize_constraints(bmap);
1531	bmap = reduce_div_coefficients(bmap);
1532	bmap = normalize_div_expressions(bmap);
1533	bmap = remove_duplicate_divs(bmap, progress: &progress);
1534	bmap = eliminate_unit_divs(bmap, progress: &progress);
1535	bmap = eliminate_divs_eq(bmap, progress: &progress);
1536	bmap = eliminate_divs_ineq(bmap, progress: &progress);
1537	bmap = isl_basic_map_gauss(bmap, progress: &progress);
1538	/* requires equalities in normal form */
1539	bmap = normalize_divs(bmap, progress: &progress);
1540	bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1541	progress: &progress, detect_divs: 1);
1542	if (bmap && progress)
1543	ISL_F_CLR(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
1544	}
1545	return bmap;
1546	}
1547
1548	__isl_give isl_basic_set *isl_basic_set_simplify(
1549	__isl_take isl_basic_set *bset)
1550	{
1551	return bset_from_bmap(bmap: isl_basic_map_simplify(bmap: bset_to_bmap(bset)));
1552	}
1553
1554
1555	isl_bool isl_basic_map_is_div_constraint(__isl_keep isl_basic_map *bmap,
1556	isl_int *constraint, unsigned div)
1557	{
1558	unsigned pos;
1559
1560	if (!bmap)
1561	return isl_bool_error;
1562
1563	pos = isl_basic_map_offset(bmap, type: isl_dim_div) + div;
1564
1565	if (isl_int_eq(constraint[pos], bmap->div[div][0])) {
1566	int neg;
1567	isl_int_sub(bmap->div[div][1],
1568	bmap->div[div][1], bmap->div[div][0]);
1569	isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
1570	neg = isl_seq_is_neg(p1: constraint, p2: bmap->div[div]+1, len: pos);
1571	isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
1572	isl_int_add(bmap->div[div][1],
1573	bmap->div[div][1], bmap->div[div][0]);
1574	if (!neg)
1575	return isl_bool_false;
1576	if (isl_seq_first_non_zero(p: constraint+pos+1,
1577	len: bmap->n_div-div-1) != -1)
1578	return isl_bool_false;
1579	} else if (isl_int_abs_eq(constraint[pos], bmap->div[div][0])) {
1580	if (!isl_seq_eq(p1: constraint, p2: bmap->div[div]+1, len: pos))
1581	return isl_bool_false;
1582	if (isl_seq_first_non_zero(p: constraint+pos+1,
1583	len: bmap->n_div-div-1) != -1)
1584	return isl_bool_false;
1585	} else
1586	return isl_bool_false;
1587
1588	return isl_bool_true;
1589	}
1590
1591	/* If the only constraints a div d=floor(f/m)
1592	* appears in are its two defining constraints
1593	*
1594	* f - m d >=0
1595	* -(f - (m - 1)) + m d >= 0
1596	*
1597	* then it can safely be removed.
1598	*/
1599	static isl_bool div_is_redundant(__isl_keep isl_basic_map *bmap, int div)
1600	{
1601	int i;
1602	isl_size v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
1603	unsigned pos = 1 + v_div + div;
1604
1605	if (v_div < 0)
1606	return isl_bool_error;
1607
1608	for (i = 0; i < bmap->n_eq; ++i)
1609	if (!isl_int_is_zero(bmap->eq[i][pos]))
1610	return isl_bool_false;
1611
1612	for (i = 0; i < bmap->n_ineq; ++i) {
1613	isl_bool red;
1614
1615	if (isl_int_is_zero(bmap->ineq[i][pos]))
1616	continue;
1617	red = isl_basic_map_is_div_constraint(bmap, constraint: bmap->ineq[i], div);
1618	if (red < 0 || !red)
1619	return red;
1620	}
1621
1622	for (i = 0; i < bmap->n_div; ++i) {
1623	if (isl_int_is_zero(bmap->div[i][0]))
1624	continue;
1625	if (!isl_int_is_zero(bmap->div[i][1+pos]))
1626	return isl_bool_false;
1627	}
1628
1629	return isl_bool_true;
1630	}
1631
1632	/*
1633	* Remove divs that don't occur in any of the constraints or other divs.
1634	* These can arise when dropping constraints from a basic map or
1635	* when the divs of a basic map have been temporarily aligned
1636	* with the divs of another basic map.
1637	*/
1638	static __isl_give isl_basic_map *remove_redundant_divs(
1639	__isl_take isl_basic_map *bmap)
1640	{
1641	int i;
1642	isl_size v_div;
1643
1644	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
1645	if (v_div < 0)
1646	return isl_basic_map_free(bmap);
1647
1648	for (i = bmap->n_div-1; i >= 0; --i) {
1649	isl_bool redundant;
1650
1651	redundant = div_is_redundant(bmap, div: i);
1652	if (redundant < 0)
1653	return isl_basic_map_free(bmap);
1654	if (!redundant)
1655	continue;
1656	bmap = isl_basic_map_drop_constraints_involving(bmap,
1657	first: v_div + i, n: 1);
1658	bmap = isl_basic_map_drop_div(bmap, div: i);
1659	}
1660	return bmap;
1661	}
1662
1663	/* Mark "bmap" as final, without checking for obviously redundant
1664	* integer divisions. This function should be used when "bmap"
1665	* is known not to involve any such integer divisions.
1666	*/
1667	__isl_give isl_basic_map *isl_basic_map_mark_final(
1668	__isl_take isl_basic_map *bmap)
1669	{
1670	if (!bmap)
1671	return NULL;
1672	ISL_F_SET(bmap, ISL_BASIC_SET_FINAL);
1673	return bmap;
1674	}
1675
1676	/* Mark "bmap" as final, after removing obviously redundant integer divisions.
1677	*/
1678	__isl_give isl_basic_map *isl_basic_map_finalize(__isl_take isl_basic_map *bmap)
1679	{
1680	bmap = remove_redundant_divs(bmap);
1681	bmap = isl_basic_map_mark_final(bmap);
1682	return bmap;
1683	}
1684
1685	__isl_give isl_basic_set *isl_basic_set_finalize(
1686	__isl_take isl_basic_set *bset)
1687	{
1688	return bset_from_bmap(bmap: isl_basic_map_finalize(bmap: bset_to_bmap(bset)));
1689	}
1690
1691	/* Remove definition of any div that is defined in terms of the given variable.
1692	* The div itself is not removed. Functions such as
1693	* eliminate_divs_ineq depend on the other divs remaining in place.
1694	*/
1695	static __isl_give isl_basic_map *remove_dependent_vars(
1696	__isl_take isl_basic_map *bmap, int pos)
1697	{
1698	int i;
1699
1700	if (!bmap)
1701	return NULL;
1702
1703	for (i = 0; i < bmap->n_div; ++i) {
1704	if (isl_int_is_zero(bmap->div[i][0]))
1705	continue;
1706	if (isl_int_is_zero(bmap->div[i][1+1+pos]))
1707	continue;
1708	bmap = isl_basic_map_mark_div_unknown(bmap, div: i);
1709	if (!bmap)
1710	return NULL;
1711	}
1712	return bmap;
1713	}
1714
1715	/* Eliminate the specified variables from the constraints using
1716	* Fourier-Motzkin. The variables themselves are not removed.
1717	*/
1718	__isl_give isl_basic_map *isl_basic_map_eliminate_vars(
1719	__isl_take isl_basic_map *bmap, unsigned pos, unsigned n)
1720	{
1721	int d;
1722	int i, j, k;
1723	isl_size total;
1724	int need_gauss = 0;
1725
1726	if (n == 0)
1727	return bmap;
1728	total = isl_basic_map_dim(bmap, type: isl_dim_all);
1729	if (total < 0)
1730	return isl_basic_map_free(bmap);
1731
1732	bmap = isl_basic_map_cow(bmap);
1733	for (d = pos + n - 1; d >= 0 && d >= pos; --d)
1734	bmap = remove_dependent_vars(bmap, pos: d);
1735	if (!bmap)
1736	return NULL;
1737
1738	for (d = pos + n - 1;
1739	d >= 0 && d >= total - bmap->n_div && d >= pos; --d)
1740	isl_seq_clr(p: bmap->div[d-(total-bmap->n_div)], len: 2+total);
1741	for (d = pos + n - 1; d >= 0 && d >= pos; --d) {
1742	int n_lower, n_upper;
1743	if (!bmap)
1744	return NULL;
1745	for (i = 0; i < bmap->n_eq; ++i) {
1746	if (isl_int_is_zero(bmap->eq[i][1+d]))
1747	continue;
1748	bmap = eliminate_var_using_equality(bmap, pos: d,
1749	eq: bmap->eq[i], keep_divs: 0, NULL);
1750	if (isl_basic_map_drop_equality(bmap, pos: i) < 0)
1751	return isl_basic_map_free(bmap);
1752	need_gauss = 1;
1753	break;
1754	}
1755	if (i < bmap->n_eq)
1756	continue;
1757	n_lower = 0;
1758	n_upper = 0;
1759	for (i = 0; i < bmap->n_ineq; ++i) {
1760	if (isl_int_is_pos(bmap->ineq[i][1+d]))
1761	n_lower++;
1762	else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1763	n_upper++;
1764	}
1765	bmap = isl_basic_map_extend_constraints(base: bmap,
1766	n_eq: 0, n_ineq: n_lower * n_upper);
1767	if (!bmap)
1768	goto error;
1769	for (i = bmap->n_ineq - 1; i >= 0; --i) {
1770	int last;
1771	if (isl_int_is_zero(bmap->ineq[i][1+d]))
1772	continue;
1773	last = -1;
1774	for (j = 0; j < i; ++j) {
1775	if (isl_int_is_zero(bmap->ineq[j][1+d]))
1776	continue;
1777	last = j;
1778	if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1779	isl_int_sgn(bmap->ineq[j][1+d]))
1780	continue;
1781	k = isl_basic_map_alloc_inequality(bmap);
1782	if (k < 0)
1783	goto error;
1784	isl_seq_cpy(dst: bmap->ineq[k], src: bmap->ineq[i],
1785	len: 1+total);
1786	isl_seq_elim(dst: bmap->ineq[k], src: bmap->ineq[j],
1787	pos: 1+d, len: 1+total, NULL);
1788	}
1789	isl_basic_map_drop_inequality(bmap, pos: i);
1790	i = last + 1;
1791	}
1792	if (n_lower > 0 && n_upper > 0) {
1793	bmap = isl_basic_map_normalize_constraints(bmap);
1794	bmap = isl_basic_map_remove_duplicate_constraints(bmap,
1795	NULL, detect_divs: 0);
1796	bmap = isl_basic_map_gauss(bmap, NULL);
1797	bmap = isl_basic_map_remove_redundancies(bmap);
1798	need_gauss = 0;
1799	if (!bmap)
1800	goto error;
1801	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1802	break;
1803	}
1804	}
1805	if (need_gauss)
1806	bmap = isl_basic_map_gauss(bmap, NULL);
1807	return bmap;
1808	error:
1809	isl_basic_map_free(bmap);
1810	return NULL;
1811	}
1812
1813	__isl_give isl_basic_set *isl_basic_set_eliminate_vars(
1814	__isl_take isl_basic_set *bset, unsigned pos, unsigned n)
1815	{
1816	return bset_from_bmap(bmap: isl_basic_map_eliminate_vars(bmap: bset_to_bmap(bset),
1817	pos, n));
1818	}
1819
1820	/* Eliminate the specified n dimensions starting at first from the
1821	* constraints, without removing the dimensions from the space.
1822	* If the set is rational, the dimensions are eliminated using Fourier-Motzkin.
1823	* Otherwise, they are projected out and the original space is restored.
1824	*/
1825	__isl_give isl_basic_map *isl_basic_map_eliminate(
1826	__isl_take isl_basic_map *bmap,
1827	enum isl_dim_type type, unsigned first, unsigned n)
1828	{
1829	isl_space *space;
1830
1831	if (!bmap)
1832	return NULL;
1833	if (n == 0)
1834	return bmap;
1835
1836	if (isl_basic_map_check_range(bmap, type, first, n) < 0)
1837	return isl_basic_map_free(bmap);
1838
1839	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
1840	first += isl_basic_map_offset(bmap, type) - 1;
1841	bmap = isl_basic_map_eliminate_vars(bmap, pos: first, n);
1842	return isl_basic_map_finalize(bmap);
1843	}
1844
1845	space = isl_basic_map_get_space(bmap);
1846	bmap = isl_basic_map_project_out(bmap, type, first, n);
1847	bmap = isl_basic_map_insert_dims(bmap, type, pos: first, n);
1848	bmap = isl_basic_map_reset_space(bmap, space);
1849	return bmap;
1850	}
1851
1852	__isl_give isl_basic_set *isl_basic_set_eliminate(
1853	__isl_take isl_basic_set *bset,
1854	enum isl_dim_type type, unsigned first, unsigned n)
1855	{
1856	return isl_basic_map_eliminate(bmap: bset, type, first, n);
1857	}
1858
1859	/* Remove all constraints from "bmap" that reference any unknown local
1860	* variables (directly or indirectly).
1861	*
1862	* Dropping all constraints on a local variable will make it redundant,
1863	* so it will get removed implicitly by
1864	* isl_basic_map_drop_constraints_involving_dims. Some other local
1865	* variables may also end up becoming redundant if they only appear
1866	* in constraints together with the unknown local variable.
1867	* Therefore, start over after calling
1868	* isl_basic_map_drop_constraints_involving_dims.
1869	*/
1870	__isl_give isl_basic_map *isl_basic_map_drop_constraints_involving_unknown_divs(
1871	__isl_take isl_basic_map *bmap)
1872	{
1873	isl_bool known;
1874	isl_size n_div;
1875	int i, o_div;
1876
1877	known = isl_basic_map_divs_known(bmap);
1878	if (known < 0)
1879	return isl_basic_map_free(bmap);
1880	if (known)
1881	return bmap;
1882
1883	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
1884	if (n_div < 0)
1885	return isl_basic_map_free(bmap);
1886	o_div = isl_basic_map_offset(bmap, type: isl_dim_div) - 1;
1887
1888	for (i = 0; i < n_div; ++i) {
1889	known = isl_basic_map_div_is_known(bmap, div: i);
1890	if (known < 0)
1891	return isl_basic_map_free(bmap);
1892	if (known)
1893	continue;
1894	bmap = remove_dependent_vars(bmap, pos: o_div + i);
1895	bmap = isl_basic_map_drop_constraints_involving_dims(bmap,
1896	type: isl_dim_div, first: i, n: 1);
1897	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
1898	if (n_div < 0)
1899	return isl_basic_map_free(bmap);
1900	i = -1;
1901	}
1902
1903	return bmap;
1904	}
1905
1906	/* Remove all constraints from "bset" that reference any unknown local
1907	* variables (directly or indirectly).
1908	*/
1909	__isl_give isl_basic_set *isl_basic_set_drop_constraints_involving_unknown_divs(
1910	__isl_take isl_basic_set *bset)
1911	{
1912	isl_basic_map *bmap;
1913
1914	bmap = bset_to_bmap(bset);
1915	bmap = isl_basic_map_drop_constraints_involving_unknown_divs(bmap);
1916	return bset_from_bmap(bmap);
1917	}
1918
1919	/* Remove all constraints from "map" that reference any unknown local
1920	* variables (directly or indirectly).
1921	*
1922	* Since constraints may get dropped from the basic maps,
1923	* they may no longer be disjoint from each other.
1924	*/
1925	__isl_give isl_map *isl_map_drop_constraints_involving_unknown_divs(
1926	__isl_take isl_map *map)
1927	{
1928	int i;
1929	isl_bool known;
1930
1931	known = isl_map_divs_known(map);
1932	if (known < 0)
1933	return isl_map_free(map);
1934	if (known)
1935	return map;
1936
1937	map = isl_map_cow(map);
1938	if (!map)
1939	return NULL;
1940
1941	for (i = 0; i < map->n; ++i) {
1942	map->p[i] =
1943	isl_basic_map_drop_constraints_involving_unknown_divs(
1944	bmap: map->p[i]);
1945	if (!map->p[i])
1946	return isl_map_free(map);
1947	}
1948
1949	if (map->n > 1)
1950	ISL_F_CLR(map, ISL_MAP_DISJOINT);
1951
1952	return map;
1953	}
1954
1955	/* Don't assume equalities are in order, because align_divs
1956	* may have changed the order of the divs.
1957	*/
1958	static void compute_elimination_index(__isl_keep isl_basic_map *bmap, int *elim,
1959	unsigned len)
1960	{
1961	int d, i;
1962
1963	for (d = 0; d < len; ++d)
1964	elim[d] = -1;
1965	for (i = 0; i < bmap->n_eq; ++i) {
1966	for (d = len - 1; d >= 0; --d) {
1967	if (isl_int_is_zero(bmap->eq[i][1+d]))
1968	continue;
1969	elim[d] = i;
1970	break;
1971	}
1972	}
1973	}
1974
1975	static void set_compute_elimination_index(__isl_keep isl_basic_set *bset,
1976	int *elim, unsigned len)
1977	{
1978	compute_elimination_index(bmap: bset_to_bmap(bset), elim, len);
1979	}
1980
1981	static int reduced_using_equalities(isl_int *dst, isl_int *src,
1982	__isl_keep isl_basic_map *bmap, int *elim, unsigned total)
1983	{
1984	int d;
1985	int copied = 0;
1986
1987	for (d = total - 1; d >= 0; --d) {
1988	if (isl_int_is_zero(src[1+d]))
1989	continue;
1990	if (elim[d] == -1)
1991	continue;
1992	if (!copied) {
1993	isl_seq_cpy(dst, src, len: 1 + total);
1994	copied = 1;
1995	}
1996	isl_seq_elim(dst, src: bmap->eq[elim[d]], pos: 1 + d, len: 1 + total, NULL);
1997	}
1998	return copied;
1999	}
2000
2001	static int set_reduced_using_equalities(isl_int *dst, isl_int *src,
2002	__isl_keep isl_basic_set *bset, int *elim, unsigned total)
2003	{
2004	return reduced_using_equalities(dst, src,
2005	bmap: bset_to_bmap(bset), elim, total);
2006	}
2007
2008	static __isl_give isl_basic_set *isl_basic_set_reduce_using_equalities(
2009	__isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
2010	{
2011	int i;
2012	int *elim;
2013	isl_size dim;
2014
2015	if (!bset || !context)
2016	goto error;
2017
2018	if (context->n_eq == 0) {
2019	isl_basic_set_free(bset: context);
2020	return bset;
2021	}
2022
2023	bset = isl_basic_set_cow(bset);
2024	dim = isl_basic_set_dim(bset, type: isl_dim_set);
2025	if (dim < 0)
2026	goto error;
2027
2028	elim = isl_alloc_array(bset->ctx, int, dim);
2029	if (!elim)
2030	goto error;
2031	set_compute_elimination_index(bset: context, elim, len: dim);
2032	for (i = 0; i < bset->n_eq; ++i)
2033	set_reduced_using_equalities(dst: bset->eq[i], src: bset->eq[i],
2034	bset: context, elim, total: dim);
2035	for (i = 0; i < bset->n_ineq; ++i)
2036	set_reduced_using_equalities(dst: bset->ineq[i], src: bset->ineq[i],
2037	bset: context, elim, total: dim);
2038	isl_basic_set_free(bset: context);
2039	free(ptr: elim);
2040	bset = isl_basic_set_simplify(bset);
2041	bset = isl_basic_set_finalize(bset);
2042	return bset;
2043	error:
2044	isl_basic_set_free(bset);
2045	isl_basic_set_free(bset: context);
2046	return NULL;
2047	}
2048
2049	/* For each inequality in "ineq" that is a shifted (more relaxed)
2050	* copy of an inequality in "context", mark the corresponding entry
2051	* in "row" with -1.
2052	* If an inequality only has a non-negative constant term, then
2053	* mark it as well.
2054	*/
2055	static isl_stat mark_shifted_constraints(__isl_keep isl_mat *ineq,
2056	__isl_keep isl_basic_set *context, int *row)
2057	{
2058	struct isl_constraint_index ci;
2059	isl_size n_ineq, cols;
2060	unsigned total;
2061	int k;
2062
2063	if (!ineq || !context)
2064	return isl_stat_error;
2065	if (context->n_ineq == 0)
2066	return isl_stat_ok;
2067	if (setup_constraint_index(ci: &ci, bset: context) < 0)
2068	return isl_stat_error;
2069
2070	n_ineq = isl_mat_rows(mat: ineq);
2071	cols = isl_mat_cols(mat: ineq);
2072	if (n_ineq < 0 || cols < 0)
2073	return isl_stat_error;
2074	total = cols - 1;
2075	for (k = 0; k < n_ineq; ++k) {
2076	int l;
2077	isl_bool redundant;
2078
2079	l = isl_seq_first_non_zero(p: ineq->row[k] + 1, len: total);
2080	if (l < 0 && isl_int_is_nonneg(ineq->row[k][0])) {
2081	row[k] = -1;
2082	continue;
2083	}
2084	redundant = constraint_index_is_redundant(ci: &ci, ineq: ineq->row[k]);
2085	if (redundant < 0)
2086	goto error;
2087	if (!redundant)
2088	continue;
2089	row[k] = -1;
2090	}
2091	constraint_index_free(ci: &ci);
2092	return isl_stat_ok;
2093	error:
2094	constraint_index_free(ci: &ci);
2095	return isl_stat_error;
2096	}
2097
2098	static __isl_give isl_basic_set *remove_shifted_constraints(
2099	__isl_take isl_basic_set *bset, __isl_keep isl_basic_set *context)
2100	{
2101	struct isl_constraint_index ci;
2102	int k;
2103
2104	if (!bset || !context)
2105	return bset;
2106
2107	if (context->n_ineq == 0)
2108	return bset;
2109	if (setup_constraint_index(ci: &ci, bset: context) < 0)
2110	return bset;
2111
2112	for (k = 0; k < bset->n_ineq; ++k) {
2113	isl_bool redundant;
2114
2115	redundant = constraint_index_is_redundant(ci: &ci, ineq: bset->ineq[k]);
2116	if (redundant < 0)
2117	goto error;
2118	if (!redundant)
2119	continue;
2120	bset = isl_basic_set_cow(bset);
2121	if (!bset)
2122	goto error;
2123	isl_basic_set_drop_inequality(bset, pos: k);
2124	--k;
2125	}
2126	constraint_index_free(ci: &ci);
2127	return bset;
2128	error:
2129	constraint_index_free(ci: &ci);
2130	return bset;
2131	}
2132
2133	/* Remove constraints from "bmap" that are identical to constraints
2134	* in "context" or that are more relaxed (greater constant term).
2135	*
2136	* We perform the test for shifted copies on the pure constraints
2137	* in remove_shifted_constraints.
2138	*/
2139	static __isl_give isl_basic_map *isl_basic_map_remove_shifted_constraints(
2140	__isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
2141	{
2142	isl_basic_set *bset, *bset_context;
2143
2144	if (!bmap || !context)
2145	goto error;
2146
2147	if (bmap->n_ineq == 0 || context->n_ineq == 0) {
2148	isl_basic_map_free(bmap: context);
2149	return bmap;
2150	}
2151
2152	bmap = isl_basic_map_order_divs(bmap);
2153	context = isl_basic_map_align_divs(dst: context, src: bmap);
2154	bmap = isl_basic_map_align_divs(dst: bmap, src: context);
2155
2156	bset = isl_basic_map_underlying_set(bmap: isl_basic_map_copy(bmap));
2157	bset_context = isl_basic_map_underlying_set(bmap: context);
2158	bset = remove_shifted_constraints(bset, context: bset_context);
2159	isl_basic_set_free(bset: bset_context);
2160
2161	bmap = isl_basic_map_overlying_set(bset, like: bmap);
2162
2163	return bmap;
2164	error:
2165	isl_basic_map_free(bmap);
2166	isl_basic_map_free(bmap: context);
2167	return NULL;
2168	}
2169
2170	/* Does the (linear part of a) constraint "c" involve any of the "len"
2171	* "relevant" dimensions?
2172	*/
2173	static int is_related(isl_int *c, int len, int *relevant)
2174	{
2175	int i;
2176
2177	for (i = 0; i < len; ++i) {
2178	if (!relevant[i])
2179	continue;
2180	if (!isl_int_is_zero(c[i]))
2181	return 1;
2182	}
2183
2184	return 0;
2185	}
2186
2187	/* Drop constraints from "bmap" that do not involve any of
2188	* the dimensions marked "relevant".
2189	*/
2190	static __isl_give isl_basic_map *drop_unrelated_constraints(
2191	__isl_take isl_basic_map *bmap, int *relevant)
2192	{
2193	int i;
2194	isl_size dim;
2195
2196	dim = isl_basic_map_dim(bmap, type: isl_dim_all);
2197	if (dim < 0)
2198	return isl_basic_map_free(bmap);
2199	for (i = 0; i < dim; ++i)
2200	if (!relevant[i])
2201	break;
2202	if (i >= dim)
2203	return bmap;
2204
2205	for (i = bmap->n_eq - 1; i >= 0; --i)
2206	if (!is_related(c: bmap->eq[i] + 1, len: dim, relevant)) {
2207	bmap = isl_basic_map_cow(bmap);
2208	if (isl_basic_map_drop_equality(bmap, pos: i) < 0)
2209	return isl_basic_map_free(bmap);
2210	}
2211
2212	for (i = bmap->n_ineq - 1; i >= 0; --i)
2213	if (!is_related(c: bmap->ineq[i] + 1, len: dim, relevant)) {
2214	bmap = isl_basic_map_cow(bmap);
2215	if (isl_basic_map_drop_inequality(bmap, pos: i) < 0)
2216	return isl_basic_map_free(bmap);
2217	}
2218
2219	return bmap;
2220	}
2221
2222	/* Update the groups in "group" based on the (linear part of a) constraint "c".
2223	*
2224	* In particular, for any variable involved in the constraint,
2225	* find the actual group id from before and replace the group
2226	* of the corresponding variable by the minimal group of all
2227	* the variables involved in the constraint considered so far
2228	* (if this minimum is smaller) or replace the minimum by this group
2229	* (if the minimum is larger).
2230	*
2231	* At the end, all the variables in "c" will (indirectly) point
2232	* to the minimal of the groups that they referred to originally.
2233	*/
2234	static void update_groups(int dim, int *group, isl_int *c)
2235	{
2236	int j;
2237	int min = dim;
2238
2239	for (j = 0; j < dim; ++j) {
2240	if (isl_int_is_zero(c[j]))
2241	continue;
2242	while (group[j] >= 0 && group[group[j]] != group[j])
2243	group[j] = group[group[j]];
2244	if (group[j] == min)
2245	continue;
2246	if (group[j] < min) {
2247	if (min >= 0 && min < dim)
2248	group[min] = group[j];
2249	min = group[j];
2250	} else
2251	group[group[j]] = min;
2252	}
2253	}
2254
2255	/* Allocate an array of groups of variables, one for each variable
2256	* in "context", initialized to zero.
2257	*/
2258	static int *alloc_groups(__isl_keep isl_basic_set *context)
2259	{
2260	isl_ctx *ctx;
2261	isl_size dim;
2262
2263	dim = isl_basic_set_dim(bset: context, type: isl_dim_set);
2264	if (dim < 0)
2265	return NULL;
2266	ctx = isl_basic_set_get_ctx(bset: context);
2267	return isl_calloc_array(ctx, int, dim);
2268	}
2269
2270	/* Drop constraints from "bmap" that only involve variables that are
2271	* not related to any of the variables marked with a "-1" in "group".
2272	*
2273	* We construct groups of variables that collect variables that
2274	* (indirectly) appear in some common constraint of "bmap".
2275	* Each group is identified by the first variable in the group,
2276	* except for the special group of variables that was already identified
2277	* in the input as -1 (or are related to those variables).
2278	* If group[i] is equal to i (or -1), then the group of i is i (or -1),
2279	* otherwise the group of i is the group of group[i].
2280	*
2281	* We first initialize groups for the remaining variables.
2282	* Then we iterate over the constraints of "bmap" and update the
2283	* group of the variables in the constraint by the smallest group.
2284	* Finally, we resolve indirect references to groups by running over
2285	* the variables.
2286	*
2287	* After computing the groups, we drop constraints that do not involve
2288	* any variables in the -1 group.
2289	*/
2290	__isl_give isl_basic_map *isl_basic_map_drop_unrelated_constraints(
2291	__isl_take isl_basic_map *bmap, __isl_take int *group)
2292	{
2293	isl_size dim;
2294	int i;
2295	int last;
2296
2297	dim = isl_basic_map_dim(bmap, type: isl_dim_all);
2298	if (dim < 0)
2299	return isl_basic_map_free(bmap);
2300
2301	last = -1;
2302	for (i = 0; i < dim; ++i)
2303	if (group[i] >= 0)
2304	last = group[i] = i;
2305	if (last < 0) {
2306	free(ptr: group);
2307	return bmap;
2308	}
2309
2310	for (i = 0; i < bmap->n_eq; ++i)
2311	update_groups(dim, group, c: bmap->eq[i] + 1);
2312	for (i = 0; i < bmap->n_ineq; ++i)
2313	update_groups(dim, group, c: bmap->ineq[i] + 1);
2314
2315	for (i = 0; i < dim; ++i)
2316	if (group[i] >= 0)
2317	group[i] = group[group[i]];
2318
2319	for (i = 0; i < dim; ++i)
2320	group[i] = group[i] == -1;
2321
2322	bmap = drop_unrelated_constraints(bmap, relevant: group);
2323
2324	free(ptr: group);
2325	return bmap;
2326	}
2327
2328	/* Drop constraints from "context" that are irrelevant for computing
2329	* the gist of "bset".
2330	*
2331	* In particular, drop constraints in variables that are not related
2332	* to any of the variables involved in the constraints of "bset"
2333	* in the sense that there is no sequence of constraints that connects them.
2334	*
2335	* We first mark all variables that appear in "bset" as belonging
2336	* to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2337	*/
2338	static __isl_give isl_basic_set *drop_irrelevant_constraints(
2339	__isl_take isl_basic_set *context, __isl_keep isl_basic_set *bset)
2340	{
2341	int *group;
2342	isl_size dim;
2343	int i, j;
2344
2345	dim = isl_basic_set_dim(bset, type: isl_dim_set);
2346	if (!context || dim < 0)
2347	return isl_basic_set_free(bset: context);
2348
2349	group = alloc_groups(context);
2350
2351	if (!group)
2352	return isl_basic_set_free(bset: context);
2353
2354	for (i = 0; i < dim; ++i) {
2355	for (j = 0; j < bset->n_eq; ++j)
2356	if (!isl_int_is_zero(bset->eq[j][1 + i]))
2357	break;
2358	if (j < bset->n_eq) {
2359	group[i] = -1;
2360	continue;
2361	}
2362	for (j = 0; j < bset->n_ineq; ++j)
2363	if (!isl_int_is_zero(bset->ineq[j][1 + i]))
2364	break;
2365	if (j < bset->n_ineq)
2366	group[i] = -1;
2367	}
2368
2369	return isl_basic_map_drop_unrelated_constraints(bmap: context, group);
2370	}
2371
2372	/* Drop constraints from "context" that are irrelevant for computing
2373	* the gist of the inequalities "ineq".
2374	* Inequalities in "ineq" for which the corresponding element of row
2375	* is set to -1 have already been marked for removal and should be ignored.
2376	*
2377	* In particular, drop constraints in variables that are not related
2378	* to any of the variables involved in "ineq"
2379	* in the sense that there is no sequence of constraints that connects them.
2380	*
2381	* We first mark all variables that appear in "bset" as belonging
2382	* to a "-1" group and then continue with group_and_drop_irrelevant_constraints.
2383	*/
2384	static __isl_give isl_basic_set *drop_irrelevant_constraints_marked(
2385	__isl_take isl_basic_set *context, __isl_keep isl_mat *ineq, int *row)
2386	{
2387	int *group;
2388	isl_size dim;
2389	int i, j;
2390	isl_size n;
2391
2392	dim = isl_basic_set_dim(bset: context, type: isl_dim_set);
2393	n = isl_mat_rows(mat: ineq);
2394	if (dim < 0 || n < 0)
2395	return isl_basic_set_free(bset: context);
2396
2397	group = alloc_groups(context);
2398
2399	if (!group)
2400	return isl_basic_set_free(bset: context);
2401
2402	for (i = 0; i < dim; ++i) {
2403	for (j = 0; j < n; ++j) {
2404	if (row[j] < 0)
2405	continue;
2406	if (!isl_int_is_zero(ineq->row[j][1 + i]))
2407	break;
2408	}
2409	if (j < n)
2410	group[i] = -1;
2411	}
2412
2413	return isl_basic_map_drop_unrelated_constraints(bmap: context, group);
2414	}
2415
2416	/* Do all "n" entries of "row" contain a negative value?
2417	*/
2418	static int all_neg(int *row, int n)
2419	{
2420	int i;
2421
2422	for (i = 0; i < n; ++i)
2423	if (row[i] >= 0)
2424	return 0;
2425
2426	return 1;
2427	}
2428
2429	/* Update the inequalities in "bset" based on the information in "row"
2430	* and "tab".
2431	*
2432	* In particular, the array "row" contains either -1, meaning that
2433	* the corresponding inequality of "bset" is redundant, or the index
2434	* of an inequality in "tab".
2435	*
2436	* If the row entry is -1, then drop the inequality.
2437	* Otherwise, if the constraint is marked redundant in the tableau,
2438	* then drop the inequality. Similarly, if it is marked as an equality
2439	* in the tableau, then turn the inequality into an equality and
2440	* perform Gaussian elimination.
2441	*/
2442	static __isl_give isl_basic_set *update_ineq(__isl_take isl_basic_set *bset,
2443	__isl_keep int *row, struct isl_tab *tab)
2444	{
2445	int i;
2446	unsigned n_ineq;
2447	unsigned n_eq;
2448	int found_equality = 0;
2449
2450	if (!bset)
2451	return NULL;
2452	if (tab && tab->empty)
2453	return isl_basic_set_set_to_empty(bset);
2454
2455	n_ineq = bset->n_ineq;
2456	for (i = n_ineq - 1; i >= 0; --i) {
2457	if (row[i] < 0) {
2458	if (isl_basic_set_drop_inequality(bset, pos: i) < 0)
2459	return isl_basic_set_free(bset);
2460	continue;
2461	}
2462	if (!tab)
2463	continue;
2464	n_eq = tab->n_eq;
2465	if (isl_tab_is_equality(tab, con: n_eq + row[i])) {
2466	isl_basic_map_inequality_to_equality(bmap: bset, pos: i);
2467	found_equality = 1;
2468	} else if (isl_tab_is_redundant(tab, con: n_eq + row[i])) {
2469	if (isl_basic_set_drop_inequality(bset, pos: i) < 0)
2470	return isl_basic_set_free(bset);
2471	}
2472	}
2473
2474	if (found_equality)
2475	bset = isl_basic_set_gauss(bset, NULL);
2476	bset = isl_basic_set_finalize(bset);
2477	return bset;
2478	}
2479
2480	/* Update the inequalities in "bset" based on the information in "row"
2481	* and "tab" and free all arguments (other than "bset").
2482	*/
2483	static __isl_give isl_basic_set *update_ineq_free(
2484	__isl_take isl_basic_set *bset, __isl_take isl_mat *ineq,
2485	__isl_take isl_basic_set *context, __isl_take int *row,
2486	struct isl_tab *tab)
2487	{
2488	isl_mat_free(mat: ineq);
2489	isl_basic_set_free(bset: context);
2490
2491	bset = update_ineq(bset, row, tab);
2492
2493	free(ptr: row);
2494	isl_tab_free(tab);
2495	return bset;
2496	}
2497
2498	/* Remove all information from bset that is redundant in the context
2499	* of context.
2500	* "ineq" contains the (possibly transformed) inequalities of "bset",
2501	* in the same order.
2502	* The (explicit) equalities of "bset" are assumed to have been taken
2503	* into account by the transformation such that only the inequalities
2504	* are relevant.
2505	* "context" is assumed not to be empty.
2506	*
2507	* "row" keeps track of the constraint index of a "bset" inequality in "tab".
2508	* A value of -1 means that the inequality is obviously redundant and may
2509	* not even appear in "tab".
2510	*
2511	* We first mark the inequalities of "bset"
2512	* that are obviously redundant with respect to some inequality in "context".
2513	* Then we remove those constraints from "context" that have become
2514	* irrelevant for computing the gist of "bset".
2515	* Note that this removal of constraints cannot be replaced by
2516	* a factorization because factors in "bset" may still be connected
2517	* to each other through constraints in "context".
2518	*
2519	* If there are any inequalities left, we construct a tableau for
2520	* the context and then add the inequalities of "bset".
2521	* Before adding these inequalities, we freeze all constraints such that
2522	* they won't be considered redundant in terms of the constraints of "bset".
2523	* Then we detect all redundant constraints (among the
2524	* constraints that weren't frozen), first by checking for redundancy in the
2525	* the tableau and then by checking if replacing a constraint by its negation
2526	* would lead to an empty set. This last step is fairly expensive
2527	* and could be optimized by more reuse of the tableau.
2528	* Finally, we update bset according to the results.
2529	*/
2530	static __isl_give isl_basic_set *uset_gist_full(__isl_take isl_basic_set *bset,
2531	__isl_take isl_mat *ineq, __isl_take isl_basic_set *context)
2532	{
2533	int i, r;
2534	int *row = NULL;
2535	isl_ctx *ctx;
2536	isl_basic_set *combined = NULL;
2537	struct isl_tab *tab = NULL;
2538	unsigned n_eq, context_ineq;
2539
2540	if (!bset || !ineq || !context)
2541	goto error;
2542
2543	if (bset->n_ineq == 0 || isl_basic_set_plain_is_universe(bset: context)) {
2544	isl_basic_set_free(bset: context);
2545	isl_mat_free(mat: ineq);
2546	return bset;
2547	}
2548
2549	ctx = isl_basic_set_get_ctx(bset: context);
2550	row = isl_calloc_array(ctx, int, bset->n_ineq);
2551	if (!row)
2552	goto error;
2553
2554	if (mark_shifted_constraints(ineq, context, row) < 0)
2555	goto error;
2556	if (all_neg(row, n: bset->n_ineq))
2557	return update_ineq_free(bset, ineq, context, row, NULL);
2558
2559	context = drop_irrelevant_constraints_marked(context, ineq, row);
2560	if (!context)
2561	goto error;
2562	if (isl_basic_set_plain_is_universe(bset: context))
2563	return update_ineq_free(bset, ineq, context, row, NULL);
2564
2565	n_eq = context->n_eq;
2566	context_ineq = context->n_ineq;
2567	combined = isl_basic_set_cow(bset: isl_basic_set_copy(bset: context));
2568	combined = isl_basic_set_extend_constraints(base: combined, n_eq: 0, n_ineq: bset->n_ineq);
2569	tab = isl_tab_from_basic_set(bset: combined, track: 0);
2570	for (i = 0; i < context_ineq; ++i)
2571	if (isl_tab_freeze_constraint(tab, con: n_eq + i) < 0)
2572	goto error;
2573	if (isl_tab_extend_cons(tab, n_new: bset->n_ineq) < 0)
2574	goto error;
2575	r = context_ineq;
2576	for (i = 0; i < bset->n_ineq; ++i) {
2577	if (row[i] < 0)
2578	continue;
2579	combined = isl_basic_set_add_ineq(bset: combined, ineq: ineq->row[i]);
2580	if (isl_tab_add_ineq(tab, ineq: ineq->row[i]) < 0)
2581	goto error;
2582	row[i] = r++;
2583	}
2584	if (isl_tab_detect_implicit_equalities(tab) < 0)
2585	goto error;
2586	if (isl_tab_detect_redundant(tab) < 0)
2587	goto error;
2588	for (i = bset->n_ineq - 1; i >= 0; --i) {
2589	isl_basic_set *test;
2590	int is_empty;
2591
2592	if (row[i] < 0)
2593	continue;
2594	r = row[i];
2595	if (tab->con[n_eq + r].is_redundant)
2596	continue;
2597	test = isl_basic_set_dup(bset: combined);
2598	test = isl_inequality_negate(bmap: test, pos: r);
2599	test = isl_basic_set_update_from_tab(bset: test, tab);
2600	is_empty = isl_basic_set_is_empty(bset: test);
2601	isl_basic_set_free(bset: test);
2602	if (is_empty < 0)
2603	goto error;
2604	if (is_empty)
2605	tab->con[n_eq + r].is_redundant = 1;
2606	}
2607	bset = update_ineq_free(bset, ineq, context, row, tab);
2608	if (bset) {
2609	ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2610	ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2611	}
2612
2613	isl_basic_set_free(bset: combined);
2614	return bset;
2615	error:
2616	free(ptr: row);
2617	isl_mat_free(mat: ineq);
2618	isl_tab_free(tab);
2619	isl_basic_set_free(bset: combined);
2620	isl_basic_set_free(bset: context);
2621	isl_basic_set_free(bset);
2622	return NULL;
2623	}
2624
2625	/* Extract the inequalities of "bset" as an isl_mat.
2626	*/
2627	static __isl_give isl_mat *extract_ineq(__isl_keep isl_basic_set *bset)
2628	{
2629	isl_size total;
2630	isl_ctx *ctx;
2631	isl_mat *ineq;
2632
2633	total = isl_basic_set_dim(bset, type: isl_dim_all);
2634	if (total < 0)
2635	return NULL;
2636
2637	ctx = isl_basic_set_get_ctx(bset);
2638	ineq = isl_mat_sub_alloc6(ctx, row: bset->ineq, first_row: 0, n_row: bset->n_ineq,
2639	first_col: 0, n_col: 1 + total);
2640
2641	return ineq;
2642	}
2643
2644	/* Remove all information from "bset" that is redundant in the context
2645	* of "context", for the case where both "bset" and "context" are
2646	* full-dimensional.
2647	*/
2648	static __isl_give isl_basic_set *uset_gist_uncompressed(
2649	__isl_take isl_basic_set *bset, __isl_take isl_basic_set *context)
2650	{
2651	isl_mat *ineq;
2652
2653	ineq = extract_ineq(bset);
2654	return uset_gist_full(bset, ineq, context);
2655	}
2656
2657	/* Replace "bset" by an empty basic set in the same space.
2658	*/
2659	static __isl_give isl_basic_set *replace_by_empty(
2660	__isl_take isl_basic_set *bset)
2661	{
2662	isl_space *space;
2663
2664	space = isl_basic_set_get_space(bset);
2665	isl_basic_set_free(bset);
2666	return isl_basic_set_empty(space);
2667	}
2668
2669	/* Remove all information from "bset" that is redundant in the context
2670	* of "context", for the case where the combined equalities of
2671	* "bset" and "context" allow for a compression that can be obtained
2672	* by preapplication of "T".
2673	* If the compression of "context" is empty, meaning that "bset" and
2674	* "context" do not intersect, then return the empty set.
2675	*
2676	* "bset" itself is not transformed by "T". Instead, the inequalities
2677	* are extracted from "bset" and those are transformed by "T".
2678	* uset_gist_full then determines which of the transformed inequalities
2679	* are redundant with respect to the transformed "context" and removes
2680	* the corresponding inequalities from "bset".
2681	*
2682	* After preapplying "T" to the inequalities, any common factor is
2683	* removed from the coefficients. If this results in a tightening
2684	* of the constant term, then the same tightening is applied to
2685	* the corresponding untransformed inequality in "bset".
2686	* That is, if after plugging in T, a constraint f(x) >= 0 is of the form
2687	*
2688	* g f'(x) + r >= 0
2689	*
2690	* with 0 <= r < g, then it is equivalent to
2691	*
2692	* f'(x) >= 0
2693	*
2694	* This means that f(x) >= 0 is equivalent to f(x) - r >= 0 in the affine
2695	* subspace compressed by T since the latter would be transformed to
2696	*
2697	* g f'(x) >= 0
2698	*/
2699	static __isl_give isl_basic_set *uset_gist_compressed(
2700	__isl_take isl_basic_set *bset, __isl_take isl_basic_set *context,
2701	__isl_take isl_mat *T)
2702	{
2703	isl_ctx *ctx;
2704	isl_mat *ineq;
2705	int i;
2706	isl_size n_row, n_col;
2707	isl_int rem;
2708
2709	ineq = extract_ineq(bset);
2710	ineq = isl_mat_product(left: ineq, right: isl_mat_copy(mat: T));
2711	context = isl_basic_set_preimage(bset: context, mat: T);
2712
2713	if (!ineq || !context)
2714	goto error;
2715	if (isl_basic_set_plain_is_empty(bset: context)) {
2716	isl_mat_free(mat: ineq);
2717	isl_basic_set_free(bset: context);
2718	return replace_by_empty(bset);
2719	}
2720
2721	ctx = isl_mat_get_ctx(mat: ineq);
2722	n_row = isl_mat_rows(mat: ineq);
2723	n_col = isl_mat_cols(mat: ineq);
2724	if (n_row < 0 || n_col < 0)
2725	goto error;
2726	isl_int_init(rem);
2727	for (i = 0; i < n_row; ++i) {
2728	isl_seq_gcd(p: ineq->row[i] + 1, len: n_col - 1, gcd: &ctx->normalize_gcd);
2729	if (isl_int_is_zero(ctx->normalize_gcd))
2730	continue;
2731	if (isl_int_is_one(ctx->normalize_gcd))
2732	continue;
2733	isl_seq_scale_down(dst: ineq->row[i] + 1, src: ineq->row[i] + 1,
2734	f: ctx->normalize_gcd, len: n_col - 1);
2735	isl_int_fdiv_r(rem, ineq->row[i][0], ctx->normalize_gcd);
2736	isl_int_fdiv_q(ineq->row[i][0],
2737	ineq->row[i][0], ctx->normalize_gcd);
2738	if (isl_int_is_zero(rem))
2739	continue;
2740	bset = isl_basic_set_cow(bset);
2741	if (!bset)
2742	break;
2743	isl_int_sub(bset->ineq[i][0], bset->ineq[i][0], rem);
2744	}
2745	isl_int_clear(rem);
2746
2747	return uset_gist_full(bset, ineq, context);
2748	error:
2749	isl_mat_free(mat: ineq);
2750	isl_basic_set_free(bset: context);
2751	isl_basic_set_free(bset);
2752	return NULL;
2753	}
2754
2755	/* Project "bset" onto the variables that are involved in "template".
2756	*/
2757	static __isl_give isl_basic_set *project_onto_involved(
2758	__isl_take isl_basic_set *bset, __isl_keep isl_basic_set *template)
2759	{
2760	int i;
2761	isl_size n;
2762
2763	n = isl_basic_set_dim(bset: template, type: isl_dim_set);
2764	if (n < 0 || !template)
2765	return isl_basic_set_free(bset);
2766
2767	for (i = 0; i < n; ++i) {
2768	isl_bool involved;
2769
2770	involved = isl_basic_set_involves_dims(bset: template,
2771	type: isl_dim_set, first: i, n: 1);
2772	if (involved < 0)
2773	return isl_basic_set_free(bset);
2774	if (involved)
2775	continue;
2776	bset = isl_basic_set_eliminate_vars(bset, pos: i, n: 1);
2777	}
2778
2779	return bset;
2780	}
2781
2782	/* Remove all information from bset that is redundant in the context
2783	* of context. In particular, equalities that are linear combinations
2784	* of those in context are removed. Then the inequalities that are
2785	* redundant in the context of the equalities and inequalities of
2786	* context are removed.
2787	*
2788	* First of all, we drop those constraints from "context"
2789	* that are irrelevant for computing the gist of "bset".
2790	* Alternatively, we could factorize the intersection of "context" and "bset".
2791	*
2792	* We first compute the intersection of the integer affine hulls
2793	* of "bset" and "context",
2794	* compute the gist inside this intersection and then reduce
2795	* the constraints with respect to the equalities of the context
2796	* that only involve variables already involved in the input.
2797	* If the intersection of the affine hulls turns out to be empty,
2798	* then return the empty set.
2799	*
2800	* If two constraints are mutually redundant, then uset_gist_full
2801	* will remove the second of those constraints. We therefore first
2802	* sort the constraints so that constraints not involving existentially
2803	* quantified variables are given precedence over those that do.
2804	* We have to perform this sorting before the variable compression,
2805	* because that may effect the order of the variables.
2806	*/
2807	static __isl_give isl_basic_set *uset_gist(__isl_take isl_basic_set *bset,
2808	__isl_take isl_basic_set *context)
2809	{
2810	isl_mat *eq;
2811	isl_mat *T;
2812	isl_basic_set *aff;
2813	isl_basic_set *aff_context;
2814	isl_size total;
2815
2816	total = isl_basic_set_dim(bset, type: isl_dim_all);
2817	if (total < 0 || !context)
2818	goto error;
2819
2820	context = drop_irrelevant_constraints(context, bset);
2821
2822	bset = isl_basic_set_detect_equalities(bset);
2823	aff = isl_basic_set_copy(bset);
2824	aff = isl_basic_set_plain_affine_hull(bset: aff);
2825	context = isl_basic_set_detect_equalities(bset: context);
2826	aff_context = isl_basic_set_copy(bset: context);
2827	aff_context = isl_basic_set_plain_affine_hull(bset: aff_context);
2828	aff = isl_basic_set_intersect(bset1: aff, bset2: aff_context);
2829	if (!aff)
2830	goto error;
2831	if (isl_basic_set_plain_is_empty(bset: aff)) {
2832	isl_basic_set_free(bset);
2833	isl_basic_set_free(bset: context);
2834	return aff;
2835	}
2836	bset = isl_basic_set_sort_constraints(bset);
2837	if (aff->n_eq == 0) {
2838	isl_basic_set_free(bset: aff);
2839	return uset_gist_uncompressed(bset, context);
2840	}
2841	eq = isl_mat_sub_alloc6(ctx: bset->ctx, row: aff->eq, first_row: 0, n_row: aff->n_eq, first_col: 0, n_col: 1 + total);
2842	eq = isl_mat_cow(mat: eq);
2843	T = isl_mat_variable_compression(B: eq, NULL);
2844	isl_basic_set_free(bset: aff);
2845	if (T && T->n_col == 0) {
2846	isl_mat_free(mat: T);
2847	isl_basic_set_free(bset: context);
2848	return replace_by_empty(bset);
2849	}
2850
2851	aff_context = isl_basic_set_affine_hull(bset: isl_basic_set_copy(bset: context));
2852	aff_context = project_onto_involved(bset: aff_context, template: bset);
2853
2854	bset = uset_gist_compressed(bset, context, T);
2855	bset = isl_basic_set_reduce_using_equalities(bset, context: aff_context);
2856
2857	if (bset) {
2858	ISL_F_SET(bset, ISL_BASIC_SET_NO_IMPLICIT);
2859	ISL_F_SET(bset, ISL_BASIC_SET_NO_REDUNDANT);
2860	}
2861
2862	return bset;
2863	error:
2864	isl_basic_set_free(bset);
2865	isl_basic_set_free(bset: context);
2866	return NULL;
2867	}
2868
2869	/* Return the number of equality constraints in "bmap" that involve
2870	* local variables. This function assumes that Gaussian elimination
2871	* has been applied to the equality constraints.
2872	*/
2873	static int n_div_eq(__isl_keep isl_basic_map *bmap)
2874	{
2875	int i;
2876	isl_size total, n_div;
2877
2878	if (!bmap)
2879	return -1;
2880
2881	if (bmap->n_eq == 0)
2882	return 0;
2883
2884	total = isl_basic_map_dim(bmap, type: isl_dim_all);
2885	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
2886	if (total < 0 || n_div < 0)
2887	return -1;
2888	total -= n_div;
2889
2890	for (i = 0; i < bmap->n_eq; ++i)
2891	if (isl_seq_first_non_zero(p: bmap->eq[i] + 1 + total,
2892	len: n_div) == -1)
2893	return i;
2894
2895	return bmap->n_eq;
2896	}
2897
2898	/* Construct a basic map in "space" defined by the equality constraints in "eq".
2899	* The constraints are assumed not to involve any local variables.
2900	*/
2901	static __isl_give isl_basic_map *basic_map_from_equalities(
2902	__isl_take isl_space *space, __isl_take isl_mat *eq)
2903	{
2904	int i, k;
2905	isl_size total;
2906	isl_basic_map *bmap = NULL;
2907
2908	total = isl_space_dim(space, type: isl_dim_all);
2909	if (total < 0 || !eq)
2910	goto error;
2911
2912	if (1 + total != eq->n_col)
2913	isl_die(isl_space_get_ctx(space), isl_error_internal,
2914	"unexpected number of columns", goto error);
2915
2916	bmap = isl_basic_map_alloc_space(space: isl_space_copy(space),
2917	extra: 0, n_eq: eq->n_row, n_ineq: 0);
2918	for (i = 0; i < eq->n_row; ++i) {
2919	k = isl_basic_map_alloc_equality(bmap);
2920	if (k < 0)
2921	goto error;
2922	isl_seq_cpy(dst: bmap->eq[k], src: eq->row[i], len: eq->n_col);
2923	}
2924
2925	isl_space_free(space);
2926	isl_mat_free(mat: eq);
2927	return bmap;
2928	error:
2929	isl_space_free(space);
2930	isl_mat_free(mat: eq);
2931	isl_basic_map_free(bmap);
2932	return NULL;
2933	}
2934
2935	/* Construct and return a variable compression based on the equality
2936	* constraints in "bmap1" and "bmap2" that do not involve the local variables.
2937	* "n1" is the number of (initial) equality constraints in "bmap1"
2938	* that do involve local variables.
2939	* "n2" is the number of (initial) equality constraints in "bmap2"
2940	* that do involve local variables.
2941	* "total" is the total number of other variables.
2942	* This function assumes that Gaussian elimination
2943	* has been applied to the equality constraints in both "bmap1" and "bmap2"
2944	* such that the equality constraints not involving local variables
2945	* are those that start at "n1" or "n2".
2946	*
2947	* If either of "bmap1" and "bmap2" does not have such equality constraints,
2948	* then simply compute the compression based on the equality constraints
2949	* in the other basic map.
2950	* Otherwise, combine the equality constraints from both into a new
2951	* basic map such that Gaussian elimination can be applied to this combination
2952	* and then construct a variable compression from the resulting
2953	* equality constraints.
2954	*/
2955	static __isl_give isl_mat *combined_variable_compression(
2956	__isl_keep isl_basic_map *bmap1, int n1,
2957	__isl_keep isl_basic_map *bmap2, int n2, int total)
2958	{
2959	isl_ctx *ctx;
2960	isl_mat *E1, *E2, *V;
2961	isl_basic_map *bmap;
2962
2963	ctx = isl_basic_map_get_ctx(bmap: bmap1);
2964	if (bmap1->n_eq == n1) {
2965	E2 = isl_mat_sub_alloc6(ctx, row: bmap2->eq,
2966	first_row: n2, n_row: bmap2->n_eq - n2, first_col: 0, n_col: 1 + total);
2967	return isl_mat_variable_compression(B: E2, NULL);
2968	}
2969	if (bmap2->n_eq == n2) {
2970	E1 = isl_mat_sub_alloc6(ctx, row: bmap1->eq,
2971	first_row: n1, n_row: bmap1->n_eq - n1, first_col: 0, n_col: 1 + total);
2972	return isl_mat_variable_compression(B: E1, NULL);
2973	}
2974	E1 = isl_mat_sub_alloc6(ctx, row: bmap1->eq,
2975	first_row: n1, n_row: bmap1->n_eq - n1, first_col: 0, n_col: 1 + total);
2976	E2 = isl_mat_sub_alloc6(ctx, row: bmap2->eq,
2977	first_row: n2, n_row: bmap2->n_eq - n2, first_col: 0, n_col: 1 + total);
2978	E1 = isl_mat_concat(top: E1, bot: E2);
2979	bmap = basic_map_from_equalities(space: isl_basic_map_get_space(bmap: bmap1), eq: E1);
2980	bmap = isl_basic_map_gauss(bmap, NULL);
2981	if (!bmap)
2982	return NULL;
2983	E1 = isl_mat_sub_alloc6(ctx, row: bmap->eq, first_row: 0, n_row: bmap->n_eq, first_col: 0, n_col: 1 + total);
2984	V = isl_mat_variable_compression(B: E1, NULL);
2985	isl_basic_map_free(bmap);
2986
2987	return V;
2988	}
2989
2990	/* Extract the stride constraints from "bmap", compressed
2991	* with respect to both the stride constraints in "context" and
2992	* the remaining equality constraints in both "bmap" and "context".
2993	* "bmap_n_eq" is the number of (initial) stride constraints in "bmap".
2994	* "context_n_eq" is the number of (initial) stride constraints in "context".
2995	*
2996	* Let x be all variables in "bmap" (and "context") other than the local
2997	* variables. First compute a variable compression
2998	*
2999	* x = V x'
3000	*
3001	* based on the non-stride equality constraints in "bmap" and "context".
3002	* Consider the stride constraints of "context",
3003	*
3004	* A(x) + B(y) = 0
3005	*
3006	* with y the local variables and plug in the variable compression,
3007	* resulting in
3008	*
3009	* A(V x') + B(y) = 0
3010	*
3011	* Use these constraints to compute a parameter compression on x'
3012	*
3013	* x' = T x''
3014	*
3015	* Now consider the stride constraints of "bmap"
3016	*
3017	* C(x) + D(y) = 0
3018	*
3019	* and plug in x = V*T x''.
3020	* That is, return A = [C*V*T D].
3021	*/
3022	static __isl_give isl_mat *extract_compressed_stride_constraints(
3023	__isl_keep isl_basic_map *bmap, int bmap_n_eq,
3024	__isl_keep isl_basic_map *context, int context_n_eq)
3025	{
3026	isl_size total, n_div;
3027	isl_ctx *ctx;
3028	isl_mat *A, *B, *T, *V;
3029
3030	total = isl_basic_map_dim(bmap: context, type: isl_dim_all);
3031	n_div = isl_basic_map_dim(bmap: context, type: isl_dim_div);
3032	if (total < 0 || n_div < 0)
3033	return NULL;
3034	total -= n_div;
3035
3036	ctx = isl_basic_map_get_ctx(bmap);
3037
3038	V = combined_variable_compression(bmap1: bmap, n1: bmap_n_eq,
3039	bmap2: context, n2: context_n_eq, total);
3040
3041	A = isl_mat_sub_alloc6(ctx, row: context->eq, first_row: 0, n_row: context_n_eq, first_col: 0, n_col: 1 + total);
3042	B = isl_mat_sub_alloc6(ctx, row: context->eq,
3043	first_row: 0, n_row: context_n_eq, first_col: 1 + total, n_col: n_div);
3044	A = isl_mat_product(left: A, right: isl_mat_copy(mat: V));
3045	T = isl_mat_parameter_compression_ext(B: A, A: B);
3046	T = isl_mat_product(left: V, right: T);
3047
3048	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
3049	if (n_div < 0)
3050	T = isl_mat_free(mat: T);
3051	else
3052	T = isl_mat_diagonal(mat1: T, mat2: isl_mat_identity(ctx, n_row: n_div));
3053
3054	A = isl_mat_sub_alloc6(ctx, row: bmap->eq,
3055	first_row: 0, n_row: bmap_n_eq, first_col: 0, n_col: 1 + total + n_div);
3056	A = isl_mat_product(left: A, right: T);
3057
3058	return A;
3059	}
3060
3061	/* Remove the prime factors from *g that have an exponent that
3062	* is strictly smaller than the exponent in "c".
3063	* All exponents in *g are known to be smaller than or equal
3064	* to those in "c".
3065	*
3066	* That is, if *g is equal to
3067	*
3068	* p_1^{e_1} p_2^{e_2} ... p_n^{e_n}
3069	*
3070	* and "c" is equal to
3071	*
3072	* p_1^{f_1} p_2^{f_2} ... p_n^{f_n}
3073	*
3074	* then update *g to
3075	*
3076	* p_1^{e_1 * (e_1 = f_1)} p_2^{e_2 * (e_2 = f_2)} ...
3077	* p_n^{e_n * (e_n = f_n)}
3078	*
3079	* If e_i = f_i, then c / *g does not have any p_i factors and therefore
3080	* neither does the gcd of *g and c / *g.
3081	* If e_i < f_i, then the gcd of *g and c / *g has a positive
3082	* power min(e_i, s_i) of p_i with s_i = f_i - e_i among its factors.
3083	* Dividing *g by this gcd therefore strictly reduces the exponent
3084	* of the prime factors that need to be removed, while leaving the
3085	* other prime factors untouched.
3086	* Repeating this process until gcd(*g, c / *g) = 1 therefore
3087	* removes all undesired factors, without removing any others.
3088	*/
3089	static void remove_incomplete_powers(isl_int *g, isl_int c)
3090	{
3091	isl_int t;
3092
3093	isl_int_init(t);
3094	for (;;) {
3095	isl_int_divexact(t, c, *g);
3096	isl_int_gcd(t, t, *g);
3097	if (isl_int_is_one(t))
3098	break;
3099	isl_int_divexact(*g, *g, t);
3100	}
3101	isl_int_clear(t);
3102	}
3103
3104	/* Reduce the "n" stride constraints in "bmap" based on a copy "A"
3105	* of the same stride constraints in a compressed space that exploits
3106	* all equalities in the context and the other equalities in "bmap".
3107	*
3108	* If the stride constraints of "bmap" are of the form
3109	*
3110	* C(x) + D(y) = 0
3111	*
3112	* then A is of the form
3113	*
3114	* B(x') + D(y) = 0
3115	*
3116	* If any of these constraints involves only a single local variable y,
3117	* then the constraint appears as
3118	*
3119	* f(x) + m y_i = 0
3120	*
3121	* in "bmap" and as
3122	*
3123	* h(x') + m y_i = 0
3124	*
3125	* in "A".
3126	*
3127	* Let g be the gcd of m and the coefficients of h.
3128	* Then, in particular, g is a divisor of the coefficients of h and
3129	*
3130	* f(x) = h(x')
3131	*
3132	* is known to be a multiple of g.
3133	* If some prime factor in m appears with the same exponent in g,
3134	* then it can be removed from m because f(x) is already known
3135	* to be a multiple of g and therefore in particular of this power
3136	* of the prime factors.
3137	* Prime factors that appear with a smaller exponent in g cannot
3138	* be removed from m.
3139	* Let g' be the divisor of g containing all prime factors that
3140	* appear with the same exponent in m and g, then
3141	*
3142	* f(x) + m y_i = 0
3143	*
3144	* can be replaced by
3145	*
3146	* f(x) + m/g' y_i' = 0
3147	*
3148	* Note that (if g' != 1) this changes the explicit representation
3149	* of y_i to that of y_i', so the integer division at position i
3150	* is marked unknown and later recomputed by a call to
3151	* isl_basic_map_gauss.
3152	*/
3153	static __isl_give isl_basic_map *reduce_stride_constraints(
3154	__isl_take isl_basic_map *bmap, int n, __isl_keep isl_mat *A)
3155	{
3156	int i;
3157	isl_size total, n_div;
3158	int any = 0;
3159	isl_int gcd;
3160
3161	total = isl_basic_map_dim(bmap, type: isl_dim_all);
3162	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
3163	if (total < 0 || n_div < 0 || !A)
3164	return isl_basic_map_free(bmap);
3165	total -= n_div;
3166
3167	isl_int_init(gcd);
3168	for (i = 0; i < n; ++i) {
3169	int div;
3170
3171	div = isl_seq_first_non_zero(p: bmap->eq[i] + 1 + total, len: n_div);
3172	if (div < 0)
3173	isl_die(isl_basic_map_get_ctx(bmap), isl_error_internal,
3174	"equality constraints modified unexpectedly",
3175	goto error);
3176	if (isl_seq_first_non_zero(p: bmap->eq[i] + 1 + total + div + 1,
3177	len: n_div - div - 1) != -1)
3178	continue;
3179	if (isl_mat_row_gcd(mat: A, row: i, gcd: &gcd) < 0)
3180	goto error;
3181	if (isl_int_is_one(gcd))
3182	continue;
3183	remove_incomplete_powers(g: &gcd, c: bmap->eq[i][1 + total + div]);
3184	if (isl_int_is_one(gcd))
3185	continue;
3186	isl_int_divexact(bmap->eq[i][1 + total + div],
3187	bmap->eq[i][1 + total + div], gcd);
3188	bmap = isl_basic_map_mark_div_unknown(bmap, div);
3189	if (!bmap)
3190	goto error;
3191	any = 1;
3192	}
3193	isl_int_clear(gcd);
3194
3195	if (any)
3196	bmap = isl_basic_map_gauss(bmap, NULL);
3197
3198	return bmap;
3199	error:
3200	isl_int_clear(gcd);
3201	isl_basic_map_free(bmap);
3202	return NULL;
3203	}
3204
3205	/* Simplify the stride constraints in "bmap" based on
3206	* the remaining equality constraints in "bmap" and all equality
3207	* constraints in "context".
3208	* Only do this if both "bmap" and "context" have stride constraints.
3209	*
3210	* First extract a copy of the stride constraints in "bmap" in a compressed
3211	* space exploiting all the other equality constraints and then
3212	* use this compressed copy to simplify the original stride constraints.
3213	*/
3214	static __isl_give isl_basic_map *gist_strides(__isl_take isl_basic_map *bmap,
3215	__isl_keep isl_basic_map *context)
3216	{
3217	int bmap_n_eq, context_n_eq;
3218	isl_mat *A;
3219
3220	if (!bmap || !context)
3221	return isl_basic_map_free(bmap);
3222
3223	bmap_n_eq = n_div_eq(bmap);
3224	context_n_eq = n_div_eq(bmap: context);
3225
3226	if (bmap_n_eq < 0 || context_n_eq < 0)
3227	return isl_basic_map_free(bmap);
3228	if (bmap_n_eq == 0 || context_n_eq == 0)
3229	return bmap;
3230
3231	A = extract_compressed_stride_constraints(bmap, bmap_n_eq,
3232	context, context_n_eq);
3233	bmap = reduce_stride_constraints(bmap, n: bmap_n_eq, A);
3234
3235	isl_mat_free(mat: A);
3236
3237	return bmap;
3238	}
3239
3240	/* Return a basic map that has the same intersection with "context" as "bmap"
3241	* and that is as "simple" as possible.
3242	*
3243	* The core computation is performed on the pure constraints.
3244	* When we add back the meaning of the integer divisions, we need
3245	* to (re)introduce the div constraints. If we happen to have
3246	* discovered that some of these integer divisions are equal to
3247	* some affine combination of other variables, then these div
3248	* constraints may end up getting simplified in terms of the equalities,
3249	* resulting in extra inequalities on the other variables that
3250	* may have been removed already or that may not even have been
3251	* part of the input. We try and remove those constraints of
3252	* this form that are most obviously redundant with respect to
3253	* the context. We also remove those div constraints that are
3254	* redundant with respect to the other constraints in the result.
3255	*
3256	* The stride constraints among the equality constraints in "bmap" are
3257	* also simplified with respecting to the other equality constraints
3258	* in "bmap" and with respect to all equality constraints in "context".
3259	*/
3260	__isl_give isl_basic_map *isl_basic_map_gist(__isl_take isl_basic_map *bmap,
3261	__isl_take isl_basic_map *context)
3262	{
3263	isl_basic_set *bset, *eq;
3264	isl_basic_map *eq_bmap;
3265	isl_size total, n_div, n_div_bmap;
3266	unsigned extra, n_eq, n_ineq;
3267
3268	if (!bmap || !context)
3269	goto error;
3270
3271	if (isl_basic_map_plain_is_universe(bmap)) {
3272	isl_basic_map_free(bmap: context);
3273	return bmap;
3274	}
3275	if (isl_basic_map_plain_is_empty(bmap: context)) {
3276	isl_space *space = isl_basic_map_get_space(bmap);
3277	isl_basic_map_free(bmap);
3278	isl_basic_map_free(bmap: context);
3279	return isl_basic_map_universe(space);
3280	}
3281	if (isl_basic_map_plain_is_empty(bmap)) {
3282	isl_basic_map_free(bmap: context);
3283	return bmap;
3284	}
3285
3286	bmap = isl_basic_map_remove_redundancies(bmap);
3287	context = isl_basic_map_remove_redundancies(bmap: context);
3288	bmap = isl_basic_map_order_divs(bmap);
3289	context = isl_basic_map_align_divs(dst: context, src: bmap);
3290
3291	n_div = isl_basic_map_dim(bmap: context, type: isl_dim_div);
3292	total = isl_basic_map_dim(bmap, type: isl_dim_all);
3293	n_div_bmap = isl_basic_map_dim(bmap, type: isl_dim_div);
3294	if (n_div < 0 || total < 0 || n_div_bmap < 0)
3295	goto error;
3296	extra = n_div - n_div_bmap;
3297
3298	bset = isl_basic_map_underlying_set(bmap: isl_basic_map_copy(bmap));
3299	bset = isl_basic_set_add_dims(bset, type: isl_dim_set, n: extra);
3300	bset = uset_gist(bset,
3301	context: isl_basic_map_underlying_set(bmap: isl_basic_map_copy(bmap: context)));
3302	bset = isl_basic_set_project_out(bset, type: isl_dim_set, first: total, n: extra);
3303
3304	if (!bset || bset->n_eq == 0 || n_div == 0 ||
3305	isl_basic_set_plain_is_empty(bset)) {
3306	isl_basic_map_free(bmap: context);
3307	return isl_basic_map_overlying_set(bset, like: bmap);
3308	}
3309
3310	n_eq = bset->n_eq;
3311	n_ineq = bset->n_ineq;
3312	eq = isl_basic_set_copy(bset);
3313	eq = isl_basic_set_cow(bset: eq);
3314	eq = isl_basic_set_free_inequality(bset: eq, n: n_ineq);
3315	bset = isl_basic_set_free_equality(bset, n: n_eq);
3316
3317	eq_bmap = isl_basic_map_overlying_set(bset: eq, like: isl_basic_map_copy(bmap));
3318	eq_bmap = gist_strides(bmap: eq_bmap, context);
3319	eq_bmap = isl_basic_map_remove_shifted_constraints(bmap: eq_bmap, context);
3320	bmap = isl_basic_map_overlying_set(bset, like: bmap);
3321	bmap = isl_basic_map_intersect(bmap1: bmap, bmap2: eq_bmap);
3322	bmap = isl_basic_map_remove_redundancies(bmap);
3323
3324	return bmap;
3325	error:
3326	isl_basic_map_free(bmap);
3327	isl_basic_map_free(bmap: context);
3328	return NULL;
3329	}
3330
3331	/*
3332	* Assumes context has no implicit divs.
3333	*/
3334	__isl_give isl_map *isl_map_gist_basic_map(__isl_take isl_map *map,
3335	__isl_take isl_basic_map *context)
3336	{
3337	int i;
3338
3339	if (!map || !context)
3340	goto error;
3341
3342	if (isl_basic_map_plain_is_empty(bmap: context)) {
3343	isl_space *space = isl_map_get_space(map);
3344	isl_map_free(map);
3345	isl_basic_map_free(bmap: context);
3346	return isl_map_universe(space);
3347	}
3348
3349	context = isl_basic_map_remove_redundancies(bmap: context);
3350	map = isl_map_cow(map);
3351	if (isl_map_basic_map_check_equal_space(map, bmap: context) < 0)
3352	goto error;
3353	map = isl_map_compute_divs(map);
3354	if (!map)
3355	goto error;
3356	for (i = map->n - 1; i >= 0; --i) {
3357	map->p[i] = isl_basic_map_gist(bmap: map->p[i],
3358	context: isl_basic_map_copy(bmap: context));
3359	if (!map->p[i])
3360	goto error;
3361	if (isl_basic_map_plain_is_empty(bmap: map->p[i])) {
3362	isl_basic_map_free(bmap: map->p[i]);
3363	if (i != map->n - 1)
3364	map->p[i] = map->p[map->n - 1];
3365	map->n--;
3366	}
3367	}
3368	isl_basic_map_free(bmap: context);
3369	ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3370	return map;
3371	error:
3372	isl_map_free(map);
3373	isl_basic_map_free(bmap: context);
3374	return NULL;
3375	}
3376
3377	/* Drop all inequalities from "bmap" that also appear in "context".
3378	* "context" is assumed to have only known local variables and
3379	* the initial local variables of "bmap" are assumed to be the same
3380	* as those of "context".
3381	* The constraints of both "bmap" and "context" are assumed
3382	* to have been sorted using isl_basic_map_sort_constraints.
3383	*
3384	* Run through the inequality constraints of "bmap" and "context"
3385	* in sorted order.
3386	* If a constraint of "bmap" involves variables not in "context",
3387	* then it cannot appear in "context".
3388	* If a matching constraint is found, it is removed from "bmap".
3389	*/
3390	static __isl_give isl_basic_map *drop_inequalities(
3391	__isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3392	{
3393	int i1, i2;
3394	isl_size total, bmap_total;
3395	unsigned extra;
3396
3397	total = isl_basic_map_dim(bmap: context, type: isl_dim_all);
3398	bmap_total = isl_basic_map_dim(bmap, type: isl_dim_all);
3399	if (total < 0 || bmap_total < 0)
3400	return isl_basic_map_free(bmap);
3401
3402	extra = bmap_total - total;
3403
3404	i1 = bmap->n_ineq - 1;
3405	i2 = context->n_ineq - 1;
3406	while (bmap && i1 >= 0 && i2 >= 0) {
3407	int cmp;
3408
3409	if (isl_seq_first_non_zero(p: bmap->ineq[i1] + 1 + total,
3410	len: extra) != -1) {
3411	--i1;
3412	continue;
3413	}
3414	cmp = isl_basic_map_constraint_cmp(bmap: context, c1: bmap->ineq[i1],
3415	c2: context->ineq[i2]);
3416	if (cmp < 0) {
3417	--i2;
3418	continue;
3419	}
3420	if (cmp > 0) {
3421	--i1;
3422	continue;
3423	}
3424	if (isl_int_eq(bmap->ineq[i1][0], context->ineq[i2][0])) {
3425	bmap = isl_basic_map_cow(bmap);
3426	if (isl_basic_map_drop_inequality(bmap, pos: i1) < 0)
3427	bmap = isl_basic_map_free(bmap);
3428	}
3429	--i1;
3430	--i2;
3431	}
3432
3433	return bmap;
3434	}
3435
3436	/* Drop all equalities from "bmap" that also appear in "context".
3437	* "context" is assumed to have only known local variables and
3438	* the initial local variables of "bmap" are assumed to be the same
3439	* as those of "context".
3440	*
3441	* Run through the equality constraints of "bmap" and "context"
3442	* in sorted order.
3443	* If a constraint of "bmap" involves variables not in "context",
3444	* then it cannot appear in "context".
3445	* If a matching constraint is found, it is removed from "bmap".
3446	*/
3447	static __isl_give isl_basic_map *drop_equalities(
3448	__isl_take isl_basic_map *bmap, __isl_keep isl_basic_map *context)
3449	{
3450	int i1, i2;
3451	isl_size total, bmap_total;
3452	unsigned extra;
3453
3454	total = isl_basic_map_dim(bmap: context, type: isl_dim_all);
3455	bmap_total = isl_basic_map_dim(bmap, type: isl_dim_all);
3456	if (total < 0 || bmap_total < 0)
3457	return isl_basic_map_free(bmap);
3458
3459	extra = bmap_total - total;
3460
3461	i1 = bmap->n_eq - 1;
3462	i2 = context->n_eq - 1;
3463
3464	while (bmap && i1 >= 0 && i2 >= 0) {
3465	int last1, last2;
3466
3467	if (isl_seq_first_non_zero(p: bmap->eq[i1] + 1 + total,
3468	len: extra) != -1)
3469	break;
3470	last1 = isl_seq_last_non_zero(p: bmap->eq[i1] + 1, len: total);
3471	last2 = isl_seq_last_non_zero(p: context->eq[i2] + 1, len: total);
3472	if (last1 > last2) {
3473	--i2;
3474	continue;
3475	}
3476	if (last1 < last2) {
3477	--i1;
3478	continue;
3479	}
3480	if (isl_seq_eq(p1: bmap->eq[i1], p2: context->eq[i2], len: 1 + total)) {
3481	bmap = isl_basic_map_cow(bmap);
3482	if (isl_basic_map_drop_equality(bmap, pos: i1) < 0)
3483	bmap = isl_basic_map_free(bmap);
3484	}
3485	--i1;
3486	--i2;
3487	}
3488
3489	return bmap;
3490	}
3491
3492	/* Remove the constraints in "context" from "bmap".
3493	* "context" is assumed to have explicit representations
3494	* for all local variables.
3495	*
3496	* First align the divs of "bmap" to those of "context" and
3497	* sort the constraints. Then drop all constraints from "bmap"
3498	* that appear in "context".
3499	*/
3500	__isl_give isl_basic_map *isl_basic_map_plain_gist(
3501	__isl_take isl_basic_map *bmap, __isl_take isl_basic_map *context)
3502	{
3503	isl_bool done, known;
3504
3505	done = isl_basic_map_plain_is_universe(bmap: context);
3506	if (done == isl_bool_false)
3507	done = isl_basic_map_plain_is_universe(bmap);
3508	if (done == isl_bool_false)
3509	done = isl_basic_map_plain_is_empty(bmap: context);
3510	if (done == isl_bool_false)
3511	done = isl_basic_map_plain_is_empty(bmap);
3512	if (done < 0)
3513	goto error;
3514	if (done) {
3515	isl_basic_map_free(bmap: context);
3516	return bmap;
3517	}
3518	known = isl_basic_map_divs_known(bmap: context);
3519	if (known < 0)
3520	goto error;
3521	if (!known)
3522	isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3523	"context has unknown divs", goto error);
3524
3525	context = isl_basic_map_order_divs(bmap: context);
3526	bmap = isl_basic_map_align_divs(dst: bmap, src: context);
3527	bmap = isl_basic_map_gauss(bmap, NULL);
3528	bmap = isl_basic_map_sort_constraints(bmap);
3529	context = isl_basic_map_sort_constraints(bmap: context);
3530
3531	bmap = drop_inequalities(bmap, context);
3532	bmap = drop_equalities(bmap, context);
3533
3534	isl_basic_map_free(bmap: context);
3535	bmap = isl_basic_map_finalize(bmap);
3536	return bmap;
3537	error:
3538	isl_basic_map_free(bmap);
3539	isl_basic_map_free(bmap: context);
3540	return NULL;
3541	}
3542
3543	/* Replace "map" by the disjunct at position "pos" and free "context".
3544	*/
3545	static __isl_give isl_map *replace_by_disjunct(__isl_take isl_map *map,
3546	int pos, __isl_take isl_basic_map *context)
3547	{
3548	isl_basic_map *bmap;
3549
3550	bmap = isl_basic_map_copy(bmap: map->p[pos]);
3551	isl_map_free(map);
3552	isl_basic_map_free(bmap: context);
3553	return isl_map_from_basic_map(bmap);
3554	}
3555
3556	/* Remove the constraints in "context" from "map".
3557	* If any of the disjuncts in the result turns out to be the universe,
3558	* then return this universe.
3559	* "context" is assumed to have explicit representations
3560	* for all local variables.
3561	*/
3562	__isl_give isl_map *isl_map_plain_gist_basic_map(__isl_take isl_map *map,
3563	__isl_take isl_basic_map *context)
3564	{
3565	int i;
3566	isl_bool univ, known;
3567
3568	univ = isl_basic_map_plain_is_universe(bmap: context);
3569	if (univ < 0)
3570	goto error;
3571	if (univ) {
3572	isl_basic_map_free(bmap: context);
3573	return map;
3574	}
3575	known = isl_basic_map_divs_known(bmap: context);
3576	if (known < 0)
3577	goto error;
3578	if (!known)
3579	isl_die(isl_map_get_ctx(map), isl_error_invalid,
3580	"context has unknown divs", goto error);
3581
3582	map = isl_map_cow(map);
3583	if (!map)
3584	goto error;
3585	for (i = 0; i < map->n; ++i) {
3586	map->p[i] = isl_basic_map_plain_gist(bmap: map->p[i],
3587	context: isl_basic_map_copy(bmap: context));
3588	univ = isl_basic_map_plain_is_universe(bmap: map->p[i]);
3589	if (univ < 0)
3590	goto error;
3591	if (univ && map->n > 1)
3592	return replace_by_disjunct(map, pos: i, context);
3593	}
3594
3595	isl_basic_map_free(bmap: context);
3596	ISL_F_CLR(map, ISL_MAP_NORMALIZED);
3597	if (map->n > 1)
3598	ISL_F_CLR(map, ISL_MAP_DISJOINT);
3599	return map;
3600	error:
3601	isl_map_free(map);
3602	isl_basic_map_free(bmap: context);
3603	return NULL;
3604	}
3605
3606	/* Remove the constraints in "context" from "set".
3607	* If any of the disjuncts in the result turns out to be the universe,
3608	* then return this universe.
3609	* "context" is assumed to have explicit representations
3610	* for all local variables.
3611	*/
3612	__isl_give isl_set *isl_set_plain_gist_basic_set(__isl_take isl_set *set,
3613	__isl_take isl_basic_set *context)
3614	{
3615	return set_from_map(isl_map_plain_gist_basic_map(map: set_to_map(set),
3616	context: bset_to_bmap(bset: context)));
3617	}
3618
3619	/* Remove the constraints in "context" from "map".
3620	* If any of the disjuncts in the result turns out to be the universe,
3621	* then return this universe.
3622	* "context" is assumed to consist of a single disjunct and
3623	* to have explicit representations for all local variables.
3624	*/
3625	__isl_give isl_map *isl_map_plain_gist(__isl_take isl_map *map,
3626	__isl_take isl_map *context)
3627	{
3628	isl_basic_map *hull;
3629
3630	hull = isl_map_unshifted_simple_hull(map: context);
3631	return isl_map_plain_gist_basic_map(map, context: hull);
3632	}
3633
3634	/* Replace "map" by a universe map in the same space and free "drop".
3635	*/
3636	static __isl_give isl_map *replace_by_universe(__isl_take isl_map *map,
3637	__isl_take isl_map *drop)
3638	{
3639	isl_map *res;
3640
3641	res = isl_map_universe(space: isl_map_get_space(map));
3642	isl_map_free(map);
3643	isl_map_free(map: drop);
3644	return res;
3645	}
3646
3647	/* Return a map that has the same intersection with "context" as "map"
3648	* and that is as "simple" as possible.
3649	*
3650	* If "map" is already the universe, then we cannot make it any simpler.
3651	* Similarly, if "context" is the universe, then we cannot exploit it
3652	* to simplify "map"
3653	* If "map" and "context" are identical to each other, then we can
3654	* return the corresponding universe.
3655	*
3656	* If either "map" or "context" consists of multiple disjuncts,
3657	* then check if "context" happens to be a subset of "map",
3658	* in which case all constraints can be removed.
3659	* In case of multiple disjuncts, the standard procedure
3660	* may not be able to detect that all constraints can be removed.
3661	*
3662	* If none of these cases apply, we have to work a bit harder.
3663	* During this computation, we make use of a single disjunct context,
3664	* so if the original context consists of more than one disjunct
3665	* then we need to approximate the context by a single disjunct set.
3666	* Simply taking the simple hull may drop constraints that are
3667	* only implicitly available in each disjunct. We therefore also
3668	* look for constraints among those defining "map" that are valid
3669	* for the context. These can then be used to simplify away
3670	* the corresponding constraints in "map".
3671	*/
3672	__isl_give isl_map *isl_map_gist(__isl_take isl_map *map,
3673	__isl_take isl_map *context)
3674	{
3675	int equal;
3676	int is_universe;
3677	isl_size n_disjunct_map, n_disjunct_context;
3678	isl_bool subset;
3679	isl_basic_map *hull;
3680
3681	is_universe = isl_map_plain_is_universe(map);
3682	if (is_universe >= 0 && !is_universe)
3683	is_universe = isl_map_plain_is_universe(map: context);
3684	if (is_universe < 0)
3685	goto error;
3686	if (is_universe) {
3687	isl_map_free(map: context);
3688	return map;
3689	}
3690
3691	isl_map_align_params_bin(map1: &map, map2: &context);
3692	equal = isl_map_plain_is_equal(map1: map, map2: context);
3693	if (equal < 0)
3694	goto error;
3695	if (equal)
3696	return replace_by_universe(map, drop: context);
3697
3698	n_disjunct_map = isl_map_n_basic_map(map);
3699	n_disjunct_context = isl_map_n_basic_map(map: context);
3700	if (n_disjunct_map < 0 || n_disjunct_context < 0)
3701	goto error;
3702	if (n_disjunct_map != 1 || n_disjunct_context != 1) {
3703	subset = isl_map_is_subset(map1: context, map2: map);
3704	if (subset < 0)
3705	goto error;
3706	if (subset)
3707	return replace_by_universe(map, drop: context);
3708	}
3709
3710	context = isl_map_compute_divs(map: context);
3711	if (!context)
3712	goto error;
3713	if (n_disjunct_context == 1) {
3714	hull = isl_map_simple_hull(map: context);
3715	} else {
3716	isl_ctx *ctx;
3717	isl_map_list *list;
3718
3719	ctx = isl_map_get_ctx(map);
3720	list = isl_map_list_alloc(ctx, n: 2);
3721	list = isl_map_list_add(list, el: isl_map_copy(map: context));
3722	list = isl_map_list_add(list, el: isl_map_copy(map));
3723	hull = isl_map_unshifted_simple_hull_from_map_list(map: context,
3724	list);
3725	}
3726	return isl_map_gist_basic_map(map, context: hull);
3727	error:
3728	isl_map_free(map);
3729	isl_map_free(map: context);
3730	return NULL;
3731	}
3732
3733	__isl_give isl_basic_set *isl_basic_set_gist(__isl_take isl_basic_set *bset,
3734	__isl_take isl_basic_set *context)
3735	{
3736	return bset_from_bmap(bmap: isl_basic_map_gist(bmap: bset_to_bmap(bset),
3737	context: bset_to_bmap(bset: context)));
3738	}
3739
3740	__isl_give isl_set *isl_set_gist_basic_set(__isl_take isl_set *set,
3741	__isl_take isl_basic_set *context)
3742	{
3743	return set_from_map(isl_map_gist_basic_map(map: set_to_map(set),
3744	context: bset_to_bmap(bset: context)));
3745	}
3746
3747	__isl_give isl_set *isl_set_gist_params_basic_set(__isl_take isl_set *set,
3748	__isl_take isl_basic_set *context)
3749	{
3750	isl_space *space = isl_set_get_space(set);
3751	isl_basic_set *dom_context = isl_basic_set_universe(space);
3752	dom_context = isl_basic_set_intersect_params(bset1: dom_context, bset2: context);
3753	return isl_set_gist_basic_set(set, context: dom_context);
3754	}
3755
3756	__isl_give isl_set *isl_set_gist(__isl_take isl_set *set,
3757	__isl_take isl_set *context)
3758	{
3759	return set_from_map(isl_map_gist(map: set_to_map(set), context: set_to_map(context)));
3760	}
3761
3762	/* Compute the gist of "bmap" with respect to the constraints "context"
3763	* on the domain.
3764	*/
3765	__isl_give isl_basic_map *isl_basic_map_gist_domain(
3766	__isl_take isl_basic_map *bmap, __isl_take isl_basic_set *context)
3767	{
3768	isl_space *space = isl_basic_map_get_space(bmap);
3769	isl_basic_map *bmap_context = isl_basic_map_universe(space);
3770
3771	bmap_context = isl_basic_map_intersect_domain(bmap: bmap_context, bset: context);
3772	return isl_basic_map_gist(bmap, context: bmap_context);
3773	}
3774
3775	__isl_give isl_map *isl_map_gist_domain(__isl_take isl_map *map,
3776	__isl_take isl_set *context)
3777	{
3778	isl_map *map_context = isl_map_universe(space: isl_map_get_space(map));
3779	map_context = isl_map_intersect_domain(map: map_context, set: context);
3780	return isl_map_gist(map, context: map_context);
3781	}
3782
3783	__isl_give isl_map *isl_map_gist_range(__isl_take isl_map *map,
3784	__isl_take isl_set *context)
3785	{
3786	isl_map *map_context = isl_map_universe(space: isl_map_get_space(map));
3787	map_context = isl_map_intersect_range(map: map_context, set: context);
3788	return isl_map_gist(map, context: map_context);
3789	}
3790
3791	__isl_give isl_map *isl_map_gist_params(__isl_take isl_map *map,
3792	__isl_take isl_set *context)
3793	{
3794	isl_map *map_context = isl_map_universe(space: isl_map_get_space(map));
3795	map_context = isl_map_intersect_params(map: map_context, params: context);
3796	return isl_map_gist(map, context: map_context);
3797	}
3798
3799	__isl_give isl_set *isl_set_gist_params(__isl_take isl_set *set,
3800	__isl_take isl_set *context)
3801	{
3802	return isl_map_gist_params(map: set, context);
3803	}
3804
3805	/* Quick check to see if two basic maps are disjoint.
3806	* In particular, we reduce the equalities and inequalities of
3807	* one basic map in the context of the equalities of the other
3808	* basic map and check if we get a contradiction.
3809	*/
3810	isl_bool isl_basic_map_plain_is_disjoint(__isl_keep isl_basic_map *bmap1,
3811	__isl_keep isl_basic_map *bmap2)
3812	{
3813	struct isl_vec *v = NULL;
3814	int *elim = NULL;
3815	isl_size total;
3816	int i;
3817
3818	if (isl_basic_map_check_equal_space(bmap1, bmap2) < 0)
3819	return isl_bool_error;
3820	if (bmap1->n_div || bmap2->n_div)
3821	return isl_bool_false;
3822	if (!bmap1->n_eq && !bmap2->n_eq)
3823	return isl_bool_false;
3824
3825	total = isl_space_dim(space: bmap1->dim, type: isl_dim_all);
3826	if (total < 0)
3827	return isl_bool_error;
3828	if (total == 0)
3829	return isl_bool_false;
3830	v = isl_vec_alloc(ctx: bmap1->ctx, size: 1 + total);
3831	if (!v)
3832	goto error;
3833	elim = isl_alloc_array(bmap1->ctx, int, total);
3834	if (!elim)
3835	goto error;
3836	compute_elimination_index(bmap: bmap1, elim, len: total);
3837	for (i = 0; i < bmap2->n_eq; ++i) {
3838	int reduced;
3839	reduced = reduced_using_equalities(dst: v->block.data, src: bmap2->eq[i],
3840	bmap: bmap1, elim, total);
3841	if (reduced && !isl_int_is_zero(v->block.data[0]) &&
3842	isl_seq_first_non_zero(p: v->block.data + 1, len: total) == -1)
3843	goto disjoint;
3844	}
3845	for (i = 0; i < bmap2->n_ineq; ++i) {
3846	int reduced;
3847	reduced = reduced_using_equalities(dst: v->block.data,
3848	src: bmap2->ineq[i], bmap: bmap1, elim, total);
3849	if (reduced && isl_int_is_neg(v->block.data[0]) &&
3850	isl_seq_first_non_zero(p: v->block.data + 1, len: total) == -1)
3851	goto disjoint;
3852	}
3853	compute_elimination_index(bmap: bmap2, elim, len: total);
3854	for (i = 0; i < bmap1->n_ineq; ++i) {
3855	int reduced;
3856	reduced = reduced_using_equalities(dst: v->block.data,
3857	src: bmap1->ineq[i], bmap: bmap2, elim, total);
3858	if (reduced && isl_int_is_neg(v->block.data[0]) &&
3859	isl_seq_first_non_zero(p: v->block.data + 1, len: total) == -1)
3860	goto disjoint;
3861	}
3862	isl_vec_free(vec: v);
3863	free(ptr: elim);
3864	return isl_bool_false;
3865	disjoint:
3866	isl_vec_free(vec: v);
3867	free(ptr: elim);
3868	return isl_bool_true;
3869	error:
3870	isl_vec_free(vec: v);
3871	free(ptr: elim);
3872	return isl_bool_error;
3873	}
3874
3875	int isl_basic_set_plain_is_disjoint(__isl_keep isl_basic_set *bset1,
3876	__isl_keep isl_basic_set *bset2)
3877	{
3878	return isl_basic_map_plain_is_disjoint(bmap1: bset_to_bmap(bset: bset1),
3879	bmap2: bset_to_bmap(bset: bset2));
3880	}
3881
3882	/* Does "test" hold for all pairs of basic maps in "map1" and "map2"?
3883	*/
3884	static isl_bool all_pairs(__isl_keep isl_map *map1, __isl_keep isl_map *map2,
3885	isl_bool (*test)(__isl_keep isl_basic_map *bmap1,
3886	__isl_keep isl_basic_map *bmap2))
3887	{
3888	int i, j;
3889
3890	if (!map1 || !map2)
3891	return isl_bool_error;
3892
3893	for (i = 0; i < map1->n; ++i) {
3894	for (j = 0; j < map2->n; ++j) {
3895	isl_bool d = test(map1->p[i], map2->p[j]);
3896	if (d != isl_bool_true)
3897	return d;
3898	}
3899	}
3900
3901	return isl_bool_true;
3902	}
3903
3904	/* Are "map1" and "map2" obviously disjoint, based on information
3905	* that can be derived without looking at the individual basic maps?
3906	*
3907	* In particular, if one of them is empty or if they live in different spaces
3908	* (ignoring parameters), then they are clearly disjoint.
3909	*/
3910	static isl_bool isl_map_plain_is_disjoint_global(__isl_keep isl_map *map1,
3911	__isl_keep isl_map *map2)
3912	{
3913	isl_bool disjoint;
3914	isl_bool match;
3915
3916	if (!map1 || !map2)
3917	return isl_bool_error;
3918
3919	disjoint = isl_map_plain_is_empty(map: map1);
3920	if (disjoint < 0 || disjoint)
3921	return disjoint;
3922
3923	disjoint = isl_map_plain_is_empty(map: map2);
3924	if (disjoint < 0 || disjoint)
3925	return disjoint;
3926
3927	match = isl_map_tuple_is_equal(map1, type1: isl_dim_in, map2, type2: isl_dim_in);
3928	if (match < 0 || !match)
3929	return match < 0 ? isl_bool_error : isl_bool_true;
3930
3931	match = isl_map_tuple_is_equal(map1, type1: isl_dim_out, map2, type2: isl_dim_out);
3932	if (match < 0 || !match)
3933	return match < 0 ? isl_bool_error : isl_bool_true;
3934
3935	return isl_bool_false;
3936	}
3937
3938	/* Are "map1" and "map2" obviously disjoint?
3939	*
3940	* If one of them is empty or if they live in different spaces (ignoring
3941	* parameters), then they are clearly disjoint.
3942	* This is checked by isl_map_plain_is_disjoint_global.
3943	*
3944	* If they have different parameters, then we skip any further tests.
3945	*
3946	* If they are obviously equal, but not obviously empty, then we will
3947	* not be able to detect if they are disjoint.
3948	*
3949	* Otherwise we check if each basic map in "map1" is obviously disjoint
3950	* from each basic map in "map2".
3951	*/
3952	isl_bool isl_map_plain_is_disjoint(__isl_keep isl_map *map1,
3953	__isl_keep isl_map *map2)
3954	{
3955	isl_bool disjoint;
3956	isl_bool intersect;
3957	isl_bool match;
3958
3959	disjoint = isl_map_plain_is_disjoint_global(map1, map2);
3960	if (disjoint < 0 || disjoint)
3961	return disjoint;
3962
3963	match = isl_map_has_equal_params(map1, map2);
3964	if (match < 0 || !match)
3965	return match < 0 ? isl_bool_error : isl_bool_false;
3966
3967	intersect = isl_map_plain_is_equal(map1, map2);
3968	if (intersect < 0 || intersect)
3969	return intersect < 0 ? isl_bool_error : isl_bool_false;
3970
3971	return all_pairs(map1, map2, test: &isl_basic_map_plain_is_disjoint);
3972	}
3973
3974	/* Are "map1" and "map2" disjoint?
3975	* The parameters are assumed to have been aligned.
3976	*
3977	* In particular, check whether all pairs of basic maps are disjoint.
3978	*/
3979	static isl_bool isl_map_is_disjoint_aligned(__isl_keep isl_map *map1,
3980	__isl_keep isl_map *map2)
3981	{
3982	return all_pairs(map1, map2, test: &isl_basic_map_is_disjoint);
3983	}
3984
3985	/* Are "map1" and "map2" disjoint?
3986	*
3987	* They are disjoint if they are "obviously disjoint" or if one of them
3988	* is empty. Otherwise, they are not disjoint if one of them is universal.
3989	* If the two inputs are (obviously) equal and not empty, then they are
3990	* not disjoint.
3991	* If none of these cases apply, then check if all pairs of basic maps
3992	* are disjoint after aligning the parameters.
3993	*/
3994	isl_bool isl_map_is_disjoint(__isl_keep isl_map *map1, __isl_keep isl_map *map2)
3995	{
3996	isl_bool disjoint;
3997	isl_bool intersect;
3998
3999	disjoint = isl_map_plain_is_disjoint_global(map1, map2);
4000	if (disjoint < 0 || disjoint)
4001	return disjoint;
4002
4003	disjoint = isl_map_is_empty(map: map1);
4004	if (disjoint < 0 || disjoint)
4005	return disjoint;
4006
4007	disjoint = isl_map_is_empty(map: map2);
4008	if (disjoint < 0 || disjoint)
4009	return disjoint;
4010
4011	intersect = isl_map_plain_is_universe(map: map1);
4012	if (intersect < 0 || intersect)
4013	return isl_bool_not(b: intersect);
4014
4015	intersect = isl_map_plain_is_universe(map: map2);
4016	if (intersect < 0 || intersect)
4017	return isl_bool_not(b: intersect);
4018
4019	intersect = isl_map_plain_is_equal(map1, map2);
4020	if (intersect < 0 || intersect)
4021	return isl_bool_not(b: intersect);
4022
4023	return isl_map_align_params_map_map_and_test(map1, map2,
4024	fn: &isl_map_is_disjoint_aligned);
4025	}
4026
4027	/* Are "bmap1" and "bmap2" disjoint?
4028	*
4029	* They are disjoint if they are "obviously disjoint" or if one of them
4030	* is empty. Otherwise, they are not disjoint if one of them is universal.
4031	* If none of these cases apply, we compute the intersection and see if
4032	* the result is empty.
4033	*/
4034	isl_bool isl_basic_map_is_disjoint(__isl_keep isl_basic_map *bmap1,
4035	__isl_keep isl_basic_map *bmap2)
4036	{
4037	isl_bool disjoint;
4038	isl_bool intersect;
4039	isl_basic_map *test;
4040
4041	disjoint = isl_basic_map_plain_is_disjoint(bmap1, bmap2);
4042	if (disjoint < 0 || disjoint)
4043	return disjoint;
4044
4045	disjoint = isl_basic_map_is_empty(bmap: bmap1);
4046	if (disjoint < 0 || disjoint)
4047	return disjoint;
4048
4049	disjoint = isl_basic_map_is_empty(bmap: bmap2);
4050	if (disjoint < 0 || disjoint)
4051	return disjoint;
4052
4053	intersect = isl_basic_map_plain_is_universe(bmap: bmap1);
4054	if (intersect < 0 || intersect)
4055	return isl_bool_not(b: intersect);
4056
4057	intersect = isl_basic_map_plain_is_universe(bmap: bmap2);
4058	if (intersect < 0 || intersect)
4059	return isl_bool_not(b: intersect);
4060
4061	test = isl_basic_map_intersect(bmap1: isl_basic_map_copy(bmap: bmap1),
4062	bmap2: isl_basic_map_copy(bmap: bmap2));
4063	disjoint = isl_basic_map_is_empty(bmap: test);
4064	isl_basic_map_free(bmap: test);
4065
4066	return disjoint;
4067	}
4068
4069	/* Are "bset1" and "bset2" disjoint?
4070	*/
4071	isl_bool isl_basic_set_is_disjoint(__isl_keep isl_basic_set *bset1,
4072	__isl_keep isl_basic_set *bset2)
4073	{
4074	return isl_basic_map_is_disjoint(bmap1: bset1, bmap2: bset2);
4075	}
4076
4077	isl_bool isl_set_plain_is_disjoint(__isl_keep isl_set *set1,
4078	__isl_keep isl_set *set2)
4079	{
4080	return isl_map_plain_is_disjoint(map1: set_to_map(set1), map2: set_to_map(set2));
4081	}
4082
4083	/* Are "set1" and "set2" disjoint?
4084	*/
4085	isl_bool isl_set_is_disjoint(__isl_keep isl_set *set1, __isl_keep isl_set *set2)
4086	{
4087	return isl_map_is_disjoint(map1: set1, map2: set2);
4088	}
4089
4090	/* Is "v" equal to 0, 1 or -1?
4091	*/
4092	static int is_zero_or_one(isl_int v)
4093	{
4094	return isl_int_is_zero(v) || isl_int_is_one(v) || isl_int_is_negone(v);
4095	}
4096
4097	/* Are the "n" coefficients starting at "first" of inequality constraints
4098	* "i" and "j" of "bmap" opposite to each other?
4099	*/
4100	static int is_opposite_part(__isl_keep isl_basic_map *bmap, int i, int j,
4101	int first, int n)
4102	{
4103	return isl_seq_is_neg(p1: bmap->ineq[i] + first, p2: bmap->ineq[j] + first, len: n);
4104	}
4105
4106	/* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4107	* apart from the constant term?
4108	*/
4109	static isl_bool is_opposite(__isl_keep isl_basic_map *bmap, int i, int j)
4110	{
4111	isl_size total;
4112
4113	total = isl_basic_map_dim(bmap, type: isl_dim_all);
4114	if (total < 0)
4115	return isl_bool_error;
4116	return is_opposite_part(bmap, i, j, first: 1, n: total);
4117	}
4118
4119	/* Check if we can combine a given div with lower bound l and upper
4120	* bound u with some other div and if so return that other div.
4121	* Otherwise, return a position beyond the integer divisions.
4122	* Return -1 on error.
4123	*
4124	* We first check that
4125	* - the bounds are opposites of each other (except for the constant
4126	* term)
4127	* - the bounds do not reference any other div
4128	* - no div is defined in terms of this div
4129	*
4130	* Let m be the size of the range allowed on the div by the bounds.
4131	* That is, the bounds are of the form
4132	*
4133	* e <= a <= e + m - 1
4134	*
4135	* with e some expression in the other variables.
4136	* We look for another div b such that no third div is defined in terms
4137	* of this second div b and such that in any constraint that contains
4138	* a (except for the given lower and upper bound), also contains b
4139	* with a coefficient that is m times that of b.
4140	* That is, all constraints (except for the lower and upper bound)
4141	* are of the form
4142	*
4143	* e + f (a + m b) >= 0
4144	*
4145	* Furthermore, in the constraints that only contain b, the coefficient
4146	* of b should be equal to 1 or -1.
4147	* If so, we return b so that "a + m b" can be replaced by
4148	* a single div "c = a + m b".
4149	*/
4150	static int div_find_coalesce(__isl_keep isl_basic_map *bmap, int *pairs,
4151	unsigned div, unsigned l, unsigned u)
4152	{
4153	int i, j;
4154	unsigned n_div;
4155	isl_size v_div;
4156	int coalesce;
4157	isl_bool opp;
4158
4159	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4160	if (n_div <= 1)
4161	return n_div;
4162	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
4163	if (v_div < 0)
4164	return -1;
4165	if (isl_seq_first_non_zero(p: bmap->ineq[l] + 1 + v_div, len: div) != -1)
4166	return n_div;
4167	if (isl_seq_first_non_zero(p: bmap->ineq[l] + 1 + v_div + div + 1,
4168	len: n_div - div - 1) != -1)
4169	return n_div;
4170	opp = is_opposite(bmap, i: l, j: u);
4171	if (opp < 0 || !opp)
4172	return opp < 0 ? -1 : n_div;
4173
4174	for (i = 0; i < n_div; ++i) {
4175	if (isl_int_is_zero(bmap->div[i][0]))
4176	continue;
4177	if (!isl_int_is_zero(bmap->div[i][1 + 1 + v_div + div]))
4178	return n_div;
4179	}
4180
4181	isl_int_add(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
4182	if (isl_int_is_neg(bmap->ineq[l][0])) {
4183	isl_int_sub(bmap->ineq[l][0],
4184	bmap->ineq[l][0], bmap->ineq[u][0]);
4185	bmap = isl_basic_map_copy(bmap);
4186	bmap = isl_basic_map_set_to_empty(bmap);
4187	isl_basic_map_free(bmap);
4188	return n_div;
4189	}
4190	isl_int_add_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
4191	coalesce = n_div;
4192	for (i = 0; i < n_div; ++i) {
4193	if (i == div)
4194	continue;
4195	if (!pairs[i])
4196	continue;
4197	for (j = 0; j < n_div; ++j) {
4198	if (isl_int_is_zero(bmap->div[j][0]))
4199	continue;
4200	if (!isl_int_is_zero(bmap->div[j][1 + 1 + v_div + i]))
4201	break;
4202	}
4203	if (j < n_div)
4204	continue;
4205	for (j = 0; j < bmap->n_ineq; ++j) {
4206	int valid;
4207	if (j == l || j == u)
4208	continue;
4209	if (isl_int_is_zero(bmap->ineq[j][1 + v_div + div])) {
4210	if (is_zero_or_one(v: bmap->ineq[j][1 + v_div + i]))
4211	continue;
4212	break;
4213	}
4214	if (isl_int_is_zero(bmap->ineq[j][1 + v_div + i]))
4215	break;
4216	isl_int_mul(bmap->ineq[j][1 + v_div + div],
4217	bmap->ineq[j][1 + v_div + div],
4218	bmap->ineq[l][0]);
4219	valid = isl_int_eq(bmap->ineq[j][1 + v_div + div],
4220	bmap->ineq[j][1 + v_div + i]);
4221	isl_int_divexact(bmap->ineq[j][1 + v_div + div],
4222	bmap->ineq[j][1 + v_div + div],
4223	bmap->ineq[l][0]);
4224	if (!valid)
4225	break;
4226	}
4227	if (j < bmap->n_ineq)
4228	continue;
4229	coalesce = i;
4230	break;
4231	}
4232	isl_int_sub_ui(bmap->ineq[l][0], bmap->ineq[l][0], 1);
4233	isl_int_sub(bmap->ineq[l][0], bmap->ineq[l][0], bmap->ineq[u][0]);
4234	return coalesce;
4235	}
4236
4237	/* Internal data structure used during the construction and/or evaluation of
4238	* an inequality that ensures that a pair of bounds always allows
4239	* for an integer value.
4240	*
4241	* "tab" is the tableau in which the inequality is evaluated. It may
4242	* be NULL until it is actually needed.
4243	* "v" contains the inequality coefficients.
4244	* "g", "fl" and "fu" are temporary scalars used during the construction and
4245	* evaluation.
4246	*/
4247	struct test_ineq_data {
4248	struct isl_tab *tab;
4249	isl_vec *v;
4250	isl_int g;
4251	isl_int fl;
4252	isl_int fu;
4253	};
4254
4255	/* Free all the memory allocated by the fields of "data".
4256	*/
4257	static void test_ineq_data_clear(struct test_ineq_data *data)
4258	{
4259	isl_tab_free(tab: data->tab);
4260	isl_vec_free(vec: data->v);
4261	isl_int_clear(data->g);
4262	isl_int_clear(data->fl);
4263	isl_int_clear(data->fu);
4264	}
4265
4266	/* Is the inequality stored in data->v satisfied by "bmap"?
4267	* That is, does it only attain non-negative values?
4268	* data->tab is a tableau corresponding to "bmap".
4269	*/
4270	static isl_bool test_ineq_is_satisfied(__isl_keep isl_basic_map *bmap,
4271	struct test_ineq_data *data)
4272	{
4273	isl_ctx *ctx;
4274	enum isl_lp_result res;
4275
4276	ctx = isl_basic_map_get_ctx(bmap);
4277	if (!data->tab)
4278	data->tab = isl_tab_from_basic_map(bmap, track: 0);
4279	res = isl_tab_min(tab: data->tab, f: data->v->el, denom: ctx->one, opt: &data->g, NULL, flags: 0);
4280	if (res == isl_lp_error)
4281	return isl_bool_error;
4282	return res == isl_lp_ok && isl_int_is_nonneg(data->g);
4283	}
4284
4285	/* Given a lower and an upper bound on div i, do they always allow
4286	* for an integer value of the given div?
4287	* Determine this property by constructing an inequality
4288	* such that the property is guaranteed when the inequality is nonnegative.
4289	* The lower bound is inequality l, while the upper bound is inequality u.
4290	* The constructed inequality is stored in data->v.
4291	*
4292	* Let the upper bound be
4293	*
4294	* -n_u a + e_u >= 0
4295	*
4296	* and the lower bound
4297	*
4298	* n_l a + e_l >= 0
4299	*
4300	* Let n_u = f_u g and n_l = f_l g, with g = gcd(n_u, n_l).
4301	* We have
4302	*
4303	* - f_u e_l <= f_u f_l g a <= f_l e_u
4304	*
4305	* Since all variables are integer valued, this is equivalent to
4306	*
4307	* - f_u e_l - (f_u - 1) <= f_u f_l g a <= f_l e_u + (f_l - 1)
4308	*
4309	* If this interval is at least f_u f_l g, then it contains at least
4310	* one integer value for a.
4311	* That is, the test constraint is
4312	*
4313	* f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 >= f_u f_l g
4314	*
4315	* or
4316	*
4317	* f_l e_u + f_u e_l + f_l - 1 + f_u - 1 + 1 - f_u f_l g >= 0
4318	*
4319	* If the coefficients of f_l e_u + f_u e_l have a common divisor g',
4320	* then the constraint can be scaled down by a factor g',
4321	* with the constant term replaced by
4322	* floor((f_l e_{u,0} + f_u e_{l,0} + f_l - 1 + f_u - 1 + 1 - f_u f_l g)/g').
4323	* Note that the result of applying Fourier-Motzkin to this pair
4324	* of constraints is
4325	*
4326	* f_l e_u + f_u e_l >= 0
4327	*
4328	* If the constant term of the scaled down version of this constraint,
4329	* i.e., floor((f_l e_{u,0} + f_u e_{l,0})/g') is equal to the constant
4330	* term of the scaled down test constraint, then the test constraint
4331	* is known to hold and no explicit evaluation is required.
4332	* This is essentially the Omega test.
4333	*
4334	* If the test constraint consists of only a constant term, then
4335	* it is sufficient to look at the sign of this constant term.
4336	*/
4337	static isl_bool int_between_bounds(__isl_keep isl_basic_map *bmap, int i,
4338	int l, int u, struct test_ineq_data *data)
4339	{
4340	unsigned offset;
4341	isl_size n_div;
4342
4343	offset = isl_basic_map_offset(bmap, type: isl_dim_div);
4344	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4345	if (n_div < 0)
4346	return isl_bool_error;
4347
4348	isl_int_gcd(data->g,
4349	bmap->ineq[l][offset + i], bmap->ineq[u][offset + i]);
4350	isl_int_divexact(data->fl, bmap->ineq[l][offset + i], data->g);
4351	isl_int_divexact(data->fu, bmap->ineq[u][offset + i], data->g);
4352	isl_int_neg(data->fu, data->fu);
4353	isl_seq_combine(dst: data->v->el, m1: data->fl, src1: bmap->ineq[u],
4354	m2: data->fu, src2: bmap->ineq[l], len: offset + n_div);
4355	isl_int_mul(data->g, data->g, data->fl);
4356	isl_int_mul(data->g, data->g, data->fu);
4357	isl_int_sub(data->g, data->g, data->fl);
4358	isl_int_sub(data->g, data->g, data->fu);
4359	isl_int_add_ui(data->g, data->g, 1);
4360	isl_int_sub(data->fl, data->v->el[0], data->g);
4361
4362	isl_seq_gcd(p: data->v->el + 1, len: offset - 1 + n_div, gcd: &data->g);
4363	if (isl_int_is_zero(data->g))
4364	return isl_int_is_nonneg(data->fl);
4365	if (isl_int_is_one(data->g)) {
4366	isl_int_set(data->v->el[0], data->fl);
4367	return test_ineq_is_satisfied(bmap, data);
4368	}
4369	isl_int_fdiv_q(data->fl, data->fl, data->g);
4370	isl_int_fdiv_q(data->v->el[0], data->v->el[0], data->g);
4371	if (isl_int_eq(data->fl, data->v->el[0]))
4372	return isl_bool_true;
4373	isl_int_set(data->v->el[0], data->fl);
4374	isl_seq_scale_down(dst: data->v->el + 1, src: data->v->el + 1, f: data->g,
4375	len: offset - 1 + n_div);
4376
4377	return test_ineq_is_satisfied(bmap, data);
4378	}
4379
4380	/* Remove more kinds of divs that are not strictly needed.
4381	* In particular, if all pairs of lower and upper bounds on a div
4382	* are such that they allow at least one integer value of the div,
4383	* then we can eliminate the div using Fourier-Motzkin without
4384	* introducing any spurious solutions.
4385	*
4386	* If at least one of the two constraints has a unit coefficient for the div,
4387	* then the presence of such a value is guaranteed so there is no need to check.
4388	* In particular, the value attained by the bound with unit coefficient
4389	* can serve as this intermediate value.
4390	*/
4391	static __isl_give isl_basic_map *drop_more_redundant_divs(
4392	__isl_take isl_basic_map *bmap, __isl_take int *pairs, int n)
4393	{
4394	isl_ctx *ctx;
4395	struct test_ineq_data data = { NULL, NULL };
4396	unsigned off;
4397	isl_size n_div;
4398	int remove = -1;
4399
4400	isl_int_init(data.g);
4401	isl_int_init(data.fl);
4402	isl_int_init(data.fu);
4403
4404	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4405	if (n_div < 0)
4406	goto error;
4407
4408	ctx = isl_basic_map_get_ctx(bmap);
4409	off = isl_basic_map_offset(bmap, type: isl_dim_div);
4410	data.v = isl_vec_alloc(ctx, size: off + n_div);
4411	if (!data.v)
4412	goto error;
4413
4414	while (n > 0) {
4415	int i, l, u;
4416	int best = -1;
4417	isl_bool has_int;
4418
4419	for (i = 0; i < n_div; ++i) {
4420	if (!pairs[i])
4421	continue;
4422	if (best >= 0 && pairs[best] <= pairs[i])
4423	continue;
4424	best = i;
4425	}
4426
4427	i = best;
4428	for (l = 0; l < bmap->n_ineq; ++l) {
4429	if (!isl_int_is_pos(bmap->ineq[l][off + i]))
4430	continue;
4431	if (isl_int_is_one(bmap->ineq[l][off + i]))
4432	continue;
4433	for (u = 0; u < bmap->n_ineq; ++u) {
4434	if (!isl_int_is_neg(bmap->ineq[u][off + i]))
4435	continue;
4436	if (isl_int_is_negone(bmap->ineq[u][off + i]))
4437	continue;
4438	has_int = int_between_bounds(bmap, i, l, u,
4439	data: &data);
4440	if (has_int < 0)
4441	goto error;
4442	if (data.tab && data.tab->empty)
4443	break;
4444	if (!has_int)
4445	break;
4446	}
4447	if (u < bmap->n_ineq)
4448	break;
4449	}
4450	if (data.tab && data.tab->empty) {
4451	bmap = isl_basic_map_set_to_empty(bmap);
4452	break;
4453	}
4454	if (l == bmap->n_ineq) {
4455	remove = i;
4456	break;
4457	}
4458	pairs[i] = 0;
4459	--n;
4460	}
4461
4462	test_ineq_data_clear(data: &data);
4463
4464	free(ptr: pairs);
4465
4466	if (remove < 0)
4467	return bmap;
4468
4469	bmap = isl_basic_map_remove_dims(bmap, type: isl_dim_div, first: remove, n: 1);
4470	return isl_basic_map_drop_redundant_divs(bmap);
4471	error:
4472	free(ptr: pairs);
4473	isl_basic_map_free(bmap);
4474	test_ineq_data_clear(data: &data);
4475	return NULL;
4476	}
4477
4478	/* Given a pair of divs div1 and div2 such that, except for the lower bound l
4479	* and the upper bound u, div1 always occurs together with div2 in the form
4480	* (div1 + m div2), where m is the constant range on the variable div1
4481	* allowed by l and u, replace the pair div1 and div2 by a single
4482	* div that is equal to div1 + m div2.
4483	*
4484	* The new div will appear in the location that contains div2.
4485	* We need to modify all constraints that contain
4486	* div2 = (div - div1) / m
4487	* The coefficient of div2 is known to be equal to 1 or -1.
4488	* (If a constraint does not contain div2, it will also not contain div1.)
4489	* If the constraint also contains div1, then we know they appear
4490	* as f (div1 + m div2) and we can simply replace (div1 + m div2) by div,
4491	* i.e., the coefficient of div is f.
4492	*
4493	* Otherwise, we first need to introduce div1 into the constraint.
4494	* Let l be
4495	*
4496	* div1 + f >=0
4497	*
4498	* and u
4499	*
4500	* -div1 + f' >= 0
4501	*
4502	* A lower bound on div2
4503	*
4504	* div2 + t >= 0
4505	*
4506	* can be replaced by
4507	*
4508	* m div2 + div1 + m t + f >= 0
4509	*
4510	* An upper bound
4511	*
4512	* -div2 + t >= 0
4513	*
4514	* can be replaced by
4515	*
4516	* -(m div2 + div1) + m t + f' >= 0
4517	*
4518	* These constraint are those that we would obtain from eliminating
4519	* div1 using Fourier-Motzkin.
4520	*
4521	* After all constraints have been modified, we drop the lower and upper
4522	* bound and then drop div1.
4523	* Since the new div is only placed in the same location that used
4524	* to store div2, but otherwise has a different meaning, any possible
4525	* explicit representation of the original div2 is removed.
4526	*/
4527	static __isl_give isl_basic_map *coalesce_divs(__isl_take isl_basic_map *bmap,
4528	unsigned div1, unsigned div2, unsigned l, unsigned u)
4529	{
4530	isl_ctx *ctx;
4531	isl_int m;
4532	isl_size v_div;
4533	unsigned total;
4534	int i;
4535
4536	ctx = isl_basic_map_get_ctx(bmap);
4537
4538	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
4539	if (v_div < 0)
4540	return isl_basic_map_free(bmap);
4541	total = 1 + v_div + bmap->n_div;
4542
4543	isl_int_init(m);
4544	isl_int_add(m, bmap->ineq[l][0], bmap->ineq[u][0]);
4545	isl_int_add_ui(m, m, 1);
4546
4547	for (i = 0; i < bmap->n_ineq; ++i) {
4548	if (i == l || i == u)
4549	continue;
4550	if (isl_int_is_zero(bmap->ineq[i][1 + v_div + div2]))
4551	continue;
4552	if (isl_int_is_zero(bmap->ineq[i][1 + v_div + div1])) {
4553	if (isl_int_is_pos(bmap->ineq[i][1 + v_div + div2]))
4554	isl_seq_combine(dst: bmap->ineq[i], m1: m, src1: bmap->ineq[i],
4555	m2: ctx->one, src2: bmap->ineq[l], len: total);
4556	else
4557	isl_seq_combine(dst: bmap->ineq[i], m1: m, src1: bmap->ineq[i],
4558	m2: ctx->one, src2: bmap->ineq[u], len: total);
4559	}
4560	isl_int_set(bmap->ineq[i][1 + v_div + div2],
4561	bmap->ineq[i][1 + v_div + div1]);
4562	isl_int_set_si(bmap->ineq[i][1 + v_div + div1], 0);
4563	}
4564
4565	isl_int_clear(m);
4566	if (l > u) {
4567	isl_basic_map_drop_inequality(bmap, pos: l);
4568	isl_basic_map_drop_inequality(bmap, pos: u);
4569	} else {
4570	isl_basic_map_drop_inequality(bmap, pos: u);
4571	isl_basic_map_drop_inequality(bmap, pos: l);
4572	}
4573	bmap = isl_basic_map_mark_div_unknown(bmap, div: div2);
4574	bmap = isl_basic_map_drop_div(bmap, div: div1);
4575	return bmap;
4576	}
4577
4578	/* First check if we can coalesce any pair of divs and
4579	* then continue with dropping more redundant divs.
4580	*
4581	* We loop over all pairs of lower and upper bounds on a div
4582	* with coefficient 1 and -1, respectively, check if there
4583	* is any other div "c" with which we can coalesce the div
4584	* and if so, perform the coalescing.
4585	*/
4586	static __isl_give isl_basic_map *coalesce_or_drop_more_redundant_divs(
4587	__isl_take isl_basic_map *bmap, int *pairs, int n)
4588	{
4589	int i, l, u;
4590	isl_size v_div;
4591	isl_size n_div;
4592
4593	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
4594	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4595	if (v_div < 0 || n_div < 0)
4596	return isl_basic_map_free(bmap);
4597
4598	for (i = 0; i < n_div; ++i) {
4599	if (!pairs[i])
4600	continue;
4601	for (l = 0; l < bmap->n_ineq; ++l) {
4602	if (!isl_int_is_one(bmap->ineq[l][1 + v_div + i]))
4603	continue;
4604	for (u = 0; u < bmap->n_ineq; ++u) {
4605	int c;
4606
4607	if (!isl_int_is_negone(bmap->ineq[u][1+v_div+i]))
4608	continue;
4609	c = div_find_coalesce(bmap, pairs, div: i, l, u);
4610	if (c < 0)
4611	goto error;
4612	if (c >= n_div)
4613	continue;
4614	free(ptr: pairs);
4615	bmap = coalesce_divs(bmap, div1: i, div2: c, l, u);
4616	return isl_basic_map_drop_redundant_divs(bmap);
4617	}
4618	}
4619	}
4620
4621	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_EMPTY)) {
4622	free(ptr: pairs);
4623	return bmap;
4624	}
4625
4626	return drop_more_redundant_divs(bmap, pairs, n);
4627	error:
4628	free(ptr: pairs);
4629	isl_basic_map_free(bmap);
4630	return NULL;
4631	}
4632
4633	/* Are the "n" coefficients starting at "first" of inequality constraints
4634	* "i" and "j" of "bmap" equal to each other?
4635	*/
4636	static int is_parallel_part(__isl_keep isl_basic_map *bmap, int i, int j,
4637	int first, int n)
4638	{
4639	return isl_seq_eq(p1: bmap->ineq[i] + first, p2: bmap->ineq[j] + first, len: n);
4640	}
4641
4642	/* Are inequality constraints "i" and "j" of "bmap" equal to each other,
4643	* apart from the constant term and the coefficient at position "pos"?
4644	*/
4645	static isl_bool is_parallel_except(__isl_keep isl_basic_map *bmap, int i, int j,
4646	int pos)
4647	{
4648	isl_size total;
4649
4650	total = isl_basic_map_dim(bmap, type: isl_dim_all);
4651	if (total < 0)
4652	return isl_bool_error;
4653	return is_parallel_part(bmap, i, j, first: 1, n: pos - 1) &&
4654	is_parallel_part(bmap, i, j, first: pos + 1, n: total - pos);
4655	}
4656
4657	/* Are inequality constraints "i" and "j" of "bmap" opposite to each other,
4658	* apart from the constant term and the coefficient at position "pos"?
4659	*/
4660	static isl_bool is_opposite_except(__isl_keep isl_basic_map *bmap, int i, int j,
4661	int pos)
4662	{
4663	isl_size total;
4664
4665	total = isl_basic_map_dim(bmap, type: isl_dim_all);
4666	if (total < 0)
4667	return isl_bool_error;
4668	return is_opposite_part(bmap, i, j, first: 1, n: pos - 1) &&
4669	is_opposite_part(bmap, i, j, first: pos + 1, n: total - pos);
4670	}
4671
4672	/* Restart isl_basic_map_drop_redundant_divs after "bmap" has
4673	* been modified, simplying it if "simplify" is set.
4674	* Free the temporary data structure "pairs" that was associated
4675	* to the old version of "bmap".
4676	*/
4677	static __isl_give isl_basic_map *drop_redundant_divs_again(
4678	__isl_take isl_basic_map *bmap, __isl_take int *pairs, int simplify)
4679	{
4680	if (simplify)
4681	bmap = isl_basic_map_simplify(bmap);
4682	free(ptr: pairs);
4683	return isl_basic_map_drop_redundant_divs(bmap);
4684	}
4685
4686	/* Is "div" the single unknown existentially quantified variable
4687	* in inequality constraint "ineq" of "bmap"?
4688	* "div" is known to have a non-zero coefficient in "ineq".
4689	*/
4690	static isl_bool single_unknown(__isl_keep isl_basic_map *bmap, int ineq,
4691	int div)
4692	{
4693	int i;
4694	isl_size n_div;
4695	unsigned o_div;
4696	isl_bool known;
4697
4698	known = isl_basic_map_div_is_known(bmap, div);
4699	if (known < 0 || known)
4700	return isl_bool_not(b: known);
4701	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4702	if (n_div < 0)
4703	return isl_bool_error;
4704	if (n_div == 1)
4705	return isl_bool_true;
4706	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
4707	for (i = 0; i < n_div; ++i) {
4708	isl_bool known;
4709
4710	if (i == div)
4711	continue;
4712	if (isl_int_is_zero(bmap->ineq[ineq][o_div + i]))
4713	continue;
4714	known = isl_basic_map_div_is_known(bmap, div: i);
4715	if (known < 0 || !known)
4716	return known;
4717	}
4718
4719	return isl_bool_true;
4720	}
4721
4722	/* Does integer division "div" have coefficient 1 in inequality constraint
4723	* "ineq" of "map"?
4724	*/
4725	static isl_bool has_coef_one(__isl_keep isl_basic_map *bmap, int div, int ineq)
4726	{
4727	unsigned o_div;
4728
4729	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
4730	if (isl_int_is_one(bmap->ineq[ineq][o_div + div]))
4731	return isl_bool_true;
4732
4733	return isl_bool_false;
4734	}
4735
4736	/* Turn inequality constraint "ineq" of "bmap" into an equality and
4737	* then try and drop redundant divs again,
4738	* freeing the temporary data structure "pairs" that was associated
4739	* to the old version of "bmap".
4740	*/
4741	static __isl_give isl_basic_map *set_eq_and_try_again(
4742	__isl_take isl_basic_map *bmap, int ineq, __isl_take int *pairs)
4743	{
4744	bmap = isl_basic_map_cow(bmap);
4745	isl_basic_map_inequality_to_equality(bmap, pos: ineq);
4746	return drop_redundant_divs_again(bmap, pairs, simplify: 1);
4747	}
4748
4749	/* Drop the integer division at position "div", along with the two
4750	* inequality constraints "ineq1" and "ineq2" in which it appears
4751	* from "bmap" and then try and drop redundant divs again,
4752	* freeing the temporary data structure "pairs" that was associated
4753	* to the old version of "bmap".
4754	*/
4755	static __isl_give isl_basic_map *drop_div_and_try_again(
4756	__isl_take isl_basic_map *bmap, int div, int ineq1, int ineq2,
4757	__isl_take int *pairs)
4758	{
4759	if (ineq1 > ineq2) {
4760	isl_basic_map_drop_inequality(bmap, pos: ineq1);
4761	isl_basic_map_drop_inequality(bmap, pos: ineq2);
4762	} else {
4763	isl_basic_map_drop_inequality(bmap, pos: ineq2);
4764	isl_basic_map_drop_inequality(bmap, pos: ineq1);
4765	}
4766	bmap = isl_basic_map_drop_div(bmap, div);
4767	return drop_redundant_divs_again(bmap, pairs, simplify: 0);
4768	}
4769
4770	/* Given two inequality constraints
4771	*
4772	* f(x) + n d + c >= 0, (ineq)
4773	*
4774	* with d the variable at position "pos", and
4775	*
4776	* f(x) + c0 >= 0, (lower)
4777	*
4778	* compute the maximal value of the lower bound ceil((-f(x) - c)/n)
4779	* determined by the first constraint.
4780	* That is, store
4781	*
4782	* ceil((c0 - c)/n)
4783	*
4784	* in *l.
4785	*/
4786	static void lower_bound_from_parallel(__isl_keep isl_basic_map *bmap,
4787	int ineq, int lower, int pos, isl_int *l)
4788	{
4789	isl_int_neg(*l, bmap->ineq[ineq][0]);
4790	isl_int_add(*l, *l, bmap->ineq[lower][0]);
4791	isl_int_cdiv_q(*l, *l, bmap->ineq[ineq][pos]);
4792	}
4793
4794	/* Given two inequality constraints
4795	*
4796	* f(x) + n d + c >= 0, (ineq)
4797	*
4798	* with d the variable at position "pos", and
4799	*
4800	* -f(x) - c0 >= 0, (upper)
4801	*
4802	* compute the minimal value of the lower bound ceil((-f(x) - c)/n)
4803	* determined by the first constraint.
4804	* That is, store
4805	*
4806	* ceil((-c1 - c)/n)
4807	*
4808	* in *u.
4809	*/
4810	static void lower_bound_from_opposite(__isl_keep isl_basic_map *bmap,
4811	int ineq, int upper, int pos, isl_int *u)
4812	{
4813	isl_int_neg(*u, bmap->ineq[ineq][0]);
4814	isl_int_sub(*u, *u, bmap->ineq[upper][0]);
4815	isl_int_cdiv_q(*u, *u, bmap->ineq[ineq][pos]);
4816	}
4817
4818	/* Given a lower bound constraint "ineq" on "div" in "bmap",
4819	* does the corresponding lower bound have a fixed value in "bmap"?
4820	*
4821	* In particular, "ineq" is of the form
4822	*
4823	* f(x) + n d + c >= 0
4824	*
4825	* with n > 0, c the constant term and
4826	* d the existentially quantified variable "div".
4827	* That is, the lower bound is
4828	*
4829	* ceil((-f(x) - c)/n)
4830	*
4831	* Look for a pair of constraints
4832	*
4833	* f(x) + c0 >= 0
4834	* -f(x) + c1 >= 0
4835	*
4836	* i.e., -c1 <= -f(x) <= c0, that fix ceil((-f(x) - c)/n) to a constant value.
4837	* That is, check that
4838	*
4839	* ceil((-c1 - c)/n) = ceil((c0 - c)/n)
4840	*
4841	* If so, return the index of inequality f(x) + c0 >= 0.
4842	* Otherwise, return bmap->n_ineq.
4843	* Return -1 on error.
4844	*/
4845	static int lower_bound_is_cst(__isl_keep isl_basic_map *bmap, int div, int ineq)
4846	{
4847	int i;
4848	int lower = -1, upper = -1;
4849	unsigned o_div;
4850	isl_int l, u;
4851	int equal;
4852
4853	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
4854	for (i = 0; i < bmap->n_ineq && (lower < 0 || upper < 0); ++i) {
4855	isl_bool par, opp;
4856
4857	if (i == ineq)
4858	continue;
4859	if (!isl_int_is_zero(bmap->ineq[i][o_div + div]))
4860	continue;
4861	par = isl_bool_false;
4862	if (lower < 0)
4863	par = is_parallel_except(bmap, i: ineq, j: i, pos: o_div + div);
4864	if (par < 0)
4865	return -1;
4866	if (par) {
4867	lower = i;
4868	continue;
4869	}
4870	opp = isl_bool_false;
4871	if (upper < 0)
4872	opp = is_opposite_except(bmap, i: ineq, j: i, pos: o_div + div);
4873	if (opp < 0)
4874	return -1;
4875	if (opp)
4876	upper = i;
4877	}
4878
4879	if (lower < 0 || upper < 0)
4880	return bmap->n_ineq;
4881
4882	isl_int_init(l);
4883	isl_int_init(u);
4884
4885	lower_bound_from_parallel(bmap, ineq, lower, pos: o_div + div, l: &l);
4886	lower_bound_from_opposite(bmap, ineq, upper, pos: o_div + div, u: &u);
4887
4888	equal = isl_int_eq(l, u);
4889
4890	isl_int_clear(l);
4891	isl_int_clear(u);
4892
4893	return equal ? lower : bmap->n_ineq;
4894	}
4895
4896	/* Given a lower bound constraint "ineq" on the existentially quantified
4897	* variable "div", such that the corresponding lower bound has
4898	* a fixed value in "bmap", assign this fixed value to the variable and
4899	* then try and drop redundant divs again,
4900	* freeing the temporary data structure "pairs" that was associated
4901	* to the old version of "bmap".
4902	* "lower" determines the constant value for the lower bound.
4903	*
4904	* In particular, "ineq" is of the form
4905	*
4906	* f(x) + n d + c >= 0,
4907	*
4908	* while "lower" is of the form
4909	*
4910	* f(x) + c0 >= 0
4911	*
4912	* The lower bound is ceil((-f(x) - c)/n) and its constant value
4913	* is ceil((c0 - c)/n).
4914	*/
4915	static __isl_give isl_basic_map *fix_cst_lower(__isl_take isl_basic_map *bmap,
4916	int div, int ineq, int lower, int *pairs)
4917	{
4918	isl_int c;
4919	unsigned o_div;
4920
4921	isl_int_init(c);
4922
4923	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
4924	lower_bound_from_parallel(bmap, ineq, lower, pos: o_div + div, l: &c);
4925	bmap = isl_basic_map_fix(bmap, type: isl_dim_div, pos: div, value: c);
4926	free(ptr: pairs);
4927
4928	isl_int_clear(c);
4929
4930	return isl_basic_map_drop_redundant_divs(bmap);
4931	}
4932
4933	/* Do any of the integer divisions of "bmap" involve integer division "div"?
4934	*
4935	* The integer division "div" could only ever appear in any later
4936	* integer division (with an explicit representation).
4937	*/
4938	static isl_bool any_div_involves_div(__isl_keep isl_basic_map *bmap, int div)
4939	{
4940	int i;
4941	isl_size v_div, n_div;
4942
4943	v_div = isl_basic_map_var_offset(bmap, type: isl_dim_div);
4944	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
4945	if (v_div < 0 || n_div < 0)
4946	return isl_bool_error;
4947
4948	for (i = div + 1; i < n_div; ++i) {
4949	isl_bool unknown;
4950
4951	unknown = isl_basic_map_div_is_marked_unknown(bmap, div: i);
4952	if (unknown < 0)
4953	return isl_bool_error;
4954	if (unknown)
4955	continue;
4956	if (!isl_int_is_zero(bmap->div[i][1 + 1 + v_div + div]))
4957	return isl_bool_true;
4958	}
4959
4960	return isl_bool_false;
4961	}
4962
4963	/* Remove divs that are not strictly needed based on the inequality
4964	* constraints.
4965	* In particular, if a div only occurs positively (or negatively)
4966	* in constraints, then it can simply be dropped.
4967	* Also, if a div occurs in only two constraints and if moreover
4968	* those two constraints are opposite to each other, except for the constant
4969	* term and if the sum of the constant terms is such that for any value
4970	* of the other values, there is always at least one integer value of the
4971	* div, i.e., if one plus this sum is greater than or equal to
4972	* the (absolute value) of the coefficient of the div in the constraints,
4973	* then we can also simply drop the div.
4974	*
4975	* If an existentially quantified variable does not have an explicit
4976	* representation, appears in only a single lower bound that does not
4977	* involve any other such existentially quantified variables and appears
4978	* in this lower bound with coefficient 1,
4979	* then fix the variable to the value of the lower bound. That is,
4980	* turn the inequality into an equality.
4981	* If for any value of the other variables, there is any value
4982	* for the existentially quantified variable satisfying the constraints,
4983	* then this lower bound also satisfies the constraints.
4984	* It is therefore safe to pick this lower bound.
4985	*
4986	* The same reasoning holds even if the coefficient is not one.
4987	* However, fixing the variable to the value of the lower bound may
4988	* in general introduce an extra integer division, in which case
4989	* it may be better to pick another value.
4990	* If this integer division has a known constant value, then plugging
4991	* in this constant value removes the existentially quantified variable
4992	* completely. In particular, if the lower bound is of the form
4993	* ceil((-f(x) - c)/n) and there are two constraints, f(x) + c0 >= 0 and
4994	* -f(x) + c1 >= 0 such that ceil((-c1 - c)/n) = ceil((c0 - c)/n),
4995	* then the existentially quantified variable can be assigned this
4996	* shared value.
4997	*
4998	* We skip divs that appear in equalities or in the definition of other divs.
4999	* Divs that appear in the definition of other divs usually occur in at least
5000	* 4 constraints, but the constraints may have been simplified.
5001	*
5002	* If any divs are left after these simple checks then we move on
5003	* to more complicated cases in drop_more_redundant_divs.
5004	*/
5005	static __isl_give isl_basic_map *isl_basic_map_drop_redundant_divs_ineq(
5006	__isl_take isl_basic_map *bmap)
5007	{
5008	int i, j;
5009	isl_size off;
5010	int *pairs = NULL;
5011	int n = 0;
5012	isl_size n_ineq;
5013
5014	if (!bmap)
5015	goto error;
5016	if (bmap->n_div == 0)
5017	return bmap;
5018
5019	off = isl_basic_map_var_offset(bmap, type: isl_dim_div);
5020	if (off < 0)
5021	return isl_basic_map_free(bmap);
5022	pairs = isl_calloc_array(bmap->ctx, int, bmap->n_div);
5023	if (!pairs)
5024	goto error;
5025
5026	n_ineq = isl_basic_map_n_inequality(bmap);
5027	if (n_ineq < 0)
5028	goto error;
5029	for (i = 0; i < bmap->n_div; ++i) {
5030	int pos, neg;
5031	int last_pos, last_neg;
5032	int redundant;
5033	int defined;
5034	isl_bool involves, opp, set_div;
5035
5036	defined = !isl_int_is_zero(bmap->div[i][0]);
5037	involves = any_div_involves_div(bmap, div: i);
5038	if (involves < 0)
5039	goto error;
5040	if (involves)
5041	continue;
5042	for (j = 0; j < bmap->n_eq; ++j)
5043	if (!isl_int_is_zero(bmap->eq[j][1 + off + i]))
5044	break;
5045	if (j < bmap->n_eq)
5046	continue;
5047	++n;
5048	pos = neg = 0;
5049	for (j = 0; j < bmap->n_ineq; ++j) {
5050	if (isl_int_is_pos(bmap->ineq[j][1 + off + i])) {
5051	last_pos = j;
5052	++pos;
5053	}
5054	if (isl_int_is_neg(bmap->ineq[j][1 + off + i])) {
5055	last_neg = j;
5056	++neg;
5057	}
5058	}
5059	pairs[i] = pos * neg;
5060	if (pairs[i] == 0) {
5061	for (j = bmap->n_ineq - 1; j >= 0; --j)
5062	if (!isl_int_is_zero(bmap->ineq[j][1+off+i]))
5063	isl_basic_map_drop_inequality(bmap, pos: j);
5064	bmap = isl_basic_map_drop_div(bmap, div: i);
5065	return drop_redundant_divs_again(bmap, pairs, simplify: 0);
5066	}
5067	if (pairs[i] != 1)
5068	opp = isl_bool_false;
5069	else
5070	opp = is_opposite(bmap, i: last_pos, j: last_neg);
5071	if (opp < 0)
5072	goto error;
5073	if (!opp) {
5074	int lower;
5075	isl_bool single, one;
5076
5077	if (pos != 1)
5078	continue;
5079	single = single_unknown(bmap, ineq: last_pos, div: i);
5080	if (single < 0)
5081	goto error;
5082	if (!single)
5083	continue;
5084	one = has_coef_one(bmap, div: i, ineq: last_pos);
5085	if (one < 0)
5086	goto error;
5087	if (one)
5088	return set_eq_and_try_again(bmap, ineq: last_pos,
5089	pairs);
5090	lower = lower_bound_is_cst(bmap, div: i, ineq: last_pos);
5091	if (lower < 0)
5092	goto error;
5093	if (lower < n_ineq)
5094	return fix_cst_lower(bmap, div: i, ineq: last_pos, lower,
5095	pairs);
5096	continue;
5097	}
5098
5099	isl_int_add(bmap->ineq[last_pos][0],
5100	bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
5101	isl_int_add_ui(bmap->ineq[last_pos][0],
5102	bmap->ineq[last_pos][0], 1);
5103	redundant = isl_int_ge(bmap->ineq[last_pos][0],
5104	bmap->ineq[last_pos][1+off+i]);
5105	isl_int_sub_ui(bmap->ineq[last_pos][0],
5106	bmap->ineq[last_pos][0], 1);
5107	isl_int_sub(bmap->ineq[last_pos][0],
5108	bmap->ineq[last_pos][0], bmap->ineq[last_neg][0]);
5109	if (redundant)
5110	return drop_div_and_try_again(bmap, div: i,
5111	ineq1: last_pos, ineq2: last_neg, pairs);
5112	if (defined)
5113	set_div = isl_bool_false;
5114	else
5115	set_div = ok_to_set_div_from_bound(bmap, div: i, ineq: last_pos);
5116	if (set_div < 0)
5117	return isl_basic_map_free(bmap);
5118	if (set_div) {
5119	bmap = set_div_from_lower_bound(bmap, div: i, ineq: last_pos);
5120	return drop_redundant_divs_again(bmap, pairs, simplify: 1);
5121	}
5122	pairs[i] = 0;
5123	--n;
5124	}
5125
5126	if (n > 0)
5127	return coalesce_or_drop_more_redundant_divs(bmap, pairs, n);
5128
5129	free(ptr: pairs);
5130	return bmap;
5131	error:
5132	free(ptr: pairs);
5133	isl_basic_map_free(bmap);
5134	return NULL;
5135	}
5136
5137	/* Consider the coefficients at "c" as a row vector and replace
5138	* them with their product with "T". "T" is assumed to be a square matrix.
5139	*/
5140	static isl_stat preimage(isl_int *c, __isl_keep isl_mat *T)
5141	{
5142	isl_size n;
5143	isl_ctx *ctx;
5144	isl_vec *v;
5145
5146	n = isl_mat_rows(mat: T);
5147	if (n < 0)
5148	return isl_stat_error;
5149	if (isl_seq_first_non_zero(p: c, len: n) == -1)
5150	return isl_stat_ok;
5151	ctx = isl_mat_get_ctx(mat: T);
5152	v = isl_vec_alloc(ctx, size: n);
5153	if (!v)
5154	return isl_stat_error;
5155	isl_seq_swp_or_cpy(dst: v->el, src: c, len: n);
5156	v = isl_vec_mat_product(vec: v, mat: isl_mat_copy(mat: T));
5157	if (!v)
5158	return isl_stat_error;
5159	isl_seq_swp_or_cpy(dst: c, src: v->el, len: n);
5160	isl_vec_free(vec: v);
5161
5162	return isl_stat_ok;
5163	}
5164
5165	/* Plug in T for the variables in "bmap" starting at "pos".
5166	* T is a linear unimodular matrix, i.e., without constant term.
5167	*/
5168	static __isl_give isl_basic_map *isl_basic_map_preimage_vars(
5169	__isl_take isl_basic_map *bmap, unsigned pos, __isl_take isl_mat *T)
5170	{
5171	int i;
5172	isl_size n_row, n_col;
5173
5174	bmap = isl_basic_map_cow(bmap);
5175	n_row = isl_mat_rows(mat: T);
5176	n_col = isl_mat_cols(mat: T);
5177	if (!bmap || n_row < 0 || n_col < 0)
5178	goto error;
5179
5180	if (n_col != n_row)
5181	isl_die(isl_mat_get_ctx(T), isl_error_invalid,
5182	"expecting square matrix", goto error);
5183
5184	if (isl_basic_map_check_range(bmap, type: isl_dim_all, first: pos, n: n_col) < 0)
5185	goto error;
5186
5187	for (i = 0; i < bmap->n_eq; ++i)
5188	if (preimage(c: bmap->eq[i] + 1 + pos, T) < 0)
5189	goto error;
5190	for (i = 0; i < bmap->n_ineq; ++i)
5191	if (preimage(c: bmap->ineq[i] + 1 + pos, T) < 0)
5192	goto error;
5193	for (i = 0; i < bmap->n_div; ++i) {
5194	if (isl_basic_map_div_is_marked_unknown(bmap, div: i))
5195	continue;
5196	if (preimage(c: bmap->div[i] + 1 + 1 + pos, T) < 0)
5197	goto error;
5198	}
5199
5200	isl_mat_free(mat: T);
5201	return bmap;
5202	error:
5203	isl_basic_map_free(bmap);
5204	isl_mat_free(mat: T);
5205	return NULL;
5206	}
5207
5208	/* Remove divs that are not strictly needed.
5209	*
5210	* First look for an equality constraint involving two or more
5211	* existentially quantified variables without an explicit
5212	* representation. Replace the combination that appears
5213	* in the equality constraint by a single existentially quantified
5214	* variable such that the equality can be used to derive
5215	* an explicit representation for the variable.
5216	* If there are no more such equality constraints, then continue
5217	* with isl_basic_map_drop_redundant_divs_ineq.
5218	*
5219	* In particular, if the equality constraint is of the form
5220	*
5221	* f(x) + \sum_i c_i a_i = 0
5222	*
5223	* with a_i existentially quantified variable without explicit
5224	* representation, then apply a transformation on the existentially
5225	* quantified variables to turn the constraint into
5226	*
5227	* f(x) + g a_1' = 0
5228	*
5229	* with g the gcd of the c_i.
5230	* In order to easily identify which existentially quantified variables
5231	* have a complete explicit representation, i.e., without being defined
5232	* in terms of other existentially quantified variables without
5233	* an explicit representation, the existentially quantified variables
5234	* are first sorted.
5235	*
5236	* The variable transformation is computed by extending the row
5237	* [c_1/g ... c_n/g] to a unimodular matrix, obtaining the transformation
5238	*
5239	* [a_1'] [c_1/g ... c_n/g] [ a_1 ]
5240	* [a_2'] [ a_2 ]
5241	* ... = U ....
5242	* [a_n'] [ a_n ]
5243	*
5244	* with [c_1/g ... c_n/g] representing the first row of U.
5245	* The inverse of U is then plugged into the original constraints.
5246	* The call to isl_basic_map_simplify makes sure the explicit
5247	* representation for a_1' is extracted from the equality constraint.
5248	*/
5249	__isl_give isl_basic_map *isl_basic_map_drop_redundant_divs(
5250	__isl_take isl_basic_map *bmap)
5251	{
5252	int first;
5253	int i;
5254	unsigned o_div;
5255	isl_size n_div;
5256	int l;
5257	isl_ctx *ctx;
5258	isl_mat *T;
5259
5260	if (!bmap)
5261	return NULL;
5262	if (isl_basic_map_divs_known(bmap))
5263	return isl_basic_map_drop_redundant_divs_ineq(bmap);
5264	if (bmap->n_eq == 0)
5265	return isl_basic_map_drop_redundant_divs_ineq(bmap);
5266	bmap = isl_basic_map_sort_divs(bmap);
5267	if (!bmap)
5268	return NULL;
5269
5270	first = isl_basic_map_first_unknown_div(bmap);
5271	if (first < 0)
5272	return isl_basic_map_free(bmap);
5273
5274	o_div = isl_basic_map_offset(bmap, type: isl_dim_div);
5275	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
5276	if (n_div < 0)
5277	return isl_basic_map_free(bmap);
5278
5279	for (i = 0; i < bmap->n_eq; ++i) {
5280	l = isl_seq_first_non_zero(p: bmap->eq[i] + o_div + first,
5281	len: n_div - (first));
5282	if (l < 0)
5283	continue;
5284	l += first;
5285	if (isl_seq_first_non_zero(p: bmap->eq[i] + o_div + l + 1,
5286	len: n_div - (l + 1)) == -1)
5287	continue;
5288	break;
5289	}
5290	if (i >= bmap->n_eq)
5291	return isl_basic_map_drop_redundant_divs_ineq(bmap);
5292
5293	ctx = isl_basic_map_get_ctx(bmap);
5294	T = isl_mat_alloc(ctx, n_row: n_div - l, n_col: n_div - l);
5295	if (!T)
5296	return isl_basic_map_free(bmap);
5297	isl_seq_cpy(dst: T->row[0], src: bmap->eq[i] + o_div + l, len: n_div - l);
5298	T = isl_mat_normalize_row(mat: T, row: 0);
5299	T = isl_mat_unimodular_complete(M: T, row: 1);
5300	T = isl_mat_right_inverse(mat: T);
5301
5302	for (i = l; i < n_div; ++i)
5303	bmap = isl_basic_map_mark_div_unknown(bmap, div: i);
5304	bmap = isl_basic_map_preimage_vars(bmap, pos: o_div - 1 + l, T);
5305	bmap = isl_basic_map_simplify(bmap);
5306
5307	return isl_basic_map_drop_redundant_divs(bmap);
5308	}
5309
5310	/* Does "bmap" satisfy any equality that involves more than 2 variables
5311	* and/or has coefficients different from -1 and 1?
5312	*/
5313	static isl_bool has_multiple_var_equality(__isl_keep isl_basic_map *bmap)
5314	{
5315	int i;
5316	isl_size total;
5317
5318	total = isl_basic_map_dim(bmap, type: isl_dim_all);
5319	if (total < 0)
5320	return isl_bool_error;
5321
5322	for (i = 0; i < bmap->n_eq; ++i) {
5323	int j, k;
5324
5325	j = isl_seq_first_non_zero(p: bmap->eq[i] + 1, len: total);
5326	if (j < 0)
5327	continue;
5328	if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
5329	!isl_int_is_negone(bmap->eq[i][1 + j]))
5330	return isl_bool_true;
5331
5332	j += 1;
5333	k = isl_seq_first_non_zero(p: bmap->eq[i] + 1 + j, len: total - j);
5334	if (k < 0)
5335	continue;
5336	j += k;
5337	if (!isl_int_is_one(bmap->eq[i][1 + j]) &&
5338	!isl_int_is_negone(bmap->eq[i][1 + j]))
5339	return isl_bool_true;
5340
5341	j += 1;
5342	k = isl_seq_first_non_zero(p: bmap->eq[i] + 1 + j, len: total - j);
5343	if (k >= 0)
5344	return isl_bool_true;
5345	}
5346
5347	return isl_bool_false;
5348	}
5349
5350	/* Remove any common factor g from the constraint coefficients in "v".
5351	* The constant term is stored in the first position and is replaced
5352	* by floor(c/g). If any common factor is removed and if this results
5353	* in a tightening of the constraint, then set *tightened.
5354	*/
5355	static __isl_give isl_vec *normalize_constraint(__isl_take isl_vec *v,
5356	int *tightened)
5357	{
5358	isl_ctx *ctx;
5359
5360	if (!v)
5361	return NULL;
5362	ctx = isl_vec_get_ctx(vec: v);
5363	isl_seq_gcd(p: v->el + 1, len: v->size - 1, gcd: &ctx->normalize_gcd);
5364	if (isl_int_is_zero(ctx->normalize_gcd))
5365	return v;
5366	if (isl_int_is_one(ctx->normalize_gcd))
5367	return v;
5368	v = isl_vec_cow(vec: v);
5369	if (!v)
5370	return NULL;
5371	if (tightened && !isl_int_is_divisible_by(v->el[0], ctx->normalize_gcd))
5372	*tightened = 1;
5373	isl_int_fdiv_q(v->el[0], v->el[0], ctx->normalize_gcd);
5374	isl_seq_scale_down(dst: v->el + 1, src: v->el + 1, f: ctx->normalize_gcd,
5375	len: v->size - 1);
5376	return v;
5377	}
5378
5379	/* If "bmap" is an integer set that satisfies any equality involving
5380	* more than 2 variables and/or has coefficients different from -1 and 1,
5381	* then use variable compression to reduce the coefficients by removing
5382	* any (hidden) common factor.
5383	* In particular, apply the variable compression to each constraint,
5384	* factor out any common factor in the non-constant coefficients and
5385	* then apply the inverse of the compression.
5386	* At the end, we mark the basic map as having reduced constants.
5387	* If this flag is still set on the next invocation of this function,
5388	* then we skip the computation.
5389	*
5390	* Removing a common factor may result in a tightening of some of
5391	* the constraints. If this happens, then we may end up with two
5392	* opposite inequalities that can be replaced by an equality.
5393	* We therefore call isl_basic_map_detect_inequality_pairs,
5394	* which checks for such pairs of inequalities as well as eliminate_divs_eq
5395	* and isl_basic_map_gauss if such a pair was found.
5396	*
5397	* Tightening may also result in some other constraints becoming
5398	* (rationally) redundant with respect to the tightened constraint
5399	* (in combination with other constraints). The basic map may
5400	* therefore no longer be assumed to have no redundant constraints.
5401	*
5402	* Note that this function may leave the result in an inconsistent state.
5403	* In particular, the constraints may not be gaussed.
5404	* Unfortunately, isl_map_coalesce actually depends on this inconsistent state
5405	* for some of the test cases to pass successfully.
5406	* Any potential modification of the representation is therefore only
5407	* performed on a single copy of the basic map.
5408	*/
5409	__isl_give isl_basic_map *isl_basic_map_reduce_coefficients(
5410	__isl_take isl_basic_map *bmap)
5411	{
5412	isl_size total;
5413	isl_bool multi;
5414	isl_ctx *ctx;
5415	isl_vec *v;
5416	isl_mat *eq, *T, *T2;
5417	int i;
5418	int tightened;
5419
5420	if (!bmap)
5421	return NULL;
5422	if (ISL_F_ISSET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS))
5423	return bmap;
5424	if (isl_basic_map_is_rational(bmap))
5425	return bmap;
5426	if (bmap->n_eq == 0)
5427	return bmap;
5428	multi = has_multiple_var_equality(bmap);
5429	if (multi < 0)
5430	return isl_basic_map_free(bmap);
5431	if (!multi)
5432	return bmap;
5433
5434	total = isl_basic_map_dim(bmap, type: isl_dim_all);
5435	if (total < 0)
5436	return isl_basic_map_free(bmap);
5437	ctx = isl_basic_map_get_ctx(bmap);
5438	v = isl_vec_alloc(ctx, size: 1 + total);
5439	if (!v)
5440	return isl_basic_map_free(bmap);
5441
5442	eq = isl_mat_sub_alloc6(ctx, row: bmap->eq, first_row: 0, n_row: bmap->n_eq, first_col: 0, n_col: 1 + total);
5443	T = isl_mat_variable_compression(B: eq, T2: &T2);
5444	if (!T || !T2)
5445	goto error;
5446	if (T->n_col == 0) {
5447	isl_mat_free(mat: T);
5448	isl_mat_free(mat: T2);
5449	isl_vec_free(vec: v);
5450	return isl_basic_map_set_to_empty(bmap);
5451	}
5452
5453	bmap = isl_basic_map_cow(bmap);
5454	if (!bmap)
5455	goto error;
5456
5457	tightened = 0;
5458	for (i = 0; i < bmap->n_ineq; ++i) {
5459	isl_seq_cpy(dst: v->el, src: bmap->ineq[i], len: 1 + total);
5460	v = isl_vec_mat_product(vec: v, mat: isl_mat_copy(mat: T));
5461	v = normalize_constraint(v, tightened: &tightened);
5462	v = isl_vec_mat_product(vec: v, mat: isl_mat_copy(mat: T2));
5463	if (!v)
5464	goto error;
5465	isl_seq_cpy(dst: bmap->ineq[i], src: v->el, len: 1 + total);
5466	}
5467
5468	isl_mat_free(mat: T);
5469	isl_mat_free(mat: T2);
5470	isl_vec_free(vec: v);
5471
5472	ISL_F_SET(bmap, ISL_BASIC_MAP_REDUCED_COEFFICIENTS);
5473
5474	if (tightened) {
5475	int progress = 0;
5476
5477	ISL_F_CLR(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
5478	bmap = isl_basic_map_detect_inequality_pairs(bmap, progress: &progress);
5479	if (progress) {
5480	bmap = eliminate_divs_eq(bmap, progress: &progress);
5481	bmap = isl_basic_map_gauss(bmap, NULL);
5482	}
5483	}
5484
5485	return bmap;
5486	error:
5487	isl_mat_free(mat: T);
5488	isl_mat_free(mat: T2);
5489	isl_vec_free(vec: v);
5490	return isl_basic_map_free(bmap);
5491	}
5492
5493	/* Shift the integer division at position "div" of "bmap"
5494	* by "shift" times the variable at position "pos".
5495	* "pos" is as determined by isl_basic_map_offset, i.e., pos == 0
5496	* corresponds to the constant term.
5497	*
5498	* That is, if the integer division has the form
5499	*
5500	* floor(f(x)/d)
5501	*
5502	* then replace it by
5503	*
5504	* floor((f(x) + shift * d * x_pos)/d) - shift * x_pos
5505	*/
5506	__isl_give isl_basic_map *isl_basic_map_shift_div(
5507	__isl_take isl_basic_map *bmap, int div, int pos, isl_int shift)
5508	{
5509	int i;
5510	isl_size total, n_div;
5511
5512	if (isl_int_is_zero(shift))
5513	return bmap;
5514	total = isl_basic_map_dim(bmap, type: isl_dim_all);
5515	n_div = isl_basic_map_dim(bmap, type: isl_dim_div);
5516	total -= n_div;
5517	if (total < 0 || n_div < 0)
5518	return isl_basic_map_free(bmap);
5519
5520	isl_int_addmul(bmap->div[div][1 + pos], shift, bmap->div[div][0]);
5521
5522	for (i = 0; i < bmap->n_eq; ++i) {
5523	if (isl_int_is_zero(bmap->eq[i][1 + total + div]))
5524	continue;
5525	isl_int_submul(bmap->eq[i][pos],
5526	shift, bmap->eq[i][1 + total + div]);
5527	}
5528	for (i = 0; i < bmap->n_ineq; ++i) {
5529	if (isl_int_is_zero(bmap->ineq[i][1 + total + div]))
5530	continue;
5531	isl_int_submul(bmap->ineq[i][pos],
5532	shift, bmap->ineq[i][1 + total + div]);
5533	}
5534	for (i = 0; i < bmap->n_div; ++i) {
5535	if (isl_int_is_zero(bmap->div[i][0]))
5536	continue;
5537	if (isl_int_is_zero(bmap->div[i][1 + 1 + total + div]))
5538	continue;
5539	isl_int_submul(bmap->div[i][1 + pos],
5540	shift, bmap->div[i][1 + 1 + total + div]);
5541	}
5542
5543	return bmap;
5544	}
5545