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