1#include <isl_ctx_private.h>
2#include <isl/constraint.h>
3#include <isl/set.h>
4#include <isl_polynomial_private.h>
5#include <isl_morph.h>
6#include <isl_range.h>
7
8struct range_data {
9	struct isl_bound	*bound;
10	int 		    	*signs;
11	int			sign;
12	int			test_monotonicity;
13	int		    	monotonicity;
14	int			tight;
15	isl_qpolynomial	    	*poly;
16	isl_pw_qpolynomial_fold *pwf;
17	isl_pw_qpolynomial_fold *pwf_tight;
18};
19
20static int propagate_on_domain(__isl_take isl_basic_set *bset,
21	__isl_take isl_qpolynomial *poly, struct range_data *data);
22
23/* Check whether the polynomial "poly" has sign "sign" over "bset",
24 * i.e., if sign == 1, check that the lower bound on the polynomial
25 * is non-negative and if sign == -1, check that the upper bound on
26 * the polynomial is non-positive.
27 */
28static int has_sign(__isl_keep isl_basic_set *bset,
29	__isl_keep isl_qpolynomial *poly, int sign, int *signs)
30{
31	struct range_data data_m;
32	unsigned nvar;
33	unsigned nparam;
34	isl_space *dim;
35	isl_qpolynomial *opt;
36	int r;
37	enum isl_fold type;
38
39	nparam = isl_basic_set_dim(bset, isl_dim_param);
40	nvar = isl_basic_set_dim(bset, isl_dim_set);
41
42	bset = isl_basic_set_copy(bset);
43	poly = isl_qpolynomial_copy(poly);
44
45	bset = isl_basic_set_move_dims(bset, isl_dim_set, 0,
46					isl_dim_param, 0, nparam);
47	poly = isl_qpolynomial_move_dims(poly, isl_dim_in, 0,
48					isl_dim_param, 0, nparam);
49
50	dim = isl_qpolynomial_get_space(poly);
51	dim = isl_space_params(dim);
52	dim = isl_space_from_domain(dim);
53	dim = isl_space_add_dims(dim, isl_dim_out, 1);
54
55	data_m.test_monotonicity = 0;
56	data_m.signs = signs;
57	data_m.sign = -sign;
58	type = data_m.sign < 0 ? isl_fold_min : isl_fold_max;
59	data_m.pwf = isl_pw_qpolynomial_fold_zero(dim, type);
60	data_m.tight = 0;
61	data_m.pwf_tight = NULL;
62
63	if (propagate_on_domain(bset, poly, &data_m) < 0)
64		goto error;
65
66	if (sign > 0)
67		opt = isl_pw_qpolynomial_fold_min(data_m.pwf);
68	else
69		opt = isl_pw_qpolynomial_fold_max(data_m.pwf);
70
71	if (!opt)
72		r = -1;
73	else if (isl_qpolynomial_is_nan(opt) ||
74		 isl_qpolynomial_is_infty(opt) ||
75		 isl_qpolynomial_is_neginfty(opt))
76		r = 0;
77	else
78		r = sign * isl_qpolynomial_sgn(opt) >= 0;
79
80	isl_qpolynomial_free(opt);
81
82	return r;
83error:
84	isl_pw_qpolynomial_fold_free(data_m.pwf);
85	return -1;
86}
87
88/* Return  1 if poly is monotonically increasing in the last set variable,
89 *        -1 if poly is monotonically decreasing in the last set variable,
90 *	   0 if no conclusion,
91 *	  -2 on error.
92 *
93 * We simply check the sign of p(x+1)-p(x)
94 */
95static int monotonicity(__isl_keep isl_basic_set *bset,
96	__isl_keep isl_qpolynomial *poly, struct range_data *data)
97{
98	isl_ctx *ctx;
99	isl_space *dim;
100	isl_qpolynomial *sub = NULL;
101	isl_qpolynomial *diff = NULL;
102	int result = 0;
103	int s;
104	unsigned nvar;
105
106	ctx = isl_qpolynomial_get_ctx(poly);
107	dim = isl_qpolynomial_get_domain_space(poly);
108
109	nvar = isl_basic_set_dim(bset, isl_dim_set);
110
111	sub = isl_qpolynomial_var_on_domain(isl_space_copy(dim), isl_dim_set, nvar - 1);
112	sub = isl_qpolynomial_add(sub,
113		isl_qpolynomial_rat_cst_on_domain(dim, ctx->one, ctx->one));
114
115	diff = isl_qpolynomial_substitute(isl_qpolynomial_copy(poly),
116			isl_dim_in, nvar - 1, 1, &sub);
117	diff = isl_qpolynomial_sub(diff, isl_qpolynomial_copy(poly));
118
119	s = has_sign(bset, diff, 1, data->signs);
120	if (s < 0)
121		goto error;
122	if (s)
123		result = 1;
124	else {
125		s = has_sign(bset, diff, -1, data->signs);
126		if (s < 0)
127			goto error;
128		if (s)
129			result = -1;
130	}
131
132	isl_qpolynomial_free(diff);
133	isl_qpolynomial_free(sub);
134
135	return result;
136error:
137	isl_qpolynomial_free(diff);
138	isl_qpolynomial_free(sub);
139	return -2;
140}
141
142static __isl_give isl_qpolynomial *bound2poly(__isl_take isl_constraint *bound,
143	__isl_take isl_space *dim, unsigned pos, int sign)
144{
145	if (!bound) {
146		if (sign > 0)
147			return isl_qpolynomial_infty_on_domain(dim);
148		else
149			return isl_qpolynomial_neginfty_on_domain(dim);
150	}
151	isl_space_free(dim);
152	return isl_qpolynomial_from_constraint(bound, isl_dim_set, pos);
153}
154
155static int bound_is_integer(__isl_take isl_constraint *bound, unsigned pos)
156{
157	isl_int c;
158	int is_int;
159
160	if (!bound)
161		return 1;
162
163	isl_int_init(c);
164	isl_constraint_get_coefficient(bound, isl_dim_set, pos, &c);
165	is_int = isl_int_is_one(c) || isl_int_is_negone(c);
166	isl_int_clear(c);
167
168	return is_int;
169}
170
171struct isl_fixed_sign_data {
172	int		*signs;
173	int		sign;
174	isl_qpolynomial	*poly;
175};
176
177/* Add term "term" to data->poly if it has sign data->sign.
178 * The sign is determined based on the signs of the parameters
179 * and variables in data->signs.  The integer divisions, if
180 * any, are assumed to be non-negative.
181 */
182static int collect_fixed_sign_terms(__isl_take isl_term *term, void *user)
183{
184	struct isl_fixed_sign_data *data = (struct isl_fixed_sign_data *)user;
185	isl_int n;
186	int i;
187	int sign;
188	unsigned nparam;
189	unsigned nvar;
190
191	if (!term)
192		return -1;
193
194	nparam = isl_term_dim(term, isl_dim_param);
195	nvar = isl_term_dim(term, isl_dim_set);
196
197	isl_int_init(n);
198
199	isl_term_get_num(term, &n);
200
201	sign = isl_int_sgn(n);
202	for (i = 0; i < nparam; ++i) {
203		if (data->signs[i] > 0)
204			continue;
205		if (isl_term_get_exp(term, isl_dim_param, i) % 2)
206			sign = -sign;
207	}
208	for (i = 0; i < nvar; ++i) {
209		if (data->signs[nparam + i] > 0)
210			continue;
211		if (isl_term_get_exp(term, isl_dim_set, i) % 2)
212			sign = -sign;
213	}
214
215	if (sign == data->sign) {
216		isl_qpolynomial *t = isl_qpolynomial_from_term(term);
217
218		data->poly = isl_qpolynomial_add(data->poly, t);
219	} else
220		isl_term_free(term);
221
222	isl_int_clear(n);
223
224	return 0;
225}
226
227/* Construct and return a polynomial that consists of the terms
228 * in "poly" that have sign "sign".  The integer divisions, if
229 * any, are assumed to be non-negative.
230 */
231__isl_give isl_qpolynomial *isl_qpolynomial_terms_of_sign(
232	__isl_keep isl_qpolynomial *poly, int *signs, int sign)
233{
234	isl_space *space;
235	struct isl_fixed_sign_data data = { signs, sign };
236
237	space = isl_qpolynomial_get_domain_space(poly);
238	data.poly = isl_qpolynomial_zero_on_domain(space);
239
240	if (isl_qpolynomial_foreach_term(poly, collect_fixed_sign_terms, &data) < 0)
241		goto error;
242
243	return data.poly;
244error:
245	isl_qpolynomial_free(data.poly);
246	return NULL;
247}
248
249/* Helper function to add a guarded polynomial to either pwf_tight or pwf,
250 * depending on whether the result has been determined to be tight.
251 */
252static int add_guarded_poly(__isl_take isl_basic_set *bset,
253	__isl_take isl_qpolynomial *poly, struct range_data *data)
254{
255	enum isl_fold type = data->sign < 0 ? isl_fold_min : isl_fold_max;
256	isl_set *set;
257	isl_qpolynomial_fold *fold;
258	isl_pw_qpolynomial_fold *pwf;
259
260	bset = isl_basic_set_params(bset);
261	poly = isl_qpolynomial_project_domain_on_params(poly);
262
263	fold = isl_qpolynomial_fold_alloc(type, poly);
264	set = isl_set_from_basic_set(bset);
265	pwf = isl_pw_qpolynomial_fold_alloc(type, set, fold);
266	if (data->tight)
267		data->pwf_tight = isl_pw_qpolynomial_fold_fold(
268						data->pwf_tight, pwf);
269	else
270		data->pwf = isl_pw_qpolynomial_fold_fold(data->pwf, pwf);
271
272	return 0;
273}
274
275/* Given a lower and upper bound on the final variable and constraints
276 * on the remaining variables where these bounds are active,
277 * eliminate the variable from data->poly based on these bounds.
278 * If the polynomial has been determined to be monotonic
279 * in the variable, then simply plug in the appropriate bound.
280 * If the current polynomial is tight and if this bound is integer,
281 * then the result is still tight.  In all other cases, the results
282 * may not be tight.
283 * Otherwise, plug in the largest bound (in absolute value) in
284 * the positive terms (if an upper bound is wanted) or the negative terms
285 * (if a lower bounded is wanted) and the other bound in the other terms.
286 *
287 * If all variables have been eliminated, then record the result.
288 * Ohterwise, recurse on the next variable.
289 */
290static int propagate_on_bound_pair(__isl_take isl_constraint *lower,
291	__isl_take isl_constraint *upper, __isl_take isl_basic_set *bset,
292	void *user)
293{
294	struct range_data *data = (struct range_data *)user;
295	int save_tight = data->tight;
296	isl_qpolynomial *poly;
297	int r;
298	unsigned nvar;
299
300	nvar = isl_basic_set_dim(bset, isl_dim_set);
301
302	if (data->monotonicity) {
303		isl_qpolynomial *sub;
304		isl_space *dim = isl_qpolynomial_get_domain_space(data->poly);
305		if (data->monotonicity * data->sign > 0) {
306			if (data->tight)
307				data->tight = bound_is_integer(upper, nvar);
308			sub = bound2poly(upper, dim, nvar, 1);
309			isl_constraint_free(lower);
310		} else {
311			if (data->tight)
312				data->tight = bound_is_integer(lower, nvar);
313			sub = bound2poly(lower, dim, nvar, -1);
314			isl_constraint_free(upper);
315		}
316		poly = isl_qpolynomial_copy(data->poly);
317		poly = isl_qpolynomial_substitute(poly, isl_dim_in, nvar, 1, &sub);
318		poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, nvar, 1);
319
320		isl_qpolynomial_free(sub);
321	} else {
322		isl_qpolynomial *l, *u;
323		isl_qpolynomial *pos, *neg;
324		isl_space *dim = isl_qpolynomial_get_domain_space(data->poly);
325		unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
326		int sign = data->sign * data->signs[nparam + nvar];
327
328		data->tight = 0;
329
330		u = bound2poly(upper, isl_space_copy(dim), nvar, 1);
331		l = bound2poly(lower, dim, nvar, -1);
332
333		pos = isl_qpolynomial_terms_of_sign(data->poly, data->signs, sign);
334		neg = isl_qpolynomial_terms_of_sign(data->poly, data->signs, -sign);
335
336		pos = isl_qpolynomial_substitute(pos, isl_dim_in, nvar, 1, &u);
337		neg = isl_qpolynomial_substitute(neg, isl_dim_in, nvar, 1, &l);
338
339		poly = isl_qpolynomial_add(pos, neg);
340		poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, nvar, 1);
341
342		isl_qpolynomial_free(u);
343		isl_qpolynomial_free(l);
344	}
345
346	if (isl_basic_set_dim(bset, isl_dim_set) == 0)
347		r = add_guarded_poly(bset, poly, data);
348	else
349		r = propagate_on_domain(bset, poly, data);
350
351	data->tight = save_tight;
352
353	return r;
354}
355
356/* Recursively perform range propagation on the polynomial "poly"
357 * defined over the basic set "bset" and collect the results in "data".
358 */
359static int propagate_on_domain(__isl_take isl_basic_set *bset,
360	__isl_take isl_qpolynomial *poly, struct range_data *data)
361{
362	isl_ctx *ctx;
363	isl_qpolynomial *save_poly = data->poly;
364	int save_monotonicity = data->monotonicity;
365	unsigned d;
366
367	if (!bset || !poly)
368		goto error;
369
370	ctx = isl_basic_set_get_ctx(bset);
371	d = isl_basic_set_dim(bset, isl_dim_set);
372	isl_assert(ctx, d >= 1, goto error);
373
374	if (isl_qpolynomial_is_cst(poly, NULL, NULL)) {
375		bset = isl_basic_set_project_out(bset, isl_dim_set, 0, d);
376		poly = isl_qpolynomial_drop_dims(poly, isl_dim_in, 0, d);
377		return add_guarded_poly(bset, poly, data);
378	}
379
380	if (data->test_monotonicity)
381		data->monotonicity = monotonicity(bset, poly, data);
382	else
383		data->monotonicity = 0;
384	if (data->monotonicity < -1)
385		goto error;
386
387	data->poly = poly;
388	if (isl_basic_set_foreach_bound_pair(bset, isl_dim_set, d - 1,
389					    &propagate_on_bound_pair, data) < 0)
390		goto error;
391
392	isl_basic_set_free(bset);
393	isl_qpolynomial_free(poly);
394	data->monotonicity = save_monotonicity;
395	data->poly = save_poly;
396
397	return 0;
398error:
399	isl_basic_set_free(bset);
400	isl_qpolynomial_free(poly);
401	data->monotonicity = save_monotonicity;
402	data->poly = save_poly;
403	return -1;
404}
405
406static int basic_guarded_poly_bound(__isl_take isl_basic_set *bset, void *user)
407{
408	struct range_data *data = (struct range_data *)user;
409	isl_ctx *ctx;
410	unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
411	unsigned dim = isl_basic_set_dim(bset, isl_dim_set);
412	int r;
413
414	data->signs = NULL;
415
416	ctx = isl_basic_set_get_ctx(bset);
417	data->signs = isl_alloc_array(ctx, int,
418					isl_basic_set_dim(bset, isl_dim_all));
419
420	if (isl_basic_set_dims_get_sign(bset, isl_dim_set, 0, dim,
421					data->signs + nparam) < 0)
422		goto error;
423	if (isl_basic_set_dims_get_sign(bset, isl_dim_param, 0, nparam,
424					data->signs) < 0)
425		goto error;
426
427	r = propagate_on_domain(bset, isl_qpolynomial_copy(data->poly), data);
428
429	free(data->signs);
430
431	return r;
432error:
433	free(data->signs);
434	isl_basic_set_free(bset);
435	return -1;
436}
437
438static int qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
439	__isl_take isl_qpolynomial *poly, struct range_data *data)
440{
441	unsigned nparam = isl_basic_set_dim(bset, isl_dim_param);
442	unsigned nvar = isl_basic_set_dim(bset, isl_dim_set);
443	isl_set *set = NULL;
444
445	if (!bset)
446		goto error;
447
448	if (nvar == 0)
449		return add_guarded_poly(bset, poly, data);
450
451	set = isl_set_from_basic_set(bset);
452	set = isl_set_split_dims(set, isl_dim_param, 0, nparam);
453	set = isl_set_split_dims(set, isl_dim_set, 0, nvar);
454
455	data->poly = poly;
456
457	data->test_monotonicity = 1;
458	if (isl_set_foreach_basic_set(set, &basic_guarded_poly_bound, data) < 0)
459		goto error;
460
461	isl_set_free(set);
462	isl_qpolynomial_free(poly);
463
464	return 0;
465error:
466	isl_set_free(set);
467	isl_qpolynomial_free(poly);
468	return -1;
469}
470
471int isl_qpolynomial_bound_on_domain_range(__isl_take isl_basic_set *bset,
472	__isl_take isl_qpolynomial *poly, struct isl_bound *bound)
473{
474	struct range_data data;
475	int r;
476
477	data.pwf = bound->pwf;
478	data.pwf_tight = bound->pwf_tight;
479	data.tight = bound->check_tight;
480	if (bound->type == isl_fold_min)
481		data.sign = -1;
482	else
483		data.sign = 1;
484
485	r = qpolynomial_bound_on_domain_range(bset, poly, &data);
486
487	bound->pwf = data.pwf;
488	bound->pwf_tight = data.pwf_tight;
489
490	return r;
491}
492