add isl_qpolynomial_from_constraint
[isl.git] / isl_fold.c
blob036c29684a8590d4086f866f8c8d120c6882c4ad
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
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 */
11 #include <isl_polynomial_private.h>
12 #include <isl_point_private.h>
13 #include <isl_dim.h>
14 #include <isl_map_private.h>
15 #include <isl_lp.h>
16 #include <isl_seq.h>
18 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
19 enum isl_fold type, __isl_take isl_dim *dim, int n)
21 isl_qpolynomial_fold *fold;
23 if (!dim)
24 goto error;
26 isl_assert(dim->ctx, n >= 0, goto error);
27 fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
28 sizeof(struct isl_qpolynomial_fold) +
29 (n - 1) * sizeof(struct isl_qpolynomial *));
30 if (!fold)
31 goto error;
33 fold->ref = 1;
34 fold->size = n;
35 fold->n = 0;
36 fold->type = type;
37 fold->dim = dim;
39 return fold;
40 error:
41 isl_dim_free(dim);
42 return NULL;
45 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
46 enum isl_dim_type type, unsigned first, unsigned n)
48 int i;
50 if (!fold)
51 return -1;
52 if (fold->n == 0 || n == 0)
53 return 0;
55 for (i = 0; i < fold->n; ++i) {
56 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
57 type, first, n);
58 if (involves < 0 || involves)
59 return involves;
61 return 0;
64 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
65 __isl_take isl_qpolynomial_fold *fold,
66 enum isl_dim_type type, unsigned first, unsigned n)
68 int i;
70 if (!fold)
71 return NULL;
72 if (n == 0)
73 return fold;
75 fold = isl_qpolynomial_fold_cow(fold);
76 if (!fold)
77 return NULL;
78 fold->dim = isl_dim_drop(fold->dim, type, first, n);
79 if (!fold->dim)
80 goto error;
82 for (i = 0; i < fold->n; ++i) {
83 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
84 type, first, n);
85 if (!fold->qp[i])
86 goto error;
89 return fold;
90 error:
91 isl_qpolynomial_fold_free(fold);
92 return NULL;
95 static int isl_qpolynomial_cst_sign(__isl_keep isl_qpolynomial *qp)
97 struct isl_upoly_cst *cst;
99 cst = isl_upoly_as_cst(qp->upoly);
100 if (!cst)
101 return 0;
103 return isl_int_sgn(cst->n) < 0 ? -1 : 1;
106 static int isl_qpolynomial_aff_sign(__isl_keep isl_set *set,
107 __isl_keep isl_qpolynomial *qp)
109 enum isl_lp_result res;
110 isl_vec *aff;
111 isl_int opt;
112 int sgn = 0;
114 aff = isl_qpolynomial_extract_affine(qp);
115 if (!aff)
116 return 0;
118 isl_int_init(opt);
120 res = isl_set_solve_lp(set, 0, aff->el + 1, aff->el[0],
121 &opt, NULL, NULL);
122 if (res == isl_lp_error)
123 goto done;
124 if (res == isl_lp_empty ||
125 (res == isl_lp_ok && !isl_int_is_neg(opt))) {
126 sgn = 1;
127 goto done;
130 res = isl_set_solve_lp(set, 1, aff->el + 1, aff->el[0],
131 &opt, NULL, NULL);
132 if (res == isl_lp_ok && !isl_int_is_pos(opt))
133 sgn = -1;
135 done:
136 isl_int_clear(opt);
137 isl_vec_free(aff);
138 return sgn;
141 /* Determine, if possible, the sign of the quasipolynomial "qp" on
142 * the domain "set".
144 * If qp is a constant, then the problem is trivial.
145 * If qp is linear, then we check if the minimum of the corresponding
146 * affine constraint is non-negative or if the maximum is non-positive.
148 * Otherwise, we check if the outermost variable "v" has a lower bound "l"
149 * in "set". If so, we write qp(v,v') as
151 * q(v,v') * (v - l) + r(v')
153 * if q(v,v') and r(v') have the same known sign, then the original
154 * quasipolynomial has the same sign as well.
156 * Return
157 * -1 if qp <= 0
158 * 1 if qp >= 0
159 * 0 if unknown
161 static int isl_qpolynomial_sign(__isl_keep isl_set *set,
162 __isl_keep isl_qpolynomial *qp)
164 int d;
165 int i;
166 int is;
167 struct isl_upoly_rec *rec;
168 isl_vec *v;
169 isl_int l;
170 enum isl_lp_result res;
171 int sgn = 0;
173 is = isl_qpolynomial_is_cst(qp, NULL, NULL);
174 if (is < 0)
175 return 0;
176 if (is)
177 return isl_qpolynomial_cst_sign(qp);
179 is = isl_qpolynomial_is_affine(qp);
180 if (is < 0)
181 return 0;
182 if (is)
183 return isl_qpolynomial_aff_sign(set, qp);
185 if (qp->div->n_row > 0)
186 return 0;
188 rec = isl_upoly_as_rec(qp->upoly);
189 if (!rec)
190 return 0;
192 d = isl_dim_total(qp->dim);
193 v = isl_vec_alloc(set->ctx, 2 + d);
194 if (!v)
195 return 0;
197 isl_seq_clr(v->el + 1, 1 + d);
198 isl_int_set_si(v->el[0], 1);
199 isl_int_set_si(v->el[2 + qp->upoly->var], 1);
201 isl_int_init(l);
203 res = isl_set_solve_lp(set, 0, v->el + 1, v->el[0], &l, NULL, NULL);
204 if (res == isl_lp_ok) {
205 isl_qpolynomial *min;
206 isl_qpolynomial *base;
207 isl_qpolynomial *r, *q;
208 isl_qpolynomial *t;
210 min = isl_qpolynomial_cst(isl_dim_copy(qp->dim), l);
211 base = isl_qpolynomial_pow(isl_dim_copy(qp->dim),
212 qp->upoly->var, 1);
214 r = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
215 isl_upoly_copy(rec->p[rec->n - 1]));
216 q = isl_qpolynomial_copy(r);
218 for (i = rec->n - 2; i >= 0; --i) {
219 r = isl_qpolynomial_mul(r, isl_qpolynomial_copy(min));
220 t = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
221 isl_upoly_copy(rec->p[i]));
222 r = isl_qpolynomial_add(r, t);
223 if (i == 0)
224 break;
225 q = isl_qpolynomial_mul(q, isl_qpolynomial_copy(base));
226 q = isl_qpolynomial_add(q, isl_qpolynomial_copy(r));
229 if (isl_qpolynomial_is_zero(q))
230 sgn = isl_qpolynomial_sign(set, r);
231 else if (isl_qpolynomial_is_zero(r))
232 sgn = isl_qpolynomial_sign(set, q);
233 else {
234 int sgn_q, sgn_r;
235 sgn_r = isl_qpolynomial_sign(set, r);
236 sgn_q = isl_qpolynomial_sign(set, q);
237 if (sgn_r == sgn_q)
238 sgn = sgn_r;
241 isl_qpolynomial_free(min);
242 isl_qpolynomial_free(base);
243 isl_qpolynomial_free(q);
244 isl_qpolynomial_free(r);
247 isl_int_clear(l);
249 isl_vec_free(v);
251 return sgn;
254 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold_on_domain(
255 __isl_keep isl_set *set,
256 __isl_take isl_qpolynomial_fold *fold1,
257 __isl_take isl_qpolynomial_fold *fold2)
259 int i, j;
260 int n1;
261 struct isl_qpolynomial_fold *res = NULL;
262 int better;
264 if (!fold1 || !fold2)
265 goto error;
267 isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
268 isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
269 goto error);
271 better = fold1->type == isl_fold_max ? -1 : 1;
273 if (isl_qpolynomial_fold_is_empty(fold1)) {
274 isl_qpolynomial_fold_free(fold1);
275 return fold2;
278 if (isl_qpolynomial_fold_is_empty(fold2)) {
279 isl_qpolynomial_fold_free(fold2);
280 return fold1;
283 res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
284 fold1->n + fold2->n);
285 if (!res)
286 goto error;
288 for (i = 0; i < fold1->n; ++i) {
289 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
290 if (!res->qp[res->n])
291 goto error;
292 res->n++;
294 n1 = res->n;
296 for (i = 0; i < fold2->n; ++i) {
297 for (j = n1 - 1; j >= 0; --j) {
298 isl_qpolynomial *d;
299 int sgn;
300 d = isl_qpolynomial_sub(
301 isl_qpolynomial_copy(res->qp[j]),
302 isl_qpolynomial_copy(fold2->qp[i]));
303 sgn = isl_qpolynomial_sign(set, d);
304 isl_qpolynomial_free(d);
305 if (sgn == 0)
306 continue;
307 if (sgn != better)
308 break;
309 isl_qpolynomial_free(res->qp[j]);
310 if (j != n1 - 1)
311 res->qp[j] = res->qp[n1 - 1];
312 n1--;
313 if (n1 != res->n - 1)
314 res->qp[n1] = res->qp[res->n - 1];
315 res->n--;
317 if (j >= 0)
318 continue;
319 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
320 if (!res->qp[res->n])
321 goto error;
322 res->n++;
325 isl_qpolynomial_fold_free(fold1);
326 isl_qpolynomial_fold_free(fold2);
328 return res;
329 error:
330 isl_qpolynomial_fold_free(res);
331 isl_qpolynomial_fold_free(fold1);
332 isl_qpolynomial_fold_free(fold2);
333 return NULL;
336 #undef PW
337 #define PW isl_pw_qpolynomial_fold
338 #undef EL
339 #define EL isl_qpolynomial_fold
340 #undef IS_ZERO
341 #define IS_ZERO is_empty
342 #undef FIELD
343 #define FIELD fold
344 #undef ADD
345 #define ADD(d,e1,e2) isl_qpolynomial_fold_fold_on_domain(d,e1,e2)
347 #include <isl_pw_templ.c>
349 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
350 __isl_take isl_dim *dim)
352 return qpolynomial_fold_alloc(type, dim, 0);
355 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
356 enum isl_fold type, __isl_take isl_qpolynomial *qp)
358 isl_qpolynomial_fold *fold;
360 if (!qp)
361 return NULL;
363 fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
364 if (!fold)
365 goto error;
367 fold->qp[0] = qp;
368 fold->n++;
370 return fold;
371 error:
372 isl_qpolynomial_fold_free(fold);
373 isl_qpolynomial_free(qp);
374 return NULL;
377 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
378 __isl_keep isl_qpolynomial_fold *fold)
380 if (!fold)
381 return NULL;
383 fold->ref++;
384 return fold;
387 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
388 __isl_keep isl_qpolynomial_fold *fold)
390 int i;
391 isl_qpolynomial_fold *dup;
393 if (!fold)
394 return NULL;
395 dup = qpolynomial_fold_alloc(fold->type,
396 isl_dim_copy(fold->dim), fold->n);
397 if (!dup)
398 return NULL;
400 for (i = 0; i < fold->n; ++i) {
401 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
402 if (!dup->qp[i])
403 goto error;
406 return dup;
407 error:
408 isl_qpolynomial_fold_free(dup);
409 return NULL;
412 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
413 __isl_take isl_qpolynomial_fold *fold)
415 if (!fold)
416 return NULL;
418 if (fold->ref == 1)
419 return fold;
420 fold->ref--;
421 return isl_qpolynomial_fold_dup(fold);
424 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
426 int i;
428 if (!fold)
429 return;
430 if (--fold->ref > 0)
431 return;
433 for (i = 0; i < fold->n; ++i)
434 isl_qpolynomial_free(fold->qp[i]);
435 isl_dim_free(fold->dim);
436 free(fold);
439 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
441 if (!fold)
442 return -1;
444 return fold->n == 0;
447 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
448 __isl_take isl_qpolynomial_fold *fold1,
449 __isl_take isl_qpolynomial_fold *fold2)
451 int i;
452 struct isl_qpolynomial_fold *res = NULL;
454 if (!fold1 || !fold2)
455 goto error;
457 isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
458 isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
459 goto error);
461 if (isl_qpolynomial_fold_is_empty(fold1)) {
462 isl_qpolynomial_fold_free(fold1);
463 return fold2;
466 if (isl_qpolynomial_fold_is_empty(fold2)) {
467 isl_qpolynomial_fold_free(fold2);
468 return fold1;
471 res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
472 fold1->n + fold2->n);
473 if (!res)
474 goto error;
476 for (i = 0; i < fold1->n; ++i) {
477 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
478 if (!res->qp[res->n])
479 goto error;
480 res->n++;
483 for (i = 0; i < fold2->n; ++i) {
484 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
485 if (!res->qp[res->n])
486 goto error;
487 res->n++;
490 isl_qpolynomial_fold_free(fold1);
491 isl_qpolynomial_fold_free(fold2);
493 return res;
494 error:
495 isl_qpolynomial_fold_free(res);
496 isl_qpolynomial_fold_free(fold1);
497 isl_qpolynomial_fold_free(fold2);
498 return NULL;
501 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
502 enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
504 int i;
505 isl_pw_qpolynomial_fold *pwf;
507 if (!pwqp)
508 return NULL;
510 pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
512 for (i = 0; i < pwqp->n; ++i)
513 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
514 isl_set_copy(pwqp->p[i].set),
515 isl_qpolynomial_fold_alloc(type,
516 isl_qpolynomial_copy(pwqp->p[i].qp)));
518 isl_pw_qpolynomial_free(pwqp);
520 return pwf;
523 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
524 __isl_keep isl_qpolynomial_fold *fold2)
526 int i;
528 if (!fold1 || !fold2)
529 return -1;
531 if (fold1->n != fold2->n)
532 return 0;
534 /* We probably want to sort the qps first... */
535 for (i = 0; i < fold1->n; ++i) {
536 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
537 if (eq < 0 || !eq)
538 return eq;
541 return 1;
544 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
545 __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
547 isl_qpolynomial *qp;
549 if (!fold || !pnt)
550 goto error;
551 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
552 isl_assert(pnt->dim->ctx,
553 fold->type == isl_fold_max || fold->type == isl_fold_min,
554 goto error);
556 if (fold->n == 0)
557 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
558 else {
559 int i;
560 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
561 isl_point_copy(pnt));
562 for (i = 1; i < fold->n; ++i) {
563 isl_qpolynomial *qp_i;
564 qp_i = isl_qpolynomial_eval(
565 isl_qpolynomial_copy(fold->qp[i]),
566 isl_point_copy(pnt));
567 if (fold->type == isl_fold_max)
568 qp = isl_qpolynomial_max_cst(qp, qp_i);
569 else
570 qp = isl_qpolynomial_min_cst(qp, qp_i);
573 isl_qpolynomial_fold_free(fold);
574 isl_point_free(pnt);
576 return qp;
577 error:
578 isl_qpolynomial_fold_free(fold);
579 isl_point_free(pnt);
580 return NULL;
583 size_t isl_pw_qpolynomial_fold_size(__isl_keep isl_pw_qpolynomial_fold *pwf)
585 int i;
586 size_t n = 0;
588 for (i = 0; i < pwf->n; ++i)
589 n += pwf->p[i].fold->n;
591 return n;
594 __isl_give isl_qpolynomial *isl_qpolynomial_fold_opt_on_domain(
595 __isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *set, int max)
597 int i;
598 isl_qpolynomial *opt;
600 if (!set || !fold)
601 goto error;
603 if (fold->n == 0) {
604 isl_dim *dim = isl_dim_copy(fold->dim);
605 isl_set_free(set);
606 isl_qpolynomial_fold_free(fold);
607 return isl_qpolynomial_zero(dim);
610 opt = isl_qpolynomial_opt_on_domain(isl_qpolynomial_copy(fold->qp[0]),
611 isl_set_copy(set), max);
612 for (i = 1; i < fold->n; ++i) {
613 isl_qpolynomial *opt_i;
614 opt_i = isl_qpolynomial_opt_on_domain(
615 isl_qpolynomial_copy(fold->qp[i]),
616 isl_set_copy(set), max);
617 if (max)
618 opt = isl_qpolynomial_max_cst(opt, opt_i);
619 else
620 opt = isl_qpolynomial_min_cst(opt, opt_i);
623 isl_set_free(set);
624 isl_qpolynomial_fold_free(fold);
626 return opt;
627 error:
628 isl_set_free(set);
629 isl_qpolynomial_fold_free(fold);
630 return NULL;