1/*
2 * Copyright 2010      INRIA Saclay
3 *
4 * Use of this software is governed by the MIT license
5 *
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
10
11#include <isl_map_private.h>
12#include <isl/set.h>
13#include <isl_space_private.h>
14#include <isl/seq.h>
15
16/*
17 * Let C be a cone and define
18 *
19 *	C' := { y | forall x in C : y x >= 0 }
20 *
21 * C' contains the coefficients of all linear constraints
22 * that are valid for C.
23 * Furthermore, C'' = C.
24 *
25 * If C is defined as { x | A x >= 0 }
26 * then any element in C' must be a non-negative combination
27 * of the rows of A, i.e., y = t A with t >= 0.  That is,
28 *
29 *	C' = { y | exists t >= 0 : y = t A }
30 *
31 * If any of the rows in A actually represents an equality, then
32 * also negative combinations of this row are allowed and so the
33 * non-negativity constraint on the corresponding element of t
34 * can be dropped.
35 *
36 * A polyhedron P = { x | b + A x >= 0 } can be represented
37 * in homogeneous coordinates by the cone
38 * C = { [z,x] | b z + A x >= and z >= 0 }
39 * The valid linear constraints on C correspond to the valid affine
40 * constraints on P.
41 * This is essentially Farkas' lemma.
42 *
43 * Let A' = [b A], then, since
44 *				  [ 1 0 ]
45 *		[ w y ] = [t_0 t] [ b A ]
46 *
47 * we have
48 *
49 *	C' = { w, y | exists t_0, t >= 0 : y = t A' and w = t_0 + t b }
50 * or
51 *
52 *	C' = { w, y | exists t >= 0 : y = t A' and w - t b >= 0 }
53 *
54 * In practice, we introduce an extra variable (w), shifting all
55 * other variables to the right, and an extra inequality
56 * (w - t b >= 0) corresponding to the positivity constraint on
57 * the homogeneous coordinate.
58 *
59 * When going back from coefficients to solutions, we immediately
60 * plug in 1 for z, which corresponds to shifting all variables
61 * to the left, with the leftmost ending up in the constant position.
62 */
63
64/* Add the given prefix to all named isl_dim_set dimensions in "dim".
65 */
66static __isl_give isl_space *isl_space_prefix(__isl_take isl_space *dim,
67	const char *prefix)
68{
69	int i;
70	isl_ctx *ctx;
71	unsigned nvar;
72	size_t prefix_len = strlen(prefix);
73
74	if (!dim)
75		return NULL;
76
77	ctx = isl_space_get_ctx(dim);
78	nvar = isl_space_dim(dim, isl_dim_set);
79
80	for (i = 0; i < nvar; ++i) {
81		const char *name;
82		char *prefix_name;
83
84		name = isl_space_get_dim_name(dim, isl_dim_set, i);
85		if (!name)
86			continue;
87
88		prefix_name = isl_alloc_array(ctx, char,
89					      prefix_len + strlen(name) + 1);
90		if (!prefix_name)
91			goto error;
92		memcpy(prefix_name, prefix, prefix_len);
93		strcpy(prefix_name + prefix_len, name);
94
95		dim = isl_space_set_dim_name(dim, isl_dim_set, i, prefix_name);
96		free(prefix_name);
97	}
98
99	return dim;
100error:
101	isl_space_free(dim);
102	return NULL;
103}
104
105/* Given a dimension specification of the solutions space, construct
106 * a dimension specification for the space of coefficients.
107 *
108 * In particular transform
109 *
110 *	[params] -> { S }
111 *
112 * to
113 *
114 *	{ coefficients[[cst, params] -> S] }
115 *
116 * and prefix each dimension name with "c_".
117 */
118static __isl_give isl_space *isl_space_coefficients(__isl_take isl_space *dim)
119{
120	isl_space *dim_param;
121	unsigned nvar;
122	unsigned nparam;
123
124	nvar = isl_space_dim(dim, isl_dim_set);
125	nparam = isl_space_dim(dim, isl_dim_param);
126	dim_param = isl_space_copy(dim);
127	dim_param = isl_space_drop_dims(dim_param, isl_dim_set, 0, nvar);
128	dim_param = isl_space_move_dims(dim_param, isl_dim_set, 0,
129				 isl_dim_param, 0, nparam);
130	dim_param = isl_space_prefix(dim_param, "c_");
131	dim_param = isl_space_insert_dims(dim_param, isl_dim_set, 0, 1);
132	dim_param = isl_space_set_dim_name(dim_param, isl_dim_set, 0, "c_cst");
133	dim = isl_space_drop_dims(dim, isl_dim_param, 0, nparam);
134	dim = isl_space_prefix(dim, "c_");
135	dim = isl_space_join(isl_space_from_domain(dim_param),
136			   isl_space_from_range(dim));
137	dim = isl_space_wrap(dim);
138	dim = isl_space_set_tuple_name(dim, isl_dim_set, "coefficients");
139
140	return dim;
141}
142
143/* Drop the given prefix from all named dimensions of type "type" in "dim".
144 */
145static __isl_give isl_space *isl_space_unprefix(__isl_take isl_space *dim,
146	enum isl_dim_type type, const char *prefix)
147{
148	int i;
149	unsigned n;
150	size_t prefix_len = strlen(prefix);
151
152	n = isl_space_dim(dim, type);
153
154	for (i = 0; i < n; ++i) {
155		const char *name;
156
157		name = isl_space_get_dim_name(dim, type, i);
158		if (!name)
159			continue;
160		if (strncmp(name, prefix, prefix_len))
161			continue;
162
163		dim = isl_space_set_dim_name(dim, type, i, name + prefix_len);
164	}
165
166	return dim;
167}
168
169/* Given a dimension specification of the space of coefficients, construct
170 * a dimension specification for the space of solutions.
171 *
172 * In particular transform
173 *
174 *	{ coefficients[[cst, params] -> S] }
175 *
176 * to
177 *
178 *	[params] -> { S }
179 *
180 * and drop the "c_" prefix from the dimension names.
181 */
182static __isl_give isl_space *isl_space_solutions(__isl_take isl_space *dim)
183{
184	unsigned nparam;
185
186	dim = isl_space_unwrap(dim);
187	dim = isl_space_drop_dims(dim, isl_dim_in, 0, 1);
188	dim = isl_space_unprefix(dim, isl_dim_in, "c_");
189	dim = isl_space_unprefix(dim, isl_dim_out, "c_");
190	nparam = isl_space_dim(dim, isl_dim_in);
191	dim = isl_space_move_dims(dim, isl_dim_param, 0, isl_dim_in, 0, nparam);
192	dim = isl_space_range(dim);
193
194	return dim;
195}
196
197/* Compute the dual of "bset" by applying Farkas' lemma.
198 * As explained above, we add an extra dimension to represent
199 * the coefficient of the constant term when going from solutions
200 * to coefficients (shift == 1) and we drop the extra dimension when going
201 * in the opposite direction (shift == -1).  "dim" is the space in which
202 * the dual should be created.
203 */
204static __isl_give isl_basic_set *farkas(__isl_take isl_space *dim,
205	__isl_take isl_basic_set *bset, int shift)
206{
207	int i, j, k;
208	isl_basic_set *dual = NULL;
209	unsigned total;
210
211	total = isl_basic_set_total_dim(bset);
212
213	dual = isl_basic_set_alloc_space(dim, bset->n_eq + bset->n_ineq,
214					total, bset->n_ineq + (shift > 0));
215	dual = isl_basic_set_set_rational(dual);
216
217	for (i = 0; i < bset->n_eq + bset->n_ineq; ++i) {
218		k = isl_basic_set_alloc_div(dual);
219		if (k < 0)
220			goto error;
221		isl_int_set_si(dual->div[k][0], 0);
222	}
223
224	for (i = 0; i < total; ++i) {
225		k = isl_basic_set_alloc_equality(dual);
226		if (k < 0)
227			goto error;
228		isl_seq_clr(dual->eq[k], 1 + shift + total);
229		isl_int_set_si(dual->eq[k][1 + shift + i], -1);
230		for (j = 0; j < bset->n_eq; ++j)
231			isl_int_set(dual->eq[k][1 + shift + total + j],
232				    bset->eq[j][1 + i]);
233		for (j = 0; j < bset->n_ineq; ++j)
234			isl_int_set(dual->eq[k][1 + shift + total + bset->n_eq + j],
235				    bset->ineq[j][1 + i]);
236	}
237
238	for (i = 0; i < bset->n_ineq; ++i) {
239		k = isl_basic_set_alloc_inequality(dual);
240		if (k < 0)
241			goto error;
242		isl_seq_clr(dual->ineq[k],
243			    1 + shift + total + bset->n_eq + bset->n_ineq);
244		isl_int_set_si(dual->ineq[k][1 + shift + total + bset->n_eq + i], 1);
245	}
246
247	if (shift > 0) {
248		k = isl_basic_set_alloc_inequality(dual);
249		if (k < 0)
250			goto error;
251		isl_seq_clr(dual->ineq[k], 2 + total);
252		isl_int_set_si(dual->ineq[k][1], 1);
253		for (j = 0; j < bset->n_eq; ++j)
254			isl_int_neg(dual->ineq[k][2 + total + j],
255				    bset->eq[j][0]);
256		for (j = 0; j < bset->n_ineq; ++j)
257			isl_int_neg(dual->ineq[k][2 + total + bset->n_eq + j],
258				    bset->ineq[j][0]);
259	}
260
261	dual = isl_basic_set_remove_divs(dual);
262	isl_basic_set_simplify(dual);
263	isl_basic_set_finalize(dual);
264
265	isl_basic_set_free(bset);
266	return dual;
267error:
268	isl_basic_set_free(bset);
269	isl_basic_set_free(dual);
270	return NULL;
271}
272
273/* Construct a basic set containing the tuples of coefficients of all
274 * valid affine constraints on the given basic set.
275 */
276__isl_give isl_basic_set *isl_basic_set_coefficients(
277	__isl_take isl_basic_set *bset)
278{
279	isl_space *dim;
280
281	if (!bset)
282		return NULL;
283	if (bset->n_div)
284		isl_die(bset->ctx, isl_error_invalid,
285			"input set not allowed to have local variables",
286			goto error);
287
288	dim = isl_basic_set_get_space(bset);
289	dim = isl_space_coefficients(dim);
290
291	return farkas(dim, bset, 1);
292error:
293	isl_basic_set_free(bset);
294	return NULL;
295}
296
297/* Construct a basic set containing the elements that satisfy all
298 * affine constraints whose coefficient tuples are
299 * contained in the given basic set.
300 */
301__isl_give isl_basic_set *isl_basic_set_solutions(
302	__isl_take isl_basic_set *bset)
303{
304	isl_space *dim;
305
306	if (!bset)
307		return NULL;
308	if (bset->n_div)
309		isl_die(bset->ctx, isl_error_invalid,
310			"input set not allowed to have local variables",
311			goto error);
312
313	dim = isl_basic_set_get_space(bset);
314	dim = isl_space_solutions(dim);
315
316	return farkas(dim, bset, -1);
317error:
318	isl_basic_set_free(bset);
319	return NULL;
320}
321
322/* Construct a basic set containing the tuples of coefficients of all
323 * valid affine constraints on the given set.
324 */
325__isl_give isl_basic_set *isl_set_coefficients(__isl_take isl_set *set)
326{
327	int i;
328	isl_basic_set *coeff;
329
330	if (!set)
331		return NULL;
332	if (set->n == 0) {
333		isl_space *dim = isl_set_get_space(set);
334		dim = isl_space_coefficients(dim);
335		coeff = isl_basic_set_universe(dim);
336		coeff = isl_basic_set_set_rational(coeff);
337		isl_set_free(set);
338		return coeff;
339	}
340
341	coeff = isl_basic_set_coefficients(isl_basic_set_copy(set->p[0]));
342
343	for (i = 1; i < set->n; ++i) {
344		isl_basic_set *bset, *coeff_i;
345		bset = isl_basic_set_copy(set->p[i]);
346		coeff_i = isl_basic_set_coefficients(bset);
347		coeff = isl_basic_set_intersect(coeff, coeff_i);
348	}
349
350	isl_set_free(set);
351	return coeff;
352}
353
354/* Construct a basic set containing the elements that satisfy all
355 * affine constraints whose coefficient tuples are
356 * contained in the given set.
357 */
358__isl_give isl_basic_set *isl_set_solutions(__isl_take isl_set *set)
359{
360	int i;
361	isl_basic_set *sol;
362
363	if (!set)
364		return NULL;
365	if (set->n == 0) {
366		isl_space *dim = isl_set_get_space(set);
367		dim = isl_space_solutions(dim);
368		sol = isl_basic_set_universe(dim);
369		sol = isl_basic_set_set_rational(sol);
370		isl_set_free(set);
371		return sol;
372	}
373
374	sol = isl_basic_set_solutions(isl_basic_set_copy(set->p[0]));
375
376	for (i = 1; i < set->n; ++i) {
377		isl_basic_set *bset, *sol_i;
378		bset = isl_basic_set_copy(set->p[i]);
379		sol_i = isl_basic_set_solutions(bset);
380		sol = isl_basic_set_intersect(sol, sol_i);
381	}
382
383	isl_set_free(set);
384	return sol;
385}
386