dump isl_ast_expr and isl_ast_node objects in block format
[isl.git] / isl_polynomial.c
blob3772abf6794189e645ff0d256fe0049bda4e5e02
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the MIT 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 <stdlib.h>
12 #define ISL_DIM_H
13 #include <isl_ctx_private.h>
14 #include <isl_map_private.h>
15 #include <isl_factorization.h>
16 #include <isl_lp_private.h>
17 #include <isl_seq.h>
18 #include <isl_union_map_private.h>
19 #include <isl_constraint_private.h>
20 #include <isl_polynomial_private.h>
21 #include <isl_point_private.h>
22 #include <isl_space_private.h>
23 #include <isl_mat_private.h>
24 #include <isl_vec_private.h>
25 #include <isl_range.h>
26 #include <isl_local.h>
27 #include <isl_local_space_private.h>
28 #include <isl_aff_private.h>
29 #include <isl_val_private.h>
30 #include <isl_config.h>
31 #include <isl/deprecated/polynomial_int.h>
33 static unsigned pos(__isl_keep isl_space *dim, enum isl_dim_type type)
35 switch (type) {
36 case isl_dim_param: return 0;
37 case isl_dim_in: return dim->nparam;
38 case isl_dim_out: return dim->nparam + dim->n_in;
39 default: return 0;
43 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
45 if (!up)
46 return -1;
48 return up->var < 0;
51 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
53 if (!up)
54 return NULL;
56 isl_assert(up->ctx, up->var < 0, return NULL);
58 return (struct isl_upoly_cst *)up;
61 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
63 if (!up)
64 return NULL;
66 isl_assert(up->ctx, up->var >= 0, return NULL);
68 return (struct isl_upoly_rec *)up;
71 /* Compare two polynomials.
73 * Return -1 if "up1" is "smaller" than "up2", 1 if "up1" is "greater"
74 * than "up2" and 0 if they are equal.
76 static int isl_upoly_plain_cmp(__isl_keep struct isl_upoly *up1,
77 __isl_keep struct isl_upoly *up2)
79 int i;
80 struct isl_upoly_rec *rec1, *rec2;
82 if (up1 == up2)
83 return 0;
84 if (!up1)
85 return -1;
86 if (!up2)
87 return 1;
88 if (up1->var != up2->var)
89 return up1->var - up2->var;
91 if (isl_upoly_is_cst(up1)) {
92 struct isl_upoly_cst *cst1, *cst2;
93 int cmp;
95 cst1 = isl_upoly_as_cst(up1);
96 cst2 = isl_upoly_as_cst(up2);
97 if (!cst1 || !cst2)
98 return 0;
99 cmp = isl_int_cmp(cst1->n, cst2->n);
100 if (cmp != 0)
101 return cmp;
102 return isl_int_cmp(cst1->d, cst2->d);
105 rec1 = isl_upoly_as_rec(up1);
106 rec2 = isl_upoly_as_rec(up2);
107 if (!rec1 || !rec2)
108 return 0;
110 if (rec1->n != rec2->n)
111 return rec1->n - rec2->n;
113 for (i = 0; i < rec1->n; ++i) {
114 int cmp = isl_upoly_plain_cmp(rec1->p[i], rec2->p[i]);
115 if (cmp != 0)
116 return cmp;
119 return 0;
122 isl_bool isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
123 __isl_keep struct isl_upoly *up2)
125 int i;
126 struct isl_upoly_rec *rec1, *rec2;
128 if (!up1 || !up2)
129 return isl_bool_error;
130 if (up1 == up2)
131 return isl_bool_true;
132 if (up1->var != up2->var)
133 return isl_bool_false;
134 if (isl_upoly_is_cst(up1)) {
135 struct isl_upoly_cst *cst1, *cst2;
136 cst1 = isl_upoly_as_cst(up1);
137 cst2 = isl_upoly_as_cst(up2);
138 if (!cst1 || !cst2)
139 return isl_bool_error;
140 return isl_int_eq(cst1->n, cst2->n) &&
141 isl_int_eq(cst1->d, cst2->d);
144 rec1 = isl_upoly_as_rec(up1);
145 rec2 = isl_upoly_as_rec(up2);
146 if (!rec1 || !rec2)
147 return isl_bool_error;
149 if (rec1->n != rec2->n)
150 return isl_bool_false;
152 for (i = 0; i < rec1->n; ++i) {
153 isl_bool eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
154 if (eq < 0 || !eq)
155 return eq;
158 return isl_bool_true;
161 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
163 struct isl_upoly_cst *cst;
165 if (!up)
166 return -1;
167 if (!isl_upoly_is_cst(up))
168 return 0;
170 cst = isl_upoly_as_cst(up);
171 if (!cst)
172 return -1;
174 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
177 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
179 struct isl_upoly_cst *cst;
181 if (!up)
182 return 0;
183 if (!isl_upoly_is_cst(up))
184 return 0;
186 cst = isl_upoly_as_cst(up);
187 if (!cst)
188 return 0;
190 return isl_int_sgn(cst->n);
193 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
195 struct isl_upoly_cst *cst;
197 if (!up)
198 return -1;
199 if (!isl_upoly_is_cst(up))
200 return 0;
202 cst = isl_upoly_as_cst(up);
203 if (!cst)
204 return -1;
206 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
209 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
211 struct isl_upoly_cst *cst;
213 if (!up)
214 return -1;
215 if (!isl_upoly_is_cst(up))
216 return 0;
218 cst = isl_upoly_as_cst(up);
219 if (!cst)
220 return -1;
222 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
225 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
227 struct isl_upoly_cst *cst;
229 if (!up)
230 return -1;
231 if (!isl_upoly_is_cst(up))
232 return 0;
234 cst = isl_upoly_as_cst(up);
235 if (!cst)
236 return -1;
238 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
241 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
243 struct isl_upoly_cst *cst;
245 if (!up)
246 return -1;
247 if (!isl_upoly_is_cst(up))
248 return 0;
250 cst = isl_upoly_as_cst(up);
251 if (!cst)
252 return -1;
254 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
257 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
259 struct isl_upoly_cst *cst;
261 if (!up)
262 return -1;
263 if (!isl_upoly_is_cst(up))
264 return 0;
266 cst = isl_upoly_as_cst(up);
267 if (!cst)
268 return -1;
270 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
273 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
275 struct isl_upoly_cst *cst;
277 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
278 if (!cst)
279 return NULL;
281 cst->up.ref = 1;
282 cst->up.ctx = ctx;
283 isl_ctx_ref(ctx);
284 cst->up.var = -1;
286 isl_int_init(cst->n);
287 isl_int_init(cst->d);
289 return cst;
292 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
294 struct isl_upoly_cst *cst;
296 cst = isl_upoly_cst_alloc(ctx);
297 if (!cst)
298 return NULL;
300 isl_int_set_si(cst->n, 0);
301 isl_int_set_si(cst->d, 1);
303 return &cst->up;
306 __isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
308 struct isl_upoly_cst *cst;
310 cst = isl_upoly_cst_alloc(ctx);
311 if (!cst)
312 return NULL;
314 isl_int_set_si(cst->n, 1);
315 isl_int_set_si(cst->d, 1);
317 return &cst->up;
320 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
322 struct isl_upoly_cst *cst;
324 cst = isl_upoly_cst_alloc(ctx);
325 if (!cst)
326 return NULL;
328 isl_int_set_si(cst->n, 1);
329 isl_int_set_si(cst->d, 0);
331 return &cst->up;
334 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
336 struct isl_upoly_cst *cst;
338 cst = isl_upoly_cst_alloc(ctx);
339 if (!cst)
340 return NULL;
342 isl_int_set_si(cst->n, -1);
343 isl_int_set_si(cst->d, 0);
345 return &cst->up;
348 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
350 struct isl_upoly_cst *cst;
352 cst = isl_upoly_cst_alloc(ctx);
353 if (!cst)
354 return NULL;
356 isl_int_set_si(cst->n, 0);
357 isl_int_set_si(cst->d, 0);
359 return &cst->up;
362 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
363 isl_int n, isl_int d)
365 struct isl_upoly_cst *cst;
367 cst = isl_upoly_cst_alloc(ctx);
368 if (!cst)
369 return NULL;
371 isl_int_set(cst->n, n);
372 isl_int_set(cst->d, d);
374 return &cst->up;
377 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
378 int var, int size)
380 struct isl_upoly_rec *rec;
382 isl_assert(ctx, var >= 0, return NULL);
383 isl_assert(ctx, size >= 0, return NULL);
384 rec = isl_calloc(ctx, struct isl_upoly_rec,
385 sizeof(struct isl_upoly_rec) +
386 size * sizeof(struct isl_upoly *));
387 if (!rec)
388 return NULL;
390 rec->up.ref = 1;
391 rec->up.ctx = ctx;
392 isl_ctx_ref(ctx);
393 rec->up.var = var;
395 rec->n = 0;
396 rec->size = size;
398 return rec;
401 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
402 __isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
404 qp = isl_qpolynomial_cow(qp);
405 if (!qp || !dim)
406 goto error;
408 isl_space_free(qp->dim);
409 qp->dim = dim;
411 return qp;
412 error:
413 isl_qpolynomial_free(qp);
414 isl_space_free(dim);
415 return NULL;
418 /* Reset the space of "qp". This function is called from isl_pw_templ.c
419 * and doesn't know if the space of an element object is represented
420 * directly or through its domain. It therefore passes along both.
422 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
423 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
424 __isl_take isl_space *domain)
426 isl_space_free(space);
427 return isl_qpolynomial_reset_domain_space(qp, domain);
430 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
432 return qp ? qp->dim->ctx : NULL;
435 __isl_give isl_space *isl_qpolynomial_get_domain_space(
436 __isl_keep isl_qpolynomial *qp)
438 return qp ? isl_space_copy(qp->dim) : NULL;
441 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
443 isl_space *space;
444 if (!qp)
445 return NULL;
446 space = isl_space_copy(qp->dim);
447 space = isl_space_from_domain(space);
448 space = isl_space_add_dims(space, isl_dim_out, 1);
449 return space;
452 /* Externally, an isl_qpolynomial has a map space, but internally, the
453 * ls field corresponds to the domain of that space.
455 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
456 enum isl_dim_type type)
458 if (!qp)
459 return 0;
460 if (type == isl_dim_out)
461 return 1;
462 if (type == isl_dim_in)
463 type = isl_dim_set;
464 return isl_space_dim(qp->dim, type);
467 isl_bool isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
469 return qp ? isl_upoly_is_zero(qp->upoly) : isl_bool_error;
472 isl_bool isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
474 return qp ? isl_upoly_is_one(qp->upoly) : isl_bool_error;
477 isl_bool isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
479 return qp ? isl_upoly_is_nan(qp->upoly) : isl_bool_error;
482 isl_bool isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
484 return qp ? isl_upoly_is_infty(qp->upoly) : isl_bool_error;
487 isl_bool isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
489 return qp ? isl_upoly_is_neginfty(qp->upoly) : isl_bool_error;
492 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
494 return qp ? isl_upoly_sgn(qp->upoly) : 0;
497 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
499 isl_int_clear(cst->n);
500 isl_int_clear(cst->d);
503 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
505 int i;
507 for (i = 0; i < rec->n; ++i)
508 isl_upoly_free(rec->p[i]);
511 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
513 if (!up)
514 return NULL;
516 up->ref++;
517 return up;
520 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
522 struct isl_upoly_cst *cst;
523 struct isl_upoly_cst *dup;
525 cst = isl_upoly_as_cst(up);
526 if (!cst)
527 return NULL;
529 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
530 if (!dup)
531 return NULL;
532 isl_int_set(dup->n, cst->n);
533 isl_int_set(dup->d, cst->d);
535 return &dup->up;
538 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
540 int i;
541 struct isl_upoly_rec *rec;
542 struct isl_upoly_rec *dup;
544 rec = isl_upoly_as_rec(up);
545 if (!rec)
546 return NULL;
548 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
549 if (!dup)
550 return NULL;
552 for (i = 0; i < rec->n; ++i) {
553 dup->p[i] = isl_upoly_copy(rec->p[i]);
554 if (!dup->p[i])
555 goto error;
556 dup->n++;
559 return &dup->up;
560 error:
561 isl_upoly_free(&dup->up);
562 return NULL;
565 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
567 if (!up)
568 return NULL;
570 if (isl_upoly_is_cst(up))
571 return isl_upoly_dup_cst(up);
572 else
573 return isl_upoly_dup_rec(up);
576 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
578 if (!up)
579 return NULL;
581 if (up->ref == 1)
582 return up;
583 up->ref--;
584 return isl_upoly_dup(up);
587 void isl_upoly_free(__isl_take struct isl_upoly *up)
589 if (!up)
590 return;
592 if (--up->ref > 0)
593 return;
595 if (up->var < 0)
596 upoly_free_cst((struct isl_upoly_cst *)up);
597 else
598 upoly_free_rec((struct isl_upoly_rec *)up);
600 isl_ctx_deref(up->ctx);
601 free(up);
604 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
606 isl_int gcd;
608 isl_int_init(gcd);
609 isl_int_gcd(gcd, cst->n, cst->d);
610 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
611 isl_int_divexact(cst->n, cst->n, gcd);
612 isl_int_divexact(cst->d, cst->d, gcd);
614 isl_int_clear(gcd);
617 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
618 __isl_take struct isl_upoly *up2)
620 struct isl_upoly_cst *cst1;
621 struct isl_upoly_cst *cst2;
623 up1 = isl_upoly_cow(up1);
624 if (!up1 || !up2)
625 goto error;
627 cst1 = isl_upoly_as_cst(up1);
628 cst2 = isl_upoly_as_cst(up2);
630 if (isl_int_eq(cst1->d, cst2->d))
631 isl_int_add(cst1->n, cst1->n, cst2->n);
632 else {
633 isl_int_mul(cst1->n, cst1->n, cst2->d);
634 isl_int_addmul(cst1->n, cst2->n, cst1->d);
635 isl_int_mul(cst1->d, cst1->d, cst2->d);
638 isl_upoly_cst_reduce(cst1);
640 isl_upoly_free(up2);
641 return up1;
642 error:
643 isl_upoly_free(up1);
644 isl_upoly_free(up2);
645 return NULL;
648 static __isl_give struct isl_upoly *replace_by_zero(
649 __isl_take struct isl_upoly *up)
651 struct isl_ctx *ctx;
653 if (!up)
654 return NULL;
655 ctx = up->ctx;
656 isl_upoly_free(up);
657 return isl_upoly_zero(ctx);
660 static __isl_give struct isl_upoly *replace_by_constant_term(
661 __isl_take struct isl_upoly *up)
663 struct isl_upoly_rec *rec;
664 struct isl_upoly *cst;
666 if (!up)
667 return NULL;
669 rec = isl_upoly_as_rec(up);
670 if (!rec)
671 goto error;
672 cst = isl_upoly_copy(rec->p[0]);
673 isl_upoly_free(up);
674 return cst;
675 error:
676 isl_upoly_free(up);
677 return NULL;
680 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
681 __isl_take struct isl_upoly *up2)
683 int i;
684 struct isl_upoly_rec *rec1, *rec2;
686 if (!up1 || !up2)
687 goto error;
689 if (isl_upoly_is_nan(up1)) {
690 isl_upoly_free(up2);
691 return up1;
694 if (isl_upoly_is_nan(up2)) {
695 isl_upoly_free(up1);
696 return up2;
699 if (isl_upoly_is_zero(up1)) {
700 isl_upoly_free(up1);
701 return up2;
704 if (isl_upoly_is_zero(up2)) {
705 isl_upoly_free(up2);
706 return up1;
709 if (up1->var < up2->var)
710 return isl_upoly_sum(up2, up1);
712 if (up2->var < up1->var) {
713 struct isl_upoly_rec *rec;
714 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
715 isl_upoly_free(up1);
716 return up2;
718 up1 = isl_upoly_cow(up1);
719 rec = isl_upoly_as_rec(up1);
720 if (!rec)
721 goto error;
722 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
723 if (rec->n == 1)
724 up1 = replace_by_constant_term(up1);
725 return up1;
728 if (isl_upoly_is_cst(up1))
729 return isl_upoly_sum_cst(up1, up2);
731 rec1 = isl_upoly_as_rec(up1);
732 rec2 = isl_upoly_as_rec(up2);
733 if (!rec1 || !rec2)
734 goto error;
736 if (rec1->n < rec2->n)
737 return isl_upoly_sum(up2, up1);
739 up1 = isl_upoly_cow(up1);
740 rec1 = isl_upoly_as_rec(up1);
741 if (!rec1)
742 goto error;
744 for (i = rec2->n - 1; i >= 0; --i) {
745 rec1->p[i] = isl_upoly_sum(rec1->p[i],
746 isl_upoly_copy(rec2->p[i]));
747 if (!rec1->p[i])
748 goto error;
749 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
750 isl_upoly_free(rec1->p[i]);
751 rec1->n--;
755 if (rec1->n == 0)
756 up1 = replace_by_zero(up1);
757 else if (rec1->n == 1)
758 up1 = replace_by_constant_term(up1);
760 isl_upoly_free(up2);
762 return up1;
763 error:
764 isl_upoly_free(up1);
765 isl_upoly_free(up2);
766 return NULL;
769 __isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
770 __isl_take struct isl_upoly *up, isl_int v)
772 struct isl_upoly_cst *cst;
774 up = isl_upoly_cow(up);
775 if (!up)
776 return NULL;
778 cst = isl_upoly_as_cst(up);
780 isl_int_addmul(cst->n, cst->d, v);
782 return up;
785 __isl_give struct isl_upoly *isl_upoly_add_isl_int(
786 __isl_take struct isl_upoly *up, isl_int v)
788 struct isl_upoly_rec *rec;
790 if (!up)
791 return NULL;
793 if (isl_upoly_is_cst(up))
794 return isl_upoly_cst_add_isl_int(up, v);
796 up = isl_upoly_cow(up);
797 rec = isl_upoly_as_rec(up);
798 if (!rec)
799 goto error;
801 rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
802 if (!rec->p[0])
803 goto error;
805 return up;
806 error:
807 isl_upoly_free(up);
808 return NULL;
811 __isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
812 __isl_take struct isl_upoly *up, isl_int v)
814 struct isl_upoly_cst *cst;
816 if (isl_upoly_is_zero(up))
817 return up;
819 up = isl_upoly_cow(up);
820 if (!up)
821 return NULL;
823 cst = isl_upoly_as_cst(up);
825 isl_int_mul(cst->n, cst->n, v);
827 return up;
830 __isl_give struct isl_upoly *isl_upoly_mul_isl_int(
831 __isl_take struct isl_upoly *up, isl_int v)
833 int i;
834 struct isl_upoly_rec *rec;
836 if (!up)
837 return NULL;
839 if (isl_upoly_is_cst(up))
840 return isl_upoly_cst_mul_isl_int(up, v);
842 up = isl_upoly_cow(up);
843 rec = isl_upoly_as_rec(up);
844 if (!rec)
845 goto error;
847 for (i = 0; i < rec->n; ++i) {
848 rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
849 if (!rec->p[i])
850 goto error;
853 return up;
854 error:
855 isl_upoly_free(up);
856 return NULL;
859 /* Multiply the constant polynomial "up" by "v".
861 static __isl_give struct isl_upoly *isl_upoly_cst_scale_val(
862 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
864 struct isl_upoly_cst *cst;
866 if (isl_upoly_is_zero(up))
867 return up;
869 up = isl_upoly_cow(up);
870 if (!up)
871 return NULL;
873 cst = isl_upoly_as_cst(up);
875 isl_int_mul(cst->n, cst->n, v->n);
876 isl_int_mul(cst->d, cst->d, v->d);
877 isl_upoly_cst_reduce(cst);
879 return up;
882 /* Multiply the polynomial "up" by "v".
884 static __isl_give struct isl_upoly *isl_upoly_scale_val(
885 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
887 int i;
888 struct isl_upoly_rec *rec;
890 if (!up)
891 return NULL;
893 if (isl_upoly_is_cst(up))
894 return isl_upoly_cst_scale_val(up, v);
896 up = isl_upoly_cow(up);
897 rec = isl_upoly_as_rec(up);
898 if (!rec)
899 goto error;
901 for (i = 0; i < rec->n; ++i) {
902 rec->p[i] = isl_upoly_scale_val(rec->p[i], v);
903 if (!rec->p[i])
904 goto error;
907 return up;
908 error:
909 isl_upoly_free(up);
910 return NULL;
913 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
914 __isl_take struct isl_upoly *up2)
916 struct isl_upoly_cst *cst1;
917 struct isl_upoly_cst *cst2;
919 up1 = isl_upoly_cow(up1);
920 if (!up1 || !up2)
921 goto error;
923 cst1 = isl_upoly_as_cst(up1);
924 cst2 = isl_upoly_as_cst(up2);
926 isl_int_mul(cst1->n, cst1->n, cst2->n);
927 isl_int_mul(cst1->d, cst1->d, cst2->d);
929 isl_upoly_cst_reduce(cst1);
931 isl_upoly_free(up2);
932 return up1;
933 error:
934 isl_upoly_free(up1);
935 isl_upoly_free(up2);
936 return NULL;
939 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
940 __isl_take struct isl_upoly *up2)
942 struct isl_upoly_rec *rec1;
943 struct isl_upoly_rec *rec2;
944 struct isl_upoly_rec *res = NULL;
945 int i, j;
946 int size;
948 rec1 = isl_upoly_as_rec(up1);
949 rec2 = isl_upoly_as_rec(up2);
950 if (!rec1 || !rec2)
951 goto error;
952 size = rec1->n + rec2->n - 1;
953 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
954 if (!res)
955 goto error;
957 for (i = 0; i < rec1->n; ++i) {
958 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
959 isl_upoly_copy(rec1->p[i]));
960 if (!res->p[i])
961 goto error;
962 res->n++;
964 for (; i < size; ++i) {
965 res->p[i] = isl_upoly_zero(up1->ctx);
966 if (!res->p[i])
967 goto error;
968 res->n++;
970 for (i = 0; i < rec1->n; ++i) {
971 for (j = 1; j < rec2->n; ++j) {
972 struct isl_upoly *up;
973 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
974 isl_upoly_copy(rec1->p[i]));
975 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
976 if (!res->p[i + j])
977 goto error;
981 isl_upoly_free(up1);
982 isl_upoly_free(up2);
984 return &res->up;
985 error:
986 isl_upoly_free(up1);
987 isl_upoly_free(up2);
988 isl_upoly_free(&res->up);
989 return NULL;
992 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
993 __isl_take struct isl_upoly *up2)
995 if (!up1 || !up2)
996 goto error;
998 if (isl_upoly_is_nan(up1)) {
999 isl_upoly_free(up2);
1000 return up1;
1003 if (isl_upoly_is_nan(up2)) {
1004 isl_upoly_free(up1);
1005 return up2;
1008 if (isl_upoly_is_zero(up1)) {
1009 isl_upoly_free(up2);
1010 return up1;
1013 if (isl_upoly_is_zero(up2)) {
1014 isl_upoly_free(up1);
1015 return up2;
1018 if (isl_upoly_is_one(up1)) {
1019 isl_upoly_free(up1);
1020 return up2;
1023 if (isl_upoly_is_one(up2)) {
1024 isl_upoly_free(up2);
1025 return up1;
1028 if (up1->var < up2->var)
1029 return isl_upoly_mul(up2, up1);
1031 if (up2->var < up1->var) {
1032 int i;
1033 struct isl_upoly_rec *rec;
1034 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
1035 isl_ctx *ctx = up1->ctx;
1036 isl_upoly_free(up1);
1037 isl_upoly_free(up2);
1038 return isl_upoly_nan(ctx);
1040 up1 = isl_upoly_cow(up1);
1041 rec = isl_upoly_as_rec(up1);
1042 if (!rec)
1043 goto error;
1045 for (i = 0; i < rec->n; ++i) {
1046 rec->p[i] = isl_upoly_mul(rec->p[i],
1047 isl_upoly_copy(up2));
1048 if (!rec->p[i])
1049 goto error;
1051 isl_upoly_free(up2);
1052 return up1;
1055 if (isl_upoly_is_cst(up1))
1056 return isl_upoly_mul_cst(up1, up2);
1058 return isl_upoly_mul_rec(up1, up2);
1059 error:
1060 isl_upoly_free(up1);
1061 isl_upoly_free(up2);
1062 return NULL;
1065 __isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
1066 unsigned power)
1068 struct isl_upoly *res;
1070 if (!up)
1071 return NULL;
1072 if (power == 1)
1073 return up;
1075 if (power % 2)
1076 res = isl_upoly_copy(up);
1077 else
1078 res = isl_upoly_one(up->ctx);
1080 while (power >>= 1) {
1081 up = isl_upoly_mul(up, isl_upoly_copy(up));
1082 if (power % 2)
1083 res = isl_upoly_mul(res, isl_upoly_copy(up));
1086 isl_upoly_free(up);
1087 return res;
1090 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
1091 unsigned n_div, __isl_take struct isl_upoly *up)
1093 struct isl_qpolynomial *qp = NULL;
1094 unsigned total;
1096 if (!dim || !up)
1097 goto error;
1099 if (!isl_space_is_set(dim))
1100 isl_die(isl_space_get_ctx(dim), isl_error_invalid,
1101 "domain of polynomial should be a set", goto error);
1103 total = isl_space_dim(dim, isl_dim_all);
1105 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
1106 if (!qp)
1107 goto error;
1109 qp->ref = 1;
1110 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
1111 if (!qp->div)
1112 goto error;
1114 qp->dim = dim;
1115 qp->upoly = up;
1117 return qp;
1118 error:
1119 isl_space_free(dim);
1120 isl_upoly_free(up);
1121 isl_qpolynomial_free(qp);
1122 return NULL;
1125 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1127 if (!qp)
1128 return NULL;
1130 qp->ref++;
1131 return qp;
1134 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1136 struct isl_qpolynomial *dup;
1138 if (!qp)
1139 return NULL;
1141 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1142 isl_upoly_copy(qp->upoly));
1143 if (!dup)
1144 return NULL;
1145 isl_mat_free(dup->div);
1146 dup->div = isl_mat_copy(qp->div);
1147 if (!dup->div)
1148 goto error;
1150 return dup;
1151 error:
1152 isl_qpolynomial_free(dup);
1153 return NULL;
1156 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1158 if (!qp)
1159 return NULL;
1161 if (qp->ref == 1)
1162 return qp;
1163 qp->ref--;
1164 return isl_qpolynomial_dup(qp);
1167 __isl_null isl_qpolynomial *isl_qpolynomial_free(
1168 __isl_take isl_qpolynomial *qp)
1170 if (!qp)
1171 return NULL;
1173 if (--qp->ref > 0)
1174 return NULL;
1176 isl_space_free(qp->dim);
1177 isl_mat_free(qp->div);
1178 isl_upoly_free(qp->upoly);
1180 free(qp);
1181 return NULL;
1184 __isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
1186 int i;
1187 struct isl_upoly_rec *rec;
1188 struct isl_upoly_cst *cst;
1190 rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1191 if (!rec)
1192 return NULL;
1193 for (i = 0; i < 1 + power; ++i) {
1194 rec->p[i] = isl_upoly_zero(ctx);
1195 if (!rec->p[i])
1196 goto error;
1197 rec->n++;
1199 cst = isl_upoly_as_cst(rec->p[power]);
1200 isl_int_set_si(cst->n, 1);
1202 return &rec->up;
1203 error:
1204 isl_upoly_free(&rec->up);
1205 return NULL;
1208 /* r array maps original positions to new positions.
1210 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
1211 int *r)
1213 int i;
1214 struct isl_upoly_rec *rec;
1215 struct isl_upoly *base;
1216 struct isl_upoly *res;
1218 if (isl_upoly_is_cst(up))
1219 return up;
1221 rec = isl_upoly_as_rec(up);
1222 if (!rec)
1223 goto error;
1225 isl_assert(up->ctx, rec->n >= 1, goto error);
1227 base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
1228 res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
1230 for (i = rec->n - 2; i >= 0; --i) {
1231 res = isl_upoly_mul(res, isl_upoly_copy(base));
1232 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
1235 isl_upoly_free(base);
1236 isl_upoly_free(up);
1238 return res;
1239 error:
1240 isl_upoly_free(up);
1241 return NULL;
1244 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
1246 int n_row, n_col;
1247 int equal;
1249 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1250 div1->n_col >= div2->n_col, return -1);
1252 if (div1->n_row == div2->n_row)
1253 return isl_mat_is_equal(div1, div2);
1255 n_row = div1->n_row;
1256 n_col = div1->n_col;
1257 div1->n_row = div2->n_row;
1258 div1->n_col = div2->n_col;
1260 equal = isl_mat_is_equal(div1, div2);
1262 div1->n_row = n_row;
1263 div1->n_col = n_col;
1265 return equal;
1268 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1270 int li, lj;
1272 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1273 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1275 if (li != lj)
1276 return li - lj;
1278 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1281 struct isl_div_sort_info {
1282 isl_mat *div;
1283 int row;
1286 static int div_sort_cmp(const void *p1, const void *p2)
1288 const struct isl_div_sort_info *i1, *i2;
1289 i1 = (const struct isl_div_sort_info *) p1;
1290 i2 = (const struct isl_div_sort_info *) p2;
1292 return cmp_row(i1->div, i1->row, i2->row);
1295 /* Sort divs and remove duplicates.
1297 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1299 int i;
1300 int skip;
1301 int len;
1302 struct isl_div_sort_info *array = NULL;
1303 int *pos = NULL, *at = NULL;
1304 int *reordering = NULL;
1305 unsigned div_pos;
1307 if (!qp)
1308 return NULL;
1309 if (qp->div->n_row <= 1)
1310 return qp;
1312 div_pos = isl_space_dim(qp->dim, isl_dim_all);
1314 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1315 qp->div->n_row);
1316 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1317 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1318 len = qp->div->n_col - 2;
1319 reordering = isl_alloc_array(qp->div->ctx, int, len);
1320 if (!array || !pos || !at || !reordering)
1321 goto error;
1323 for (i = 0; i < qp->div->n_row; ++i) {
1324 array[i].div = qp->div;
1325 array[i].row = i;
1326 pos[i] = i;
1327 at[i] = i;
1330 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1331 div_sort_cmp);
1333 for (i = 0; i < div_pos; ++i)
1334 reordering[i] = i;
1336 for (i = 0; i < qp->div->n_row; ++i) {
1337 if (pos[array[i].row] == i)
1338 continue;
1339 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1340 pos[at[i]] = pos[array[i].row];
1341 at[pos[array[i].row]] = at[i];
1342 at[i] = array[i].row;
1343 pos[array[i].row] = i;
1346 skip = 0;
1347 for (i = 0; i < len - div_pos; ++i) {
1348 if (i > 0 &&
1349 isl_seq_eq(qp->div->row[i - skip - 1],
1350 qp->div->row[i - skip], qp->div->n_col)) {
1351 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1352 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1353 2 + div_pos + i - skip);
1354 qp->div = isl_mat_drop_cols(qp->div,
1355 2 + div_pos + i - skip, 1);
1356 skip++;
1358 reordering[div_pos + array[i].row] = div_pos + i - skip;
1361 qp->upoly = reorder(qp->upoly, reordering);
1363 if (!qp->upoly || !qp->div)
1364 goto error;
1366 free(at);
1367 free(pos);
1368 free(array);
1369 free(reordering);
1371 return qp;
1372 error:
1373 free(at);
1374 free(pos);
1375 free(array);
1376 free(reordering);
1377 isl_qpolynomial_free(qp);
1378 return NULL;
1381 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1382 int *exp, int first)
1384 int i;
1385 struct isl_upoly_rec *rec;
1387 if (isl_upoly_is_cst(up))
1388 return up;
1390 if (up->var < first)
1391 return up;
1393 if (exp[up->var - first] == up->var - first)
1394 return up;
1396 up = isl_upoly_cow(up);
1397 if (!up)
1398 goto error;
1400 up->var = exp[up->var - first] + first;
1402 rec = isl_upoly_as_rec(up);
1403 if (!rec)
1404 goto error;
1406 for (i = 0; i < rec->n; ++i) {
1407 rec->p[i] = expand(rec->p[i], exp, first);
1408 if (!rec->p[i])
1409 goto error;
1412 return up;
1413 error:
1414 isl_upoly_free(up);
1415 return NULL;
1418 static __isl_give isl_qpolynomial *with_merged_divs(
1419 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1420 __isl_take isl_qpolynomial *qp2),
1421 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1423 int *exp1 = NULL;
1424 int *exp2 = NULL;
1425 isl_mat *div = NULL;
1426 int n_div1, n_div2;
1428 qp1 = isl_qpolynomial_cow(qp1);
1429 qp2 = isl_qpolynomial_cow(qp2);
1431 if (!qp1 || !qp2)
1432 goto error;
1434 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1435 qp1->div->n_col >= qp2->div->n_col, goto error);
1437 n_div1 = qp1->div->n_row;
1438 n_div2 = qp2->div->n_row;
1439 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1);
1440 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2);
1441 if ((n_div1 && !exp1) || (n_div2 && !exp2))
1442 goto error;
1444 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1445 if (!div)
1446 goto error;
1448 isl_mat_free(qp1->div);
1449 qp1->div = isl_mat_copy(div);
1450 isl_mat_free(qp2->div);
1451 qp2->div = isl_mat_copy(div);
1453 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1454 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1456 if (!qp1->upoly || !qp2->upoly)
1457 goto error;
1459 isl_mat_free(div);
1460 free(exp1);
1461 free(exp2);
1463 return fn(qp1, qp2);
1464 error:
1465 isl_mat_free(div);
1466 free(exp1);
1467 free(exp2);
1468 isl_qpolynomial_free(qp1);
1469 isl_qpolynomial_free(qp2);
1470 return NULL;
1473 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1474 __isl_take isl_qpolynomial *qp2)
1476 qp1 = isl_qpolynomial_cow(qp1);
1478 if (!qp1 || !qp2)
1479 goto error;
1481 if (qp1->div->n_row < qp2->div->n_row)
1482 return isl_qpolynomial_add(qp2, qp1);
1484 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1485 if (!compatible_divs(qp1->div, qp2->div))
1486 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1488 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1489 if (!qp1->upoly)
1490 goto error;
1492 isl_qpolynomial_free(qp2);
1494 return qp1;
1495 error:
1496 isl_qpolynomial_free(qp1);
1497 isl_qpolynomial_free(qp2);
1498 return NULL;
1501 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1502 __isl_keep isl_set *dom,
1503 __isl_take isl_qpolynomial *qp1,
1504 __isl_take isl_qpolynomial *qp2)
1506 qp1 = isl_qpolynomial_add(qp1, qp2);
1507 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1508 return qp1;
1511 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1512 __isl_take isl_qpolynomial *qp2)
1514 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1517 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1518 __isl_take isl_qpolynomial *qp, isl_int v)
1520 if (isl_int_is_zero(v))
1521 return qp;
1523 qp = isl_qpolynomial_cow(qp);
1524 if (!qp)
1525 return NULL;
1527 qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1528 if (!qp->upoly)
1529 goto error;
1531 return qp;
1532 error:
1533 isl_qpolynomial_free(qp);
1534 return NULL;
1538 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1540 if (!qp)
1541 return NULL;
1543 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1546 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1547 __isl_take isl_qpolynomial *qp, isl_int v)
1549 if (isl_int_is_one(v))
1550 return qp;
1552 if (qp && isl_int_is_zero(v)) {
1553 isl_qpolynomial *zero;
1554 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1555 isl_qpolynomial_free(qp);
1556 return zero;
1559 qp = isl_qpolynomial_cow(qp);
1560 if (!qp)
1561 return NULL;
1563 qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1564 if (!qp->upoly)
1565 goto error;
1567 return qp;
1568 error:
1569 isl_qpolynomial_free(qp);
1570 return NULL;
1573 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1574 __isl_take isl_qpolynomial *qp, isl_int v)
1576 return isl_qpolynomial_mul_isl_int(qp, v);
1579 /* Multiply "qp" by "v".
1581 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1582 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1584 if (!qp || !v)
1585 goto error;
1587 if (!isl_val_is_rat(v))
1588 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1589 "expecting rational factor", goto error);
1591 if (isl_val_is_one(v)) {
1592 isl_val_free(v);
1593 return qp;
1596 if (isl_val_is_zero(v)) {
1597 isl_space *space;
1599 space = isl_qpolynomial_get_domain_space(qp);
1600 isl_qpolynomial_free(qp);
1601 isl_val_free(v);
1602 return isl_qpolynomial_zero_on_domain(space);
1605 qp = isl_qpolynomial_cow(qp);
1606 if (!qp)
1607 goto error;
1609 qp->upoly = isl_upoly_scale_val(qp->upoly, v);
1610 if (!qp->upoly)
1611 qp = isl_qpolynomial_free(qp);
1613 isl_val_free(v);
1614 return qp;
1615 error:
1616 isl_val_free(v);
1617 isl_qpolynomial_free(qp);
1618 return NULL;
1621 /* Divide "qp" by "v".
1623 __isl_give isl_qpolynomial *isl_qpolynomial_scale_down_val(
1624 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1626 if (!qp || !v)
1627 goto error;
1629 if (!isl_val_is_rat(v))
1630 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1631 "expecting rational factor", goto error);
1632 if (isl_val_is_zero(v))
1633 isl_die(isl_val_get_ctx(v), isl_error_invalid,
1634 "cannot scale down by zero", goto error);
1636 return isl_qpolynomial_scale_val(qp, isl_val_inv(v));
1637 error:
1638 isl_val_free(v);
1639 isl_qpolynomial_free(qp);
1640 return NULL;
1643 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1644 __isl_take isl_qpolynomial *qp2)
1646 qp1 = isl_qpolynomial_cow(qp1);
1648 if (!qp1 || !qp2)
1649 goto error;
1651 if (qp1->div->n_row < qp2->div->n_row)
1652 return isl_qpolynomial_mul(qp2, qp1);
1654 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1655 if (!compatible_divs(qp1->div, qp2->div))
1656 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1658 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1659 if (!qp1->upoly)
1660 goto error;
1662 isl_qpolynomial_free(qp2);
1664 return qp1;
1665 error:
1666 isl_qpolynomial_free(qp1);
1667 isl_qpolynomial_free(qp2);
1668 return NULL;
1671 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1672 unsigned power)
1674 qp = isl_qpolynomial_cow(qp);
1676 if (!qp)
1677 return NULL;
1679 qp->upoly = isl_upoly_pow(qp->upoly, power);
1680 if (!qp->upoly)
1681 goto error;
1683 return qp;
1684 error:
1685 isl_qpolynomial_free(qp);
1686 return NULL;
1689 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1690 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1692 int i;
1694 if (power == 1)
1695 return pwqp;
1697 pwqp = isl_pw_qpolynomial_cow(pwqp);
1698 if (!pwqp)
1699 return NULL;
1701 for (i = 0; i < pwqp->n; ++i) {
1702 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1703 if (!pwqp->p[i].qp)
1704 return isl_pw_qpolynomial_free(pwqp);
1707 return pwqp;
1710 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1711 __isl_take isl_space *dim)
1713 if (!dim)
1714 return NULL;
1715 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1718 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1719 __isl_take isl_space *dim)
1721 if (!dim)
1722 return NULL;
1723 return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1726 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1727 __isl_take isl_space *dim)
1729 if (!dim)
1730 return NULL;
1731 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1734 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1735 __isl_take isl_space *dim)
1737 if (!dim)
1738 return NULL;
1739 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1742 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1743 __isl_take isl_space *dim)
1745 if (!dim)
1746 return NULL;
1747 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1750 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1751 __isl_take isl_space *dim,
1752 isl_int v)
1754 struct isl_qpolynomial *qp;
1755 struct isl_upoly_cst *cst;
1757 if (!dim)
1758 return NULL;
1760 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1761 if (!qp)
1762 return NULL;
1764 cst = isl_upoly_as_cst(qp->upoly);
1765 isl_int_set(cst->n, v);
1767 return qp;
1770 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1771 isl_int *n, isl_int *d)
1773 struct isl_upoly_cst *cst;
1775 if (!qp)
1776 return -1;
1778 if (!isl_upoly_is_cst(qp->upoly))
1779 return 0;
1781 cst = isl_upoly_as_cst(qp->upoly);
1782 if (!cst)
1783 return -1;
1785 if (n)
1786 isl_int_set(*n, cst->n);
1787 if (d)
1788 isl_int_set(*d, cst->d);
1790 return 1;
1793 /* Return the constant term of "up".
1795 static __isl_give isl_val *isl_upoly_get_constant_val(
1796 __isl_keep struct isl_upoly *up)
1798 struct isl_upoly_cst *cst;
1800 if (!up)
1801 return NULL;
1803 while (!isl_upoly_is_cst(up)) {
1804 struct isl_upoly_rec *rec;
1806 rec = isl_upoly_as_rec(up);
1807 if (!rec)
1808 return NULL;
1809 up = rec->p[0];
1812 cst = isl_upoly_as_cst(up);
1813 if (!cst)
1814 return NULL;
1815 return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1818 /* Return the constant term of "qp".
1820 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1821 __isl_keep isl_qpolynomial *qp)
1823 if (!qp)
1824 return NULL;
1826 return isl_upoly_get_constant_val(qp->upoly);
1829 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1831 int is_cst;
1832 struct isl_upoly_rec *rec;
1834 if (!up)
1835 return -1;
1837 if (up->var < 0)
1838 return 1;
1840 rec = isl_upoly_as_rec(up);
1841 if (!rec)
1842 return -1;
1844 if (rec->n > 2)
1845 return 0;
1847 isl_assert(up->ctx, rec->n > 1, return -1);
1849 is_cst = isl_upoly_is_cst(rec->p[1]);
1850 if (is_cst < 0)
1851 return -1;
1852 if (!is_cst)
1853 return 0;
1855 return isl_upoly_is_affine(rec->p[0]);
1858 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1860 if (!qp)
1861 return -1;
1863 if (qp->div->n_row > 0)
1864 return 0;
1866 return isl_upoly_is_affine(qp->upoly);
1869 static void update_coeff(__isl_keep isl_vec *aff,
1870 __isl_keep struct isl_upoly_cst *cst, int pos)
1872 isl_int gcd;
1873 isl_int f;
1875 if (isl_int_is_zero(cst->n))
1876 return;
1878 isl_int_init(gcd);
1879 isl_int_init(f);
1880 isl_int_gcd(gcd, cst->d, aff->el[0]);
1881 isl_int_divexact(f, cst->d, gcd);
1882 isl_int_divexact(gcd, aff->el[0], gcd);
1883 isl_seq_scale(aff->el, aff->el, f, aff->size);
1884 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1885 isl_int_clear(gcd);
1886 isl_int_clear(f);
1889 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1890 __isl_keep isl_vec *aff)
1892 struct isl_upoly_cst *cst;
1893 struct isl_upoly_rec *rec;
1895 if (!up || !aff)
1896 return -1;
1898 if (up->var < 0) {
1899 struct isl_upoly_cst *cst;
1901 cst = isl_upoly_as_cst(up);
1902 if (!cst)
1903 return -1;
1904 update_coeff(aff, cst, 0);
1905 return 0;
1908 rec = isl_upoly_as_rec(up);
1909 if (!rec)
1910 return -1;
1911 isl_assert(up->ctx, rec->n == 2, return -1);
1913 cst = isl_upoly_as_cst(rec->p[1]);
1914 if (!cst)
1915 return -1;
1916 update_coeff(aff, cst, 1 + up->var);
1918 return isl_upoly_update_affine(rec->p[0], aff);
1921 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1922 __isl_keep isl_qpolynomial *qp)
1924 isl_vec *aff;
1925 unsigned d;
1927 if (!qp)
1928 return NULL;
1930 d = isl_space_dim(qp->dim, isl_dim_all);
1931 aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1932 if (!aff)
1933 return NULL;
1935 isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1936 isl_int_set_si(aff->el[0], 1);
1938 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1939 goto error;
1941 return aff;
1942 error:
1943 isl_vec_free(aff);
1944 return NULL;
1947 /* Compare two quasi-polynomials.
1949 * Return -1 if "qp1" is "smaller" than "qp2", 1 if "qp1" is "greater"
1950 * than "qp2" and 0 if they are equal.
1952 int isl_qpolynomial_plain_cmp(__isl_keep isl_qpolynomial *qp1,
1953 __isl_keep isl_qpolynomial *qp2)
1955 int cmp;
1957 if (qp1 == qp2)
1958 return 0;
1959 if (!qp1)
1960 return -1;
1961 if (!qp2)
1962 return 1;
1964 cmp = isl_space_cmp(qp1->dim, qp2->dim);
1965 if (cmp != 0)
1966 return cmp;
1968 cmp = isl_local_cmp(qp1->div, qp2->div);
1969 if (cmp != 0)
1970 return cmp;
1972 return isl_upoly_plain_cmp(qp1->upoly, qp2->upoly);
1975 /* Is "qp1" obviously equal to "qp2"?
1977 * NaN is not equal to anything, not even to another NaN.
1979 isl_bool isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
1980 __isl_keep isl_qpolynomial *qp2)
1982 isl_bool equal;
1984 if (!qp1 || !qp2)
1985 return isl_bool_error;
1987 if (isl_qpolynomial_is_nan(qp1) || isl_qpolynomial_is_nan(qp2))
1988 return isl_bool_false;
1990 equal = isl_space_is_equal(qp1->dim, qp2->dim);
1991 if (equal < 0 || !equal)
1992 return equal;
1994 equal = isl_mat_is_equal(qp1->div, qp2->div);
1995 if (equal < 0 || !equal)
1996 return equal;
1998 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
2001 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
2003 int i;
2004 struct isl_upoly_rec *rec;
2006 if (isl_upoly_is_cst(up)) {
2007 struct isl_upoly_cst *cst;
2008 cst = isl_upoly_as_cst(up);
2009 if (!cst)
2010 return;
2011 isl_int_lcm(*d, *d, cst->d);
2012 return;
2015 rec = isl_upoly_as_rec(up);
2016 if (!rec)
2017 return;
2019 for (i = 0; i < rec->n; ++i)
2020 upoly_update_den(rec->p[i], d);
2023 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
2025 isl_int_set_si(*d, 1);
2026 if (!qp)
2027 return;
2028 upoly_update_den(qp->upoly, d);
2031 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
2032 __isl_take isl_space *dim, int pos, int power)
2034 struct isl_ctx *ctx;
2036 if (!dim)
2037 return NULL;
2039 ctx = dim->ctx;
2041 return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
2044 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
2045 enum isl_dim_type type, unsigned pos)
2047 if (!dim)
2048 return NULL;
2050 isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
2051 isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
2053 if (type == isl_dim_set)
2054 pos += isl_space_dim(dim, isl_dim_param);
2056 return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
2057 error:
2058 isl_space_free(dim);
2059 return NULL;
2062 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
2063 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
2065 int i;
2066 struct isl_upoly_rec *rec;
2067 struct isl_upoly *base, *res;
2069 if (!up)
2070 return NULL;
2072 if (isl_upoly_is_cst(up))
2073 return up;
2075 if (up->var < first)
2076 return up;
2078 rec = isl_upoly_as_rec(up);
2079 if (!rec)
2080 goto error;
2082 isl_assert(up->ctx, rec->n >= 1, goto error);
2084 if (up->var >= first + n)
2085 base = isl_upoly_var_pow(up->ctx, up->var, 1);
2086 else
2087 base = isl_upoly_copy(subs[up->var - first]);
2089 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
2090 for (i = rec->n - 2; i >= 0; --i) {
2091 struct isl_upoly *t;
2092 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
2093 res = isl_upoly_mul(res, isl_upoly_copy(base));
2094 res = isl_upoly_sum(res, t);
2097 isl_upoly_free(base);
2098 isl_upoly_free(up);
2100 return res;
2101 error:
2102 isl_upoly_free(up);
2103 return NULL;
2106 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
2107 isl_int denom, unsigned len)
2109 int i;
2110 struct isl_upoly *up;
2112 isl_assert(ctx, len >= 1, return NULL);
2114 up = isl_upoly_rat_cst(ctx, f[0], denom);
2115 for (i = 0; i < len - 1; ++i) {
2116 struct isl_upoly *t;
2117 struct isl_upoly *c;
2119 if (isl_int_is_zero(f[1 + i]))
2120 continue;
2122 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
2123 t = isl_upoly_var_pow(ctx, i, 1);
2124 t = isl_upoly_mul(c, t);
2125 up = isl_upoly_sum(up, t);
2128 return up;
2131 /* Remove common factor of non-constant terms and denominator.
2133 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2135 isl_ctx *ctx = qp->div->ctx;
2136 unsigned total = qp->div->n_col - 2;
2138 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2139 isl_int_gcd(ctx->normalize_gcd,
2140 ctx->normalize_gcd, qp->div->row[div][0]);
2141 if (isl_int_is_one(ctx->normalize_gcd))
2142 return;
2144 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2145 ctx->normalize_gcd, total);
2146 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2147 ctx->normalize_gcd);
2148 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2149 ctx->normalize_gcd);
2152 /* Replace the integer division identified by "div" by the polynomial "s".
2153 * The integer division is assumed not to appear in the definition
2154 * of any other integer divisions.
2156 static __isl_give isl_qpolynomial *substitute_div(
2157 __isl_take isl_qpolynomial *qp,
2158 int div, __isl_take struct isl_upoly *s)
2160 int i;
2161 int total;
2162 int *reordering;
2164 if (!qp || !s)
2165 goto error;
2167 qp = isl_qpolynomial_cow(qp);
2168 if (!qp)
2169 goto error;
2171 total = isl_space_dim(qp->dim, isl_dim_all);
2172 qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
2173 if (!qp->upoly)
2174 goto error;
2176 reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
2177 if (!reordering)
2178 goto error;
2179 for (i = 0; i < total + div; ++i)
2180 reordering[i] = i;
2181 for (i = total + div + 1; i < total + qp->div->n_row; ++i)
2182 reordering[i] = i - 1;
2183 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2184 qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
2185 qp->upoly = reorder(qp->upoly, reordering);
2186 free(reordering);
2188 if (!qp->upoly || !qp->div)
2189 goto error;
2191 isl_upoly_free(s);
2192 return qp;
2193 error:
2194 isl_qpolynomial_free(qp);
2195 isl_upoly_free(s);
2196 return NULL;
2199 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2200 * divisions because d is equal to 1 by their definition, i.e., e.
2202 static __isl_give isl_qpolynomial *substitute_non_divs(
2203 __isl_take isl_qpolynomial *qp)
2205 int i, j;
2206 int total;
2207 struct isl_upoly *s;
2209 if (!qp)
2210 return NULL;
2212 total = isl_space_dim(qp->dim, isl_dim_all);
2213 for (i = 0; qp && i < qp->div->n_row; ++i) {
2214 if (!isl_int_is_one(qp->div->row[i][0]))
2215 continue;
2216 for (j = i + 1; j < qp->div->n_row; ++j) {
2217 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2218 continue;
2219 isl_seq_combine(qp->div->row[j] + 1,
2220 qp->div->ctx->one, qp->div->row[j] + 1,
2221 qp->div->row[j][2 + total + i],
2222 qp->div->row[i] + 1, 1 + total + i);
2223 isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2224 normalize_div(qp, j);
2226 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2227 qp->div->row[i][0], qp->div->n_col - 1);
2228 qp = substitute_div(qp, i, s);
2229 --i;
2232 return qp;
2235 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2236 * with d the denominator. When replacing the coefficient e of x by
2237 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2238 * inside the division, so we need to add floor(e/d) * x outside.
2239 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2240 * to adjust the coefficient of x in each later div that depends on the
2241 * current div "div" and also in the affine expression "aff"
2242 * (if it too depends on "div").
2244 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2245 __isl_keep isl_vec *aff)
2247 int i, j;
2248 isl_int v;
2249 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2251 isl_int_init(v);
2252 for (i = 0; i < 1 + total + div; ++i) {
2253 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2254 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2255 continue;
2256 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2257 isl_int_fdiv_r(qp->div->row[div][1 + i],
2258 qp->div->row[div][1 + i], qp->div->row[div][0]);
2259 if (!isl_int_is_zero(aff->el[1 + total + div]))
2260 isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
2261 for (j = div + 1; j < qp->div->n_row; ++j) {
2262 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2263 continue;
2264 isl_int_addmul(qp->div->row[j][1 + i],
2265 v, qp->div->row[j][2 + total + div]);
2268 isl_int_clear(v);
2271 /* Check if the last non-zero coefficient is bigger that half of the
2272 * denominator. If so, we will invert the div to further reduce the number
2273 * of distinct divs that may appear.
2274 * If the last non-zero coefficient is exactly half the denominator,
2275 * then we continue looking for earlier coefficients that are bigger
2276 * than half the denominator.
2278 static int needs_invert(__isl_keep isl_mat *div, int row)
2280 int i;
2281 int cmp;
2283 for (i = div->n_col - 1; i >= 1; --i) {
2284 if (isl_int_is_zero(div->row[row][i]))
2285 continue;
2286 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2287 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2288 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2289 if (cmp)
2290 return cmp > 0;
2291 if (i == 1)
2292 return 1;
2295 return 0;
2298 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2299 * We only invert the coefficients of e (and the coefficient of q in
2300 * later divs and in "aff"). After calling this function, the
2301 * coefficients of e should be reduced again.
2303 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2304 __isl_keep isl_vec *aff)
2306 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2308 isl_seq_neg(qp->div->row[div] + 1,
2309 qp->div->row[div] + 1, qp->div->n_col - 1);
2310 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2311 isl_int_add(qp->div->row[div][1],
2312 qp->div->row[div][1], qp->div->row[div][0]);
2313 if (!isl_int_is_zero(aff->el[1 + total + div]))
2314 isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
2315 isl_mat_col_mul(qp->div, 2 + total + div,
2316 qp->div->ctx->negone, 2 + total + div);
2319 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2320 * in the interval [0, d-1], with d the denominator and such that the
2321 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2323 * After the reduction, some divs may have become redundant or identical,
2324 * so we call substitute_non_divs and sort_divs. If these functions
2325 * eliminate divs or merge two or more divs into one, the coefficients
2326 * of the enclosing divs may have to be reduced again, so we call
2327 * ourselves recursively if the number of divs decreases.
2329 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2331 int i;
2332 isl_vec *aff = NULL;
2333 struct isl_upoly *s;
2334 unsigned n_div;
2336 if (!qp)
2337 return NULL;
2339 aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
2340 aff = isl_vec_clr(aff);
2341 if (!aff)
2342 goto error;
2344 isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
2346 for (i = 0; i < qp->div->n_row; ++i) {
2347 normalize_div(qp, i);
2348 reduce_div(qp, i, aff);
2349 if (needs_invert(qp->div, i)) {
2350 invert_div(qp, i, aff);
2351 reduce_div(qp, i, aff);
2355 s = isl_upoly_from_affine(qp->div->ctx, aff->el,
2356 qp->div->ctx->one, aff->size);
2357 qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
2358 isl_upoly_free(s);
2359 if (!qp->upoly)
2360 goto error;
2362 isl_vec_free(aff);
2364 n_div = qp->div->n_row;
2365 qp = substitute_non_divs(qp);
2366 qp = sort_divs(qp);
2367 if (qp && qp->div->n_row < n_div)
2368 return reduce_divs(qp);
2370 return qp;
2371 error:
2372 isl_qpolynomial_free(qp);
2373 isl_vec_free(aff);
2374 return NULL;
2377 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2378 __isl_take isl_space *dim, const isl_int n, const isl_int d)
2380 struct isl_qpolynomial *qp;
2381 struct isl_upoly_cst *cst;
2383 if (!dim)
2384 return NULL;
2386 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2387 if (!qp)
2388 return NULL;
2390 cst = isl_upoly_as_cst(qp->upoly);
2391 isl_int_set(cst->n, n);
2392 isl_int_set(cst->d, d);
2394 return qp;
2397 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2399 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2400 __isl_take isl_space *domain, __isl_take isl_val *val)
2402 isl_qpolynomial *qp;
2403 struct isl_upoly_cst *cst;
2405 if (!domain || !val)
2406 goto error;
2408 qp = isl_qpolynomial_alloc(isl_space_copy(domain), 0,
2409 isl_upoly_zero(domain->ctx));
2410 if (!qp)
2411 goto error;
2413 cst = isl_upoly_as_cst(qp->upoly);
2414 isl_int_set(cst->n, val->n);
2415 isl_int_set(cst->d, val->d);
2417 isl_space_free(domain);
2418 isl_val_free(val);
2419 return qp;
2420 error:
2421 isl_space_free(domain);
2422 isl_val_free(val);
2423 return NULL;
2426 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2428 struct isl_upoly_rec *rec;
2429 int i;
2431 if (!up)
2432 return -1;
2434 if (isl_upoly_is_cst(up))
2435 return 0;
2437 if (up->var < d)
2438 active[up->var] = 1;
2440 rec = isl_upoly_as_rec(up);
2441 for (i = 0; i < rec->n; ++i)
2442 if (up_set_active(rec->p[i], active, d) < 0)
2443 return -1;
2445 return 0;
2448 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2450 int i, j;
2451 int d = isl_space_dim(qp->dim, isl_dim_all);
2453 if (!qp || !active)
2454 return -1;
2456 for (i = 0; i < d; ++i)
2457 for (j = 0; j < qp->div->n_row; ++j) {
2458 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2459 continue;
2460 active[i] = 1;
2461 break;
2464 return up_set_active(qp->upoly, active, d);
2467 isl_bool isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2468 enum isl_dim_type type, unsigned first, unsigned n)
2470 int i;
2471 int *active = NULL;
2472 isl_bool involves = isl_bool_false;
2474 if (!qp)
2475 return isl_bool_error;
2476 if (n == 0)
2477 return isl_bool_false;
2479 isl_assert(qp->dim->ctx,
2480 first + n <= isl_qpolynomial_dim(qp, type),
2481 return isl_bool_error);
2482 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2483 type == isl_dim_in, return isl_bool_error);
2485 active = isl_calloc_array(qp->dim->ctx, int,
2486 isl_space_dim(qp->dim, isl_dim_all));
2487 if (set_active(qp, active) < 0)
2488 goto error;
2490 if (type == isl_dim_in)
2491 first += isl_space_dim(qp->dim, isl_dim_param);
2492 for (i = 0; i < n; ++i)
2493 if (active[first + i]) {
2494 involves = isl_bool_true;
2495 break;
2498 free(active);
2500 return involves;
2501 error:
2502 free(active);
2503 return isl_bool_error;
2506 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2507 * of the divs that do appear in the quasi-polynomial.
2509 static __isl_give isl_qpolynomial *remove_redundant_divs(
2510 __isl_take isl_qpolynomial *qp)
2512 int i, j;
2513 int d;
2514 int len;
2515 int skip;
2516 int *active = NULL;
2517 int *reordering = NULL;
2518 int redundant = 0;
2519 int n_div;
2520 isl_ctx *ctx;
2522 if (!qp)
2523 return NULL;
2524 if (qp->div->n_row == 0)
2525 return qp;
2527 d = isl_space_dim(qp->dim, isl_dim_all);
2528 len = qp->div->n_col - 2;
2529 ctx = isl_qpolynomial_get_ctx(qp);
2530 active = isl_calloc_array(ctx, int, len);
2531 if (!active)
2532 goto error;
2534 if (up_set_active(qp->upoly, active, len) < 0)
2535 goto error;
2537 for (i = qp->div->n_row - 1; i >= 0; --i) {
2538 if (!active[d + i]) {
2539 redundant = 1;
2540 continue;
2542 for (j = 0; j < i; ++j) {
2543 if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2544 continue;
2545 active[d + j] = 1;
2546 break;
2550 if (!redundant) {
2551 free(active);
2552 return qp;
2555 reordering = isl_alloc_array(qp->div->ctx, int, len);
2556 if (!reordering)
2557 goto error;
2559 for (i = 0; i < d; ++i)
2560 reordering[i] = i;
2562 skip = 0;
2563 n_div = qp->div->n_row;
2564 for (i = 0; i < n_div; ++i) {
2565 if (!active[d + i]) {
2566 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2567 qp->div = isl_mat_drop_cols(qp->div,
2568 2 + d + i - skip, 1);
2569 skip++;
2571 reordering[d + i] = d + i - skip;
2574 qp->upoly = reorder(qp->upoly, reordering);
2576 if (!qp->upoly || !qp->div)
2577 goto error;
2579 free(active);
2580 free(reordering);
2582 return qp;
2583 error:
2584 free(active);
2585 free(reordering);
2586 isl_qpolynomial_free(qp);
2587 return NULL;
2590 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2591 unsigned first, unsigned n)
2593 int i;
2594 struct isl_upoly_rec *rec;
2596 if (!up)
2597 return NULL;
2598 if (n == 0 || up->var < 0 || up->var < first)
2599 return up;
2600 if (up->var < first + n) {
2601 up = replace_by_constant_term(up);
2602 return isl_upoly_drop(up, first, n);
2604 up = isl_upoly_cow(up);
2605 if (!up)
2606 return NULL;
2607 up->var -= n;
2608 rec = isl_upoly_as_rec(up);
2609 if (!rec)
2610 goto error;
2612 for (i = 0; i < rec->n; ++i) {
2613 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2614 if (!rec->p[i])
2615 goto error;
2618 return up;
2619 error:
2620 isl_upoly_free(up);
2621 return NULL;
2624 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2625 __isl_take isl_qpolynomial *qp,
2626 enum isl_dim_type type, unsigned pos, const char *s)
2628 qp = isl_qpolynomial_cow(qp);
2629 if (!qp)
2630 return NULL;
2631 qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2632 if (!qp->dim)
2633 goto error;
2634 return qp;
2635 error:
2636 isl_qpolynomial_free(qp);
2637 return NULL;
2640 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2641 __isl_take isl_qpolynomial *qp,
2642 enum isl_dim_type type, unsigned first, unsigned n)
2644 if (!qp)
2645 return NULL;
2646 if (type == isl_dim_out)
2647 isl_die(qp->dim->ctx, isl_error_invalid,
2648 "cannot drop output/set dimension",
2649 goto error);
2650 if (type == isl_dim_in)
2651 type = isl_dim_set;
2652 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2653 return qp;
2655 qp = isl_qpolynomial_cow(qp);
2656 if (!qp)
2657 return NULL;
2659 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2660 goto error);
2661 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2662 type == isl_dim_set, goto error);
2664 qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2665 if (!qp->dim)
2666 goto error;
2668 if (type == isl_dim_set)
2669 first += isl_space_dim(qp->dim, isl_dim_param);
2671 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2672 if (!qp->div)
2673 goto error;
2675 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2676 if (!qp->upoly)
2677 goto error;
2679 return qp;
2680 error:
2681 isl_qpolynomial_free(qp);
2682 return NULL;
2685 /* Project the domain of the quasi-polynomial onto its parameter space.
2686 * The quasi-polynomial may not involve any of the domain dimensions.
2688 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2689 __isl_take isl_qpolynomial *qp)
2691 isl_space *space;
2692 unsigned n;
2693 int involves;
2695 n = isl_qpolynomial_dim(qp, isl_dim_in);
2696 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2697 if (involves < 0)
2698 return isl_qpolynomial_free(qp);
2699 if (involves)
2700 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2701 "polynomial involves some of the domain dimensions",
2702 return isl_qpolynomial_free(qp));
2703 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2704 space = isl_qpolynomial_get_domain_space(qp);
2705 space = isl_space_params(space);
2706 qp = isl_qpolynomial_reset_domain_space(qp, space);
2707 return qp;
2710 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2711 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2713 int i, j, k;
2714 isl_int denom;
2715 unsigned total;
2716 unsigned n_div;
2717 struct isl_upoly *up;
2719 if (!eq)
2720 goto error;
2721 if (eq->n_eq == 0) {
2722 isl_basic_set_free(eq);
2723 return qp;
2726 qp = isl_qpolynomial_cow(qp);
2727 if (!qp)
2728 goto error;
2729 qp->div = isl_mat_cow(qp->div);
2730 if (!qp->div)
2731 goto error;
2733 total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2734 n_div = eq->n_div;
2735 isl_int_init(denom);
2736 for (i = 0; i < eq->n_eq; ++i) {
2737 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2738 if (j < 0 || j == 0 || j >= total)
2739 continue;
2741 for (k = 0; k < qp->div->n_row; ++k) {
2742 if (isl_int_is_zero(qp->div->row[k][1 + j]))
2743 continue;
2744 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2745 &qp->div->row[k][0]);
2746 normalize_div(qp, k);
2749 if (isl_int_is_pos(eq->eq[i][j]))
2750 isl_seq_neg(eq->eq[i], eq->eq[i], total);
2751 isl_int_abs(denom, eq->eq[i][j]);
2752 isl_int_set_si(eq->eq[i][j], 0);
2754 up = isl_upoly_from_affine(qp->dim->ctx,
2755 eq->eq[i], denom, total);
2756 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2757 isl_upoly_free(up);
2759 isl_int_clear(denom);
2761 if (!qp->upoly)
2762 goto error;
2764 isl_basic_set_free(eq);
2766 qp = substitute_non_divs(qp);
2767 qp = sort_divs(qp);
2769 return qp;
2770 error:
2771 isl_basic_set_free(eq);
2772 isl_qpolynomial_free(qp);
2773 return NULL;
2776 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2778 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2779 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2781 if (!qp || !eq)
2782 goto error;
2783 if (qp->div->n_row > 0)
2784 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2785 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2786 error:
2787 isl_basic_set_free(eq);
2788 isl_qpolynomial_free(qp);
2789 return NULL;
2792 static __isl_give isl_basic_set *add_div_constraints(
2793 __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2795 int i;
2796 unsigned total;
2798 if (!bset || !div)
2799 goto error;
2801 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2802 if (!bset)
2803 goto error;
2804 total = isl_basic_set_total_dim(bset);
2805 for (i = 0; i < div->n_row; ++i)
2806 if (isl_basic_set_add_div_constraints_var(bset,
2807 total - div->n_row + i, div->row[i]) < 0)
2808 goto error;
2810 isl_mat_free(div);
2811 return bset;
2812 error:
2813 isl_mat_free(div);
2814 isl_basic_set_free(bset);
2815 return NULL;
2818 /* Look for equalities among the variables shared by context and qp
2819 * and the integer divisions of qp, if any.
2820 * The equalities are then used to eliminate variables and/or integer
2821 * divisions from qp.
2823 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2824 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2826 isl_basic_set *aff;
2828 if (!qp)
2829 goto error;
2830 if (qp->div->n_row > 0) {
2831 isl_basic_set *bset;
2832 context = isl_set_add_dims(context, isl_dim_set,
2833 qp->div->n_row);
2834 bset = isl_basic_set_universe(isl_set_get_space(context));
2835 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2836 context = isl_set_intersect(context,
2837 isl_set_from_basic_set(bset));
2840 aff = isl_set_affine_hull(context);
2841 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2842 error:
2843 isl_qpolynomial_free(qp);
2844 isl_set_free(context);
2845 return NULL;
2848 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2849 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2851 isl_space *space = isl_qpolynomial_get_domain_space(qp);
2852 isl_set *dom_context = isl_set_universe(space);
2853 dom_context = isl_set_intersect_params(dom_context, context);
2854 return isl_qpolynomial_gist(qp, dom_context);
2857 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2858 __isl_take isl_qpolynomial *qp)
2860 isl_set *dom;
2862 if (!qp)
2863 return NULL;
2864 if (isl_qpolynomial_is_zero(qp)) {
2865 isl_space *dim = isl_qpolynomial_get_space(qp);
2866 isl_qpolynomial_free(qp);
2867 return isl_pw_qpolynomial_zero(dim);
2870 dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2871 return isl_pw_qpolynomial_alloc(dom, qp);
2874 #undef PW
2875 #define PW isl_pw_qpolynomial
2876 #undef EL
2877 #define EL isl_qpolynomial
2878 #undef EL_IS_ZERO
2879 #define EL_IS_ZERO is_zero
2880 #undef ZERO
2881 #define ZERO zero
2882 #undef IS_ZERO
2883 #define IS_ZERO is_zero
2884 #undef FIELD
2885 #define FIELD qp
2886 #undef DEFAULT_IS_ZERO
2887 #define DEFAULT_IS_ZERO 1
2889 #define NO_PULLBACK
2891 #include <isl_pw_templ.c>
2893 #undef UNION
2894 #define UNION isl_union_pw_qpolynomial
2895 #undef PART
2896 #define PART isl_pw_qpolynomial
2897 #undef PARTS
2898 #define PARTS pw_qpolynomial
2900 #include <isl_union_single.c>
2901 #include <isl_union_eval.c>
2902 #include <isl_union_neg.c>
2904 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2906 if (!pwqp)
2907 return -1;
2909 if (pwqp->n != -1)
2910 return 0;
2912 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2913 return 0;
2915 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2918 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2919 __isl_take isl_pw_qpolynomial *pwqp1,
2920 __isl_take isl_pw_qpolynomial *pwqp2)
2922 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2925 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2926 __isl_take isl_pw_qpolynomial *pwqp1,
2927 __isl_take isl_pw_qpolynomial *pwqp2)
2929 int i, j, n;
2930 struct isl_pw_qpolynomial *res;
2932 if (!pwqp1 || !pwqp2)
2933 goto error;
2935 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2936 goto error);
2938 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2939 isl_pw_qpolynomial_free(pwqp2);
2940 return pwqp1;
2943 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2944 isl_pw_qpolynomial_free(pwqp1);
2945 return pwqp2;
2948 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2949 isl_pw_qpolynomial_free(pwqp1);
2950 return pwqp2;
2953 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2954 isl_pw_qpolynomial_free(pwqp2);
2955 return pwqp1;
2958 n = pwqp1->n * pwqp2->n;
2959 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2961 for (i = 0; i < pwqp1->n; ++i) {
2962 for (j = 0; j < pwqp2->n; ++j) {
2963 struct isl_set *common;
2964 struct isl_qpolynomial *prod;
2965 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2966 isl_set_copy(pwqp2->p[j].set));
2967 if (isl_set_plain_is_empty(common)) {
2968 isl_set_free(common);
2969 continue;
2972 prod = isl_qpolynomial_mul(
2973 isl_qpolynomial_copy(pwqp1->p[i].qp),
2974 isl_qpolynomial_copy(pwqp2->p[j].qp));
2976 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2980 isl_pw_qpolynomial_free(pwqp1);
2981 isl_pw_qpolynomial_free(pwqp2);
2983 return res;
2984 error:
2985 isl_pw_qpolynomial_free(pwqp1);
2986 isl_pw_qpolynomial_free(pwqp2);
2987 return NULL;
2990 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
2991 __isl_take isl_vec *vec)
2993 int i;
2994 struct isl_upoly_rec *rec;
2995 isl_val *res;
2996 isl_val *base;
2998 if (isl_upoly_is_cst(up)) {
2999 isl_vec_free(vec);
3000 res = isl_upoly_get_constant_val(up);
3001 isl_upoly_free(up);
3002 return res;
3005 rec = isl_upoly_as_rec(up);
3006 if (!rec)
3007 goto error;
3009 isl_assert(up->ctx, rec->n >= 1, goto error);
3011 base = isl_val_rat_from_isl_int(up->ctx,
3012 vec->el[1 + up->var], vec->el[0]);
3014 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
3015 isl_vec_copy(vec));
3017 for (i = rec->n - 2; i >= 0; --i) {
3018 res = isl_val_mul(res, isl_val_copy(base));
3019 res = isl_val_add(res,
3020 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
3021 isl_vec_copy(vec)));
3024 isl_val_free(base);
3025 isl_upoly_free(up);
3026 isl_vec_free(vec);
3027 return res;
3028 error:
3029 isl_upoly_free(up);
3030 isl_vec_free(vec);
3031 return NULL;
3034 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
3035 __isl_take isl_point *pnt)
3037 isl_vec *ext;
3038 isl_val *v;
3040 if (!qp || !pnt)
3041 goto error;
3042 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
3044 if (qp->div->n_row == 0)
3045 ext = isl_vec_copy(pnt->vec);
3046 else {
3047 int i;
3048 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
3049 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
3050 if (!ext)
3051 goto error;
3053 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
3054 for (i = 0; i < qp->div->n_row; ++i) {
3055 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
3056 1 + dim + i, &ext->el[1+dim+i]);
3057 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
3058 qp->div->row[i][0]);
3062 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
3064 isl_qpolynomial_free(qp);
3065 isl_point_free(pnt);
3067 return v;
3068 error:
3069 isl_qpolynomial_free(qp);
3070 isl_point_free(pnt);
3071 return NULL;
3074 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
3075 __isl_keep struct isl_upoly_cst *cst2)
3077 int cmp;
3078 isl_int t;
3079 isl_int_init(t);
3080 isl_int_mul(t, cst1->n, cst2->d);
3081 isl_int_submul(t, cst2->n, cst1->d);
3082 cmp = isl_int_sgn(t);
3083 isl_int_clear(t);
3084 return cmp;
3087 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
3088 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
3089 unsigned first, unsigned n)
3091 unsigned total;
3092 unsigned g_pos;
3093 int *exp;
3095 if (!qp)
3096 return NULL;
3097 if (type == isl_dim_out)
3098 isl_die(qp->div->ctx, isl_error_invalid,
3099 "cannot insert output/set dimensions",
3100 goto error);
3101 if (type == isl_dim_in)
3102 type = isl_dim_set;
3103 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3104 return qp;
3106 qp = isl_qpolynomial_cow(qp);
3107 if (!qp)
3108 return NULL;
3110 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
3111 goto error);
3113 g_pos = pos(qp->dim, type) + first;
3115 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3116 if (!qp->div)
3117 goto error;
3119 total = qp->div->n_col - 2;
3120 if (total > g_pos) {
3121 int i;
3122 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3123 if (!exp)
3124 goto error;
3125 for (i = 0; i < total - g_pos; ++i)
3126 exp[i] = i + n;
3127 qp->upoly = expand(qp->upoly, exp, g_pos);
3128 free(exp);
3129 if (!qp->upoly)
3130 goto error;
3133 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3134 if (!qp->dim)
3135 goto error;
3137 return qp;
3138 error:
3139 isl_qpolynomial_free(qp);
3140 return NULL;
3143 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3144 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3146 unsigned pos;
3148 pos = isl_qpolynomial_dim(qp, type);
3150 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3153 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3154 __isl_take isl_pw_qpolynomial *pwqp,
3155 enum isl_dim_type type, unsigned n)
3157 unsigned pos;
3159 pos = isl_pw_qpolynomial_dim(pwqp, type);
3161 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3164 static int *reordering_move(isl_ctx *ctx,
3165 unsigned len, unsigned dst, unsigned src, unsigned n)
3167 int i;
3168 int *reordering;
3170 reordering = isl_alloc_array(ctx, int, len);
3171 if (!reordering)
3172 return NULL;
3174 if (dst <= src) {
3175 for (i = 0; i < dst; ++i)
3176 reordering[i] = i;
3177 for (i = 0; i < n; ++i)
3178 reordering[src + i] = dst + i;
3179 for (i = 0; i < src - dst; ++i)
3180 reordering[dst + i] = dst + n + i;
3181 for (i = 0; i < len - src - n; ++i)
3182 reordering[src + n + i] = src + n + i;
3183 } else {
3184 for (i = 0; i < src; ++i)
3185 reordering[i] = i;
3186 for (i = 0; i < n; ++i)
3187 reordering[src + i] = dst + i;
3188 for (i = 0; i < dst - src; ++i)
3189 reordering[src + n + i] = src + i;
3190 for (i = 0; i < len - dst - n; ++i)
3191 reordering[dst + n + i] = dst + n + i;
3194 return reordering;
3197 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3198 __isl_take isl_qpolynomial *qp,
3199 enum isl_dim_type dst_type, unsigned dst_pos,
3200 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3202 unsigned g_dst_pos;
3203 unsigned g_src_pos;
3204 int *reordering;
3206 if (n == 0)
3207 return qp;
3209 qp = isl_qpolynomial_cow(qp);
3210 if (!qp)
3211 return NULL;
3213 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3214 isl_die(qp->dim->ctx, isl_error_invalid,
3215 "cannot move output/set dimension",
3216 goto error);
3217 if (dst_type == isl_dim_in)
3218 dst_type = isl_dim_set;
3219 if (src_type == isl_dim_in)
3220 src_type = isl_dim_set;
3222 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3223 goto error);
3225 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3226 g_src_pos = pos(qp->dim, src_type) + src_pos;
3227 if (dst_type > src_type)
3228 g_dst_pos -= n;
3230 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3231 if (!qp->div)
3232 goto error;
3233 qp = sort_divs(qp);
3234 if (!qp)
3235 goto error;
3237 reordering = reordering_move(qp->dim->ctx,
3238 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3239 if (!reordering)
3240 goto error;
3242 qp->upoly = reorder(qp->upoly, reordering);
3243 free(reordering);
3244 if (!qp->upoly)
3245 goto error;
3247 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3248 if (!qp->dim)
3249 goto error;
3251 return qp;
3252 error:
3253 isl_qpolynomial_free(qp);
3254 return NULL;
3257 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3258 isl_int *f, isl_int denom)
3260 struct isl_upoly *up;
3262 dim = isl_space_domain(dim);
3263 if (!dim)
3264 return NULL;
3266 up = isl_upoly_from_affine(dim->ctx, f, denom,
3267 1 + isl_space_dim(dim, isl_dim_all));
3269 return isl_qpolynomial_alloc(dim, 0, up);
3272 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3274 isl_ctx *ctx;
3275 struct isl_upoly *up;
3276 isl_qpolynomial *qp;
3278 if (!aff)
3279 return NULL;
3281 ctx = isl_aff_get_ctx(aff);
3282 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3283 aff->v->size - 1);
3285 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3286 aff->ls->div->n_row, up);
3287 if (!qp)
3288 goto error;
3290 isl_mat_free(qp->div);
3291 qp->div = isl_mat_copy(aff->ls->div);
3292 qp->div = isl_mat_cow(qp->div);
3293 if (!qp->div)
3294 goto error;
3296 isl_aff_free(aff);
3297 qp = reduce_divs(qp);
3298 qp = remove_redundant_divs(qp);
3299 return qp;
3300 error:
3301 isl_aff_free(aff);
3302 return isl_qpolynomial_free(qp);
3305 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3306 __isl_take isl_pw_aff *pwaff)
3308 int i;
3309 isl_pw_qpolynomial *pwqp;
3311 if (!pwaff)
3312 return NULL;
3314 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3315 pwaff->n);
3317 for (i = 0; i < pwaff->n; ++i) {
3318 isl_set *dom;
3319 isl_qpolynomial *qp;
3321 dom = isl_set_copy(pwaff->p[i].set);
3322 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3323 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3326 isl_pw_aff_free(pwaff);
3327 return pwqp;
3330 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3331 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3333 isl_aff *aff;
3335 aff = isl_constraint_get_bound(c, type, pos);
3336 isl_constraint_free(c);
3337 return isl_qpolynomial_from_aff(aff);
3340 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3341 * in "qp" by subs[i].
3343 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3344 __isl_take isl_qpolynomial *qp,
3345 enum isl_dim_type type, unsigned first, unsigned n,
3346 __isl_keep isl_qpolynomial **subs)
3348 int i;
3349 struct isl_upoly **ups;
3351 if (n == 0)
3352 return qp;
3354 qp = isl_qpolynomial_cow(qp);
3355 if (!qp)
3356 return NULL;
3358 if (type == isl_dim_out)
3359 isl_die(qp->dim->ctx, isl_error_invalid,
3360 "cannot substitute output/set dimension",
3361 goto error);
3362 if (type == isl_dim_in)
3363 type = isl_dim_set;
3365 for (i = 0; i < n; ++i)
3366 if (!subs[i])
3367 goto error;
3369 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3370 goto error);
3372 for (i = 0; i < n; ++i)
3373 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3374 goto error);
3376 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3377 for (i = 0; i < n; ++i)
3378 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3380 first += pos(qp->dim, type);
3382 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3383 if (!ups)
3384 goto error;
3385 for (i = 0; i < n; ++i)
3386 ups[i] = subs[i]->upoly;
3388 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3390 free(ups);
3392 if (!qp->upoly)
3393 goto error;
3395 return qp;
3396 error:
3397 isl_qpolynomial_free(qp);
3398 return NULL;
3401 /* Extend "bset" with extra set dimensions for each integer division
3402 * in "qp" and then call "fn" with the extended bset and the polynomial
3403 * that results from replacing each of the integer divisions by the
3404 * corresponding extra set dimension.
3406 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3407 __isl_keep isl_basic_set *bset,
3408 int (*fn)(__isl_take isl_basic_set *bset,
3409 __isl_take isl_qpolynomial *poly, void *user), void *user)
3411 isl_space *dim;
3412 isl_mat *div;
3413 isl_qpolynomial *poly;
3415 if (!qp || !bset)
3416 goto error;
3417 if (qp->div->n_row == 0)
3418 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3419 user);
3421 div = isl_mat_copy(qp->div);
3422 dim = isl_space_copy(qp->dim);
3423 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3424 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3425 bset = isl_basic_set_copy(bset);
3426 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3427 bset = add_div_constraints(bset, div);
3429 return fn(bset, poly, user);
3430 error:
3431 return -1;
3434 /* Return total degree in variables first (inclusive) up to last (exclusive).
3436 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3438 int deg = -1;
3439 int i;
3440 struct isl_upoly_rec *rec;
3442 if (!up)
3443 return -2;
3444 if (isl_upoly_is_zero(up))
3445 return -1;
3446 if (isl_upoly_is_cst(up) || up->var < first)
3447 return 0;
3449 rec = isl_upoly_as_rec(up);
3450 if (!rec)
3451 return -2;
3453 for (i = 0; i < rec->n; ++i) {
3454 int d;
3456 if (isl_upoly_is_zero(rec->p[i]))
3457 continue;
3458 d = isl_upoly_degree(rec->p[i], first, last);
3459 if (up->var < last)
3460 d += i;
3461 if (d > deg)
3462 deg = d;
3465 return deg;
3468 /* Return total degree in set variables.
3470 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3472 unsigned ovar;
3473 unsigned nvar;
3475 if (!poly)
3476 return -2;
3478 ovar = isl_space_offset(poly->dim, isl_dim_set);
3479 nvar = isl_space_dim(poly->dim, isl_dim_set);
3480 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3483 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3484 unsigned pos, int deg)
3486 int i;
3487 struct isl_upoly_rec *rec;
3489 if (!up)
3490 return NULL;
3492 if (isl_upoly_is_cst(up) || up->var < pos) {
3493 if (deg == 0)
3494 return isl_upoly_copy(up);
3495 else
3496 return isl_upoly_zero(up->ctx);
3499 rec = isl_upoly_as_rec(up);
3500 if (!rec)
3501 return NULL;
3503 if (up->var == pos) {
3504 if (deg < rec->n)
3505 return isl_upoly_copy(rec->p[deg]);
3506 else
3507 return isl_upoly_zero(up->ctx);
3510 up = isl_upoly_copy(up);
3511 up = isl_upoly_cow(up);
3512 rec = isl_upoly_as_rec(up);
3513 if (!rec)
3514 goto error;
3516 for (i = 0; i < rec->n; ++i) {
3517 struct isl_upoly *t;
3518 t = isl_upoly_coeff(rec->p[i], pos, deg);
3519 if (!t)
3520 goto error;
3521 isl_upoly_free(rec->p[i]);
3522 rec->p[i] = t;
3525 return up;
3526 error:
3527 isl_upoly_free(up);
3528 return NULL;
3531 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3533 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3534 __isl_keep isl_qpolynomial *qp,
3535 enum isl_dim_type type, unsigned t_pos, int deg)
3537 unsigned g_pos;
3538 struct isl_upoly *up;
3539 isl_qpolynomial *c;
3541 if (!qp)
3542 return NULL;
3544 if (type == isl_dim_out)
3545 isl_die(qp->div->ctx, isl_error_invalid,
3546 "output/set dimension does not have a coefficient",
3547 return NULL);
3548 if (type == isl_dim_in)
3549 type = isl_dim_set;
3551 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3552 return NULL);
3554 g_pos = pos(qp->dim, type) + t_pos;
3555 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3557 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3558 if (!c)
3559 return NULL;
3560 isl_mat_free(c->div);
3561 c->div = isl_mat_copy(qp->div);
3562 if (!c->div)
3563 goto error;
3564 return c;
3565 error:
3566 isl_qpolynomial_free(c);
3567 return NULL;
3570 /* Homogenize the polynomial in the variables first (inclusive) up to
3571 * last (exclusive) by inserting powers of variable first.
3572 * Variable first is assumed not to appear in the input.
3574 __isl_give struct isl_upoly *isl_upoly_homogenize(
3575 __isl_take struct isl_upoly *up, int deg, int target,
3576 int first, int last)
3578 int i;
3579 struct isl_upoly_rec *rec;
3581 if (!up)
3582 return NULL;
3583 if (isl_upoly_is_zero(up))
3584 return up;
3585 if (deg == target)
3586 return up;
3587 if (isl_upoly_is_cst(up) || up->var < first) {
3588 struct isl_upoly *hom;
3590 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3591 if (!hom)
3592 goto error;
3593 rec = isl_upoly_as_rec(hom);
3594 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3596 return hom;
3599 up = isl_upoly_cow(up);
3600 rec = isl_upoly_as_rec(up);
3601 if (!rec)
3602 goto error;
3604 for (i = 0; i < rec->n; ++i) {
3605 if (isl_upoly_is_zero(rec->p[i]))
3606 continue;
3607 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3608 up->var < last ? deg + i : i, target,
3609 first, last);
3610 if (!rec->p[i])
3611 goto error;
3614 return up;
3615 error:
3616 isl_upoly_free(up);
3617 return NULL;
3620 /* Homogenize the polynomial in the set variables by introducing
3621 * powers of an extra set variable at position 0.
3623 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3624 __isl_take isl_qpolynomial *poly)
3626 unsigned ovar;
3627 unsigned nvar;
3628 int deg = isl_qpolynomial_degree(poly);
3630 if (deg < -1)
3631 goto error;
3633 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3634 poly = isl_qpolynomial_cow(poly);
3635 if (!poly)
3636 goto error;
3638 ovar = isl_space_offset(poly->dim, isl_dim_set);
3639 nvar = isl_space_dim(poly->dim, isl_dim_set);
3640 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3641 ovar, ovar + nvar);
3642 if (!poly->upoly)
3643 goto error;
3645 return poly;
3646 error:
3647 isl_qpolynomial_free(poly);
3648 return NULL;
3651 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3652 __isl_take isl_mat *div)
3654 isl_term *term;
3655 int n;
3657 if (!dim || !div)
3658 goto error;
3660 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3662 term = isl_calloc(dim->ctx, struct isl_term,
3663 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3664 if (!term)
3665 goto error;
3667 term->ref = 1;
3668 term->dim = dim;
3669 term->div = div;
3670 isl_int_init(term->n);
3671 isl_int_init(term->d);
3673 return term;
3674 error:
3675 isl_space_free(dim);
3676 isl_mat_free(div);
3677 return NULL;
3680 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3682 if (!term)
3683 return NULL;
3685 term->ref++;
3686 return term;
3689 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3691 int i;
3692 isl_term *dup;
3693 unsigned total;
3695 if (!term)
3696 return NULL;
3698 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3700 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3701 if (!dup)
3702 return NULL;
3704 isl_int_set(dup->n, term->n);
3705 isl_int_set(dup->d, term->d);
3707 for (i = 0; i < total; ++i)
3708 dup->pow[i] = term->pow[i];
3710 return dup;
3713 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3715 if (!term)
3716 return NULL;
3718 if (term->ref == 1)
3719 return term;
3720 term->ref--;
3721 return isl_term_dup(term);
3724 void isl_term_free(__isl_take isl_term *term)
3726 if (!term)
3727 return;
3729 if (--term->ref > 0)
3730 return;
3732 isl_space_free(term->dim);
3733 isl_mat_free(term->div);
3734 isl_int_clear(term->n);
3735 isl_int_clear(term->d);
3736 free(term);
3739 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3741 if (!term)
3742 return 0;
3744 switch (type) {
3745 case isl_dim_param:
3746 case isl_dim_in:
3747 case isl_dim_out: return isl_space_dim(term->dim, type);
3748 case isl_dim_div: return term->div->n_row;
3749 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3750 term->div->n_row;
3751 default: return 0;
3755 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3757 return term ? term->dim->ctx : NULL;
3760 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3762 if (!term)
3763 return;
3764 isl_int_set(*n, term->n);
3767 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3769 if (!term)
3770 return;
3771 isl_int_set(*d, term->d);
3774 /* Return the coefficient of the term "term".
3776 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3778 if (!term)
3779 return NULL;
3781 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3782 term->n, term->d);
3785 int isl_term_get_exp(__isl_keep isl_term *term,
3786 enum isl_dim_type type, unsigned pos)
3788 if (!term)
3789 return -1;
3791 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3793 if (type >= isl_dim_set)
3794 pos += isl_space_dim(term->dim, isl_dim_param);
3795 if (type >= isl_dim_div)
3796 pos += isl_space_dim(term->dim, isl_dim_set);
3798 return term->pow[pos];
3801 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3803 isl_local_space *ls;
3804 isl_aff *aff;
3806 if (!term)
3807 return NULL;
3809 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3810 return NULL);
3812 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3813 isl_mat_copy(term->div));
3814 aff = isl_aff_alloc(ls);
3815 if (!aff)
3816 return NULL;
3818 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3820 aff = isl_aff_normalize(aff);
3822 return aff;
3825 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3826 isl_stat (*fn)(__isl_take isl_term *term, void *user),
3827 __isl_take isl_term *term, void *user)
3829 int i;
3830 struct isl_upoly_rec *rec;
3832 if (!up || !term)
3833 goto error;
3835 if (isl_upoly_is_zero(up))
3836 return term;
3838 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3839 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3840 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3842 if (isl_upoly_is_cst(up)) {
3843 struct isl_upoly_cst *cst;
3844 cst = isl_upoly_as_cst(up);
3845 if (!cst)
3846 goto error;
3847 term = isl_term_cow(term);
3848 if (!term)
3849 goto error;
3850 isl_int_set(term->n, cst->n);
3851 isl_int_set(term->d, cst->d);
3852 if (fn(isl_term_copy(term), user) < 0)
3853 goto error;
3854 return term;
3857 rec = isl_upoly_as_rec(up);
3858 if (!rec)
3859 goto error;
3861 for (i = 0; i < rec->n; ++i) {
3862 term = isl_term_cow(term);
3863 if (!term)
3864 goto error;
3865 term->pow[up->var] = i;
3866 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3867 if (!term)
3868 goto error;
3870 term->pow[up->var] = 0;
3872 return term;
3873 error:
3874 isl_term_free(term);
3875 return NULL;
3878 isl_stat isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3879 isl_stat (*fn)(__isl_take isl_term *term, void *user), void *user)
3881 isl_term *term;
3883 if (!qp)
3884 return isl_stat_error;
3886 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3887 if (!term)
3888 return isl_stat_error;
3890 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3892 isl_term_free(term);
3894 return term ? isl_stat_ok : isl_stat_error;
3897 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3899 struct isl_upoly *up;
3900 isl_qpolynomial *qp;
3901 int i, n;
3903 if (!term)
3904 return NULL;
3906 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3908 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3909 for (i = 0; i < n; ++i) {
3910 if (!term->pow[i])
3911 continue;
3912 up = isl_upoly_mul(up,
3913 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3916 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3917 if (!qp)
3918 goto error;
3919 isl_mat_free(qp->div);
3920 qp->div = isl_mat_copy(term->div);
3921 if (!qp->div)
3922 goto error;
3924 isl_term_free(term);
3925 return qp;
3926 error:
3927 isl_qpolynomial_free(qp);
3928 isl_term_free(term);
3929 return NULL;
3932 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3933 __isl_take isl_space *dim)
3935 int i;
3936 int extra;
3937 unsigned total;
3939 if (!qp || !dim)
3940 goto error;
3942 if (isl_space_is_equal(qp->dim, dim)) {
3943 isl_space_free(dim);
3944 return qp;
3947 qp = isl_qpolynomial_cow(qp);
3948 if (!qp)
3949 goto error;
3951 extra = isl_space_dim(dim, isl_dim_set) -
3952 isl_space_dim(qp->dim, isl_dim_set);
3953 total = isl_space_dim(qp->dim, isl_dim_all);
3954 if (qp->div->n_row) {
3955 int *exp;
3957 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3958 if (!exp)
3959 goto error;
3960 for (i = 0; i < qp->div->n_row; ++i)
3961 exp[i] = extra + i;
3962 qp->upoly = expand(qp->upoly, exp, total);
3963 free(exp);
3964 if (!qp->upoly)
3965 goto error;
3967 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3968 if (!qp->div)
3969 goto error;
3970 for (i = 0; i < qp->div->n_row; ++i)
3971 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3973 isl_space_free(qp->dim);
3974 qp->dim = dim;
3976 return qp;
3977 error:
3978 isl_space_free(dim);
3979 isl_qpolynomial_free(qp);
3980 return NULL;
3983 /* For each parameter or variable that does not appear in qp,
3984 * first eliminate the variable from all constraints and then set it to zero.
3986 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3987 __isl_keep isl_qpolynomial *qp)
3989 int *active = NULL;
3990 int i;
3991 int d;
3992 unsigned nparam;
3993 unsigned nvar;
3995 if (!set || !qp)
3996 goto error;
3998 d = isl_space_dim(set->dim, isl_dim_all);
3999 active = isl_calloc_array(set->ctx, int, d);
4000 if (set_active(qp, active) < 0)
4001 goto error;
4003 for (i = 0; i < d; ++i)
4004 if (!active[i])
4005 break;
4007 if (i == d) {
4008 free(active);
4009 return set;
4012 nparam = isl_space_dim(set->dim, isl_dim_param);
4013 nvar = isl_space_dim(set->dim, isl_dim_set);
4014 for (i = 0; i < nparam; ++i) {
4015 if (active[i])
4016 continue;
4017 set = isl_set_eliminate(set, isl_dim_param, i, 1);
4018 set = isl_set_fix_si(set, isl_dim_param, i, 0);
4020 for (i = 0; i < nvar; ++i) {
4021 if (active[nparam + i])
4022 continue;
4023 set = isl_set_eliminate(set, isl_dim_set, i, 1);
4024 set = isl_set_fix_si(set, isl_dim_set, i, 0);
4027 free(active);
4029 return set;
4030 error:
4031 free(active);
4032 isl_set_free(set);
4033 return NULL;
4036 struct isl_opt_data {
4037 isl_qpolynomial *qp;
4038 int first;
4039 isl_val *opt;
4040 int max;
4043 static isl_stat opt_fn(__isl_take isl_point *pnt, void *user)
4045 struct isl_opt_data *data = (struct isl_opt_data *)user;
4046 isl_val *val;
4048 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
4049 if (data->first) {
4050 data->first = 0;
4051 data->opt = val;
4052 } else if (data->max) {
4053 data->opt = isl_val_max(data->opt, val);
4054 } else {
4055 data->opt = isl_val_min(data->opt, val);
4058 return isl_stat_ok;
4061 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
4062 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
4064 struct isl_opt_data data = { NULL, 1, NULL, max };
4066 if (!set || !qp)
4067 goto error;
4069 if (isl_upoly_is_cst(qp->upoly)) {
4070 isl_set_free(set);
4071 data.opt = isl_qpolynomial_get_constant_val(qp);
4072 isl_qpolynomial_free(qp);
4073 return data.opt;
4076 set = fix_inactive(set, qp);
4078 data.qp = qp;
4079 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
4080 goto error;
4082 if (data.first)
4083 data.opt = isl_val_zero(isl_set_get_ctx(set));
4085 isl_set_free(set);
4086 isl_qpolynomial_free(qp);
4087 return data.opt;
4088 error:
4089 isl_set_free(set);
4090 isl_qpolynomial_free(qp);
4091 isl_val_free(data.opt);
4092 return NULL;
4095 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
4096 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
4098 int i;
4099 int n_sub;
4100 isl_ctx *ctx;
4101 struct isl_upoly **subs;
4102 isl_mat *mat, *diag;
4104 qp = isl_qpolynomial_cow(qp);
4105 if (!qp || !morph)
4106 goto error;
4108 ctx = qp->dim->ctx;
4109 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
4111 n_sub = morph->inv->n_row - 1;
4112 if (morph->inv->n_row != morph->inv->n_col)
4113 n_sub += qp->div->n_row;
4114 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
4115 if (n_sub && !subs)
4116 goto error;
4118 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4119 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4120 morph->inv->row[0][0], morph->inv->n_col);
4121 if (morph->inv->n_row != morph->inv->n_col)
4122 for (i = 0; i < qp->div->n_row; ++i)
4123 subs[morph->inv->n_row - 1 + i] =
4124 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4126 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4128 for (i = 0; i < n_sub; ++i)
4129 isl_upoly_free(subs[i]);
4130 free(subs);
4132 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4133 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4134 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4135 mat = isl_mat_diagonal(mat, diag);
4136 qp->div = isl_mat_product(qp->div, mat);
4137 isl_space_free(qp->dim);
4138 qp->dim = isl_space_copy(morph->ran->dim);
4140 if (!qp->upoly || !qp->div || !qp->dim)
4141 goto error;
4143 isl_morph_free(morph);
4145 return qp;
4146 error:
4147 isl_qpolynomial_free(qp);
4148 isl_morph_free(morph);
4149 return NULL;
4152 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4153 __isl_take isl_union_pw_qpolynomial *upwqp1,
4154 __isl_take isl_union_pw_qpolynomial *upwqp2)
4156 return isl_union_pw_qpolynomial_match_bin_op(upwqp1, upwqp2,
4157 &isl_pw_qpolynomial_mul);
4160 /* Reorder the columns of the given div definitions according to the
4161 * given reordering.
4163 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4164 __isl_take isl_reordering *r)
4166 int i, j;
4167 isl_mat *mat;
4168 int extra;
4170 if (!div || !r)
4171 goto error;
4173 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4174 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4175 if (!mat)
4176 goto error;
4178 for (i = 0; i < div->n_row; ++i) {
4179 isl_seq_cpy(mat->row[i], div->row[i], 2);
4180 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4181 for (j = 0; j < r->len; ++j)
4182 isl_int_set(mat->row[i][2 + r->pos[j]],
4183 div->row[i][2 + j]);
4186 isl_reordering_free(r);
4187 isl_mat_free(div);
4188 return mat;
4189 error:
4190 isl_reordering_free(r);
4191 isl_mat_free(div);
4192 return NULL;
4195 /* Reorder the dimension of "qp" according to the given reordering.
4197 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4198 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4200 qp = isl_qpolynomial_cow(qp);
4201 if (!qp)
4202 goto error;
4204 r = isl_reordering_extend(r, qp->div->n_row);
4205 if (!r)
4206 goto error;
4208 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4209 if (!qp->div)
4210 goto error;
4212 qp->upoly = reorder(qp->upoly, r->pos);
4213 if (!qp->upoly)
4214 goto error;
4216 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4218 isl_reordering_free(r);
4219 return qp;
4220 error:
4221 isl_qpolynomial_free(qp);
4222 isl_reordering_free(r);
4223 return NULL;
4226 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4227 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4229 if (!qp || !model)
4230 goto error;
4232 if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4233 isl_reordering *exp;
4235 model = isl_space_drop_dims(model, isl_dim_in,
4236 0, isl_space_dim(model, isl_dim_in));
4237 model = isl_space_drop_dims(model, isl_dim_out,
4238 0, isl_space_dim(model, isl_dim_out));
4239 exp = isl_parameter_alignment_reordering(qp->dim, model);
4240 exp = isl_reordering_extend_space(exp,
4241 isl_qpolynomial_get_domain_space(qp));
4242 qp = isl_qpolynomial_realign_domain(qp, exp);
4245 isl_space_free(model);
4246 return qp;
4247 error:
4248 isl_space_free(model);
4249 isl_qpolynomial_free(qp);
4250 return NULL;
4253 struct isl_split_periods_data {
4254 int max_periods;
4255 isl_pw_qpolynomial *res;
4258 /* Create a slice where the integer division "div" has the fixed value "v".
4259 * In particular, if "div" refers to floor(f/m), then create a slice
4261 * m v <= f <= m v + (m - 1)
4263 * or
4265 * f - m v >= 0
4266 * -f + m v + (m - 1) >= 0
4268 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4269 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4271 int total;
4272 isl_basic_set *bset = NULL;
4273 int k;
4275 if (!dim || !qp)
4276 goto error;
4278 total = isl_space_dim(dim, isl_dim_all);
4279 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4281 k = isl_basic_set_alloc_inequality(bset);
4282 if (k < 0)
4283 goto error;
4284 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4285 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4287 k = isl_basic_set_alloc_inequality(bset);
4288 if (k < 0)
4289 goto error;
4290 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4291 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4292 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4293 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4295 isl_space_free(dim);
4296 return isl_set_from_basic_set(bset);
4297 error:
4298 isl_basic_set_free(bset);
4299 isl_space_free(dim);
4300 return NULL;
4303 static isl_stat split_periods(__isl_take isl_set *set,
4304 __isl_take isl_qpolynomial *qp, void *user);
4306 /* Create a slice of the domain "set" such that integer division "div"
4307 * has the fixed value "v" and add the results to data->res,
4308 * replacing the integer division by "v" in "qp".
4310 static isl_stat set_div(__isl_take isl_set *set,
4311 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4312 struct isl_split_periods_data *data)
4314 int i;
4315 int total;
4316 isl_set *slice;
4317 struct isl_upoly *cst;
4319 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4320 set = isl_set_intersect(set, slice);
4322 if (!qp)
4323 goto error;
4325 total = isl_space_dim(qp->dim, isl_dim_all);
4327 for (i = div + 1; i < qp->div->n_row; ++i) {
4328 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4329 continue;
4330 isl_int_addmul(qp->div->row[i][1],
4331 qp->div->row[i][2 + total + div], v);
4332 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4335 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4336 qp = substitute_div(qp, div, cst);
4338 return split_periods(set, qp, data);
4339 error:
4340 isl_set_free(set);
4341 isl_qpolynomial_free(qp);
4342 return -1;
4345 /* Split the domain "set" such that integer division "div"
4346 * has a fixed value (ranging from "min" to "max") on each slice
4347 * and add the results to data->res.
4349 static isl_stat split_div(__isl_take isl_set *set,
4350 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4351 struct isl_split_periods_data *data)
4353 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4354 isl_set *set_i = isl_set_copy(set);
4355 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4357 if (set_div(set_i, qp_i, div, min, data) < 0)
4358 goto error;
4360 isl_set_free(set);
4361 isl_qpolynomial_free(qp);
4362 return isl_stat_ok;
4363 error:
4364 isl_set_free(set);
4365 isl_qpolynomial_free(qp);
4366 return isl_stat_error;
4369 /* If "qp" refers to any integer division
4370 * that can only attain "max_periods" distinct values on "set"
4371 * then split the domain along those distinct values.
4372 * Add the results (or the original if no splitting occurs)
4373 * to data->res.
4375 static isl_stat split_periods(__isl_take isl_set *set,
4376 __isl_take isl_qpolynomial *qp, void *user)
4378 int i;
4379 isl_pw_qpolynomial *pwqp;
4380 struct isl_split_periods_data *data;
4381 isl_int min, max;
4382 int total;
4383 isl_stat r = isl_stat_ok;
4385 data = (struct isl_split_periods_data *)user;
4387 if (!set || !qp)
4388 goto error;
4390 if (qp->div->n_row == 0) {
4391 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4392 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4393 return isl_stat_ok;
4396 isl_int_init(min);
4397 isl_int_init(max);
4398 total = isl_space_dim(qp->dim, isl_dim_all);
4399 for (i = 0; i < qp->div->n_row; ++i) {
4400 enum isl_lp_result lp_res;
4402 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4403 qp->div->n_row) != -1)
4404 continue;
4406 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4407 set->ctx->one, &min, NULL, NULL);
4408 if (lp_res == isl_lp_error)
4409 goto error2;
4410 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4411 continue;
4412 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4414 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4415 set->ctx->one, &max, NULL, NULL);
4416 if (lp_res == isl_lp_error)
4417 goto error2;
4418 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4419 continue;
4420 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4422 isl_int_sub(max, max, min);
4423 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4424 isl_int_add(max, max, min);
4425 break;
4429 if (i < qp->div->n_row) {
4430 r = split_div(set, qp, i, min, max, data);
4431 } else {
4432 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4433 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4436 isl_int_clear(max);
4437 isl_int_clear(min);
4439 return r;
4440 error2:
4441 isl_int_clear(max);
4442 isl_int_clear(min);
4443 error:
4444 isl_set_free(set);
4445 isl_qpolynomial_free(qp);
4446 return isl_stat_error;
4449 /* If any quasi-polynomial in pwqp refers to any integer division
4450 * that can only attain "max_periods" distinct values on its domain
4451 * then split the domain along those distinct values.
4453 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4454 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4456 struct isl_split_periods_data data;
4458 data.max_periods = max_periods;
4459 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4461 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4462 goto error;
4464 isl_pw_qpolynomial_free(pwqp);
4466 return data.res;
4467 error:
4468 isl_pw_qpolynomial_free(data.res);
4469 isl_pw_qpolynomial_free(pwqp);
4470 return NULL;
4473 /* Construct a piecewise quasipolynomial that is constant on the given
4474 * domain. In particular, it is
4475 * 0 if cst == 0
4476 * 1 if cst == 1
4477 * infinity if cst == -1
4479 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4480 __isl_take isl_basic_set *bset, int cst)
4482 isl_space *dim;
4483 isl_qpolynomial *qp;
4485 if (!bset)
4486 return NULL;
4488 bset = isl_basic_set_params(bset);
4489 dim = isl_basic_set_get_space(bset);
4490 if (cst < 0)
4491 qp = isl_qpolynomial_infty_on_domain(dim);
4492 else if (cst == 0)
4493 qp = isl_qpolynomial_zero_on_domain(dim);
4494 else
4495 qp = isl_qpolynomial_one_on_domain(dim);
4496 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4499 /* Factor bset, call fn on each of the factors and return the product.
4501 * If no factors can be found, simply call fn on the input.
4502 * Otherwise, construct the factors based on the factorizer,
4503 * call fn on each factor and compute the product.
4505 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4506 __isl_take isl_basic_set *bset,
4507 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4509 int i, n;
4510 isl_space *dim;
4511 isl_set *set;
4512 isl_factorizer *f;
4513 isl_qpolynomial *qp;
4514 isl_pw_qpolynomial *pwqp;
4515 unsigned nparam;
4516 unsigned nvar;
4518 f = isl_basic_set_factorizer(bset);
4519 if (!f)
4520 goto error;
4521 if (f->n_group == 0) {
4522 isl_factorizer_free(f);
4523 return fn(bset);
4526 nparam = isl_basic_set_dim(bset, isl_dim_param);
4527 nvar = isl_basic_set_dim(bset, isl_dim_set);
4529 dim = isl_basic_set_get_space(bset);
4530 dim = isl_space_domain(dim);
4531 set = isl_set_universe(isl_space_copy(dim));
4532 qp = isl_qpolynomial_one_on_domain(dim);
4533 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4535 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4537 for (i = 0, n = 0; i < f->n_group; ++i) {
4538 isl_basic_set *bset_i;
4539 isl_pw_qpolynomial *pwqp_i;
4541 bset_i = isl_basic_set_copy(bset);
4542 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4543 nparam + n + f->len[i], nvar - n - f->len[i]);
4544 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4545 nparam, n);
4546 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4547 n + f->len[i], nvar - n - f->len[i]);
4548 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4550 pwqp_i = fn(bset_i);
4551 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4553 n += f->len[i];
4556 isl_basic_set_free(bset);
4557 isl_factorizer_free(f);
4559 return pwqp;
4560 error:
4561 isl_basic_set_free(bset);
4562 return NULL;
4565 /* Factor bset, call fn on each of the factors and return the product.
4566 * The function is assumed to evaluate to zero on empty domains,
4567 * to one on zero-dimensional domains and to infinity on unbounded domains
4568 * and will not be called explicitly on zero-dimensional or unbounded domains.
4570 * We first check for some special cases and remove all equalities.
4571 * Then we hand over control to compressed_multiplicative_call.
4573 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4574 __isl_take isl_basic_set *bset,
4575 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4577 int bounded;
4578 isl_morph *morph;
4579 isl_pw_qpolynomial *pwqp;
4581 if (!bset)
4582 return NULL;
4584 if (isl_basic_set_plain_is_empty(bset))
4585 return constant_on_domain(bset, 0);
4587 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4588 return constant_on_domain(bset, 1);
4590 bounded = isl_basic_set_is_bounded(bset);
4591 if (bounded < 0)
4592 goto error;
4593 if (!bounded)
4594 return constant_on_domain(bset, -1);
4596 if (bset->n_eq == 0)
4597 return compressed_multiplicative_call(bset, fn);
4599 morph = isl_basic_set_full_compression(bset);
4600 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4602 pwqp = compressed_multiplicative_call(bset, fn);
4604 morph = isl_morph_dom_params(morph);
4605 morph = isl_morph_ran_params(morph);
4606 morph = isl_morph_inverse(morph);
4608 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4610 return pwqp;
4611 error:
4612 isl_basic_set_free(bset);
4613 return NULL;
4616 /* Drop all floors in "qp", turning each integer division [a/m] into
4617 * a rational division a/m. If "down" is set, then the integer division
4618 * is replaced by (a-(m-1))/m instead.
4620 static __isl_give isl_qpolynomial *qp_drop_floors(
4621 __isl_take isl_qpolynomial *qp, int down)
4623 int i;
4624 struct isl_upoly *s;
4626 if (!qp)
4627 return NULL;
4628 if (qp->div->n_row == 0)
4629 return qp;
4631 qp = isl_qpolynomial_cow(qp);
4632 if (!qp)
4633 return NULL;
4635 for (i = qp->div->n_row - 1; i >= 0; --i) {
4636 if (down) {
4637 isl_int_sub(qp->div->row[i][1],
4638 qp->div->row[i][1], qp->div->row[i][0]);
4639 isl_int_add_ui(qp->div->row[i][1],
4640 qp->div->row[i][1], 1);
4642 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4643 qp->div->row[i][0], qp->div->n_col - 1);
4644 qp = substitute_div(qp, i, s);
4645 if (!qp)
4646 return NULL;
4649 return qp;
4652 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4653 * a rational division a/m.
4655 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4656 __isl_take isl_pw_qpolynomial *pwqp)
4658 int i;
4660 if (!pwqp)
4661 return NULL;
4663 if (isl_pw_qpolynomial_is_zero(pwqp))
4664 return pwqp;
4666 pwqp = isl_pw_qpolynomial_cow(pwqp);
4667 if (!pwqp)
4668 return NULL;
4670 for (i = 0; i < pwqp->n; ++i) {
4671 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4672 if (!pwqp->p[i].qp)
4673 goto error;
4676 return pwqp;
4677 error:
4678 isl_pw_qpolynomial_free(pwqp);
4679 return NULL;
4682 /* Adjust all the integer divisions in "qp" such that they are at least
4683 * one over the given orthant (identified by "signs"). This ensures
4684 * that they will still be non-negative even after subtracting (m-1)/m.
4686 * In particular, f is replaced by f' + v, changing f = [a/m]
4687 * to f' = [(a - m v)/m].
4688 * If the constant term k in a is smaller than m,
4689 * the constant term of v is set to floor(k/m) - 1.
4690 * For any other term, if the coefficient c and the variable x have
4691 * the same sign, then no changes are needed.
4692 * Otherwise, if the variable is positive (and c is negative),
4693 * then the coefficient of x in v is set to floor(c/m).
4694 * If the variable is negative (and c is positive),
4695 * then the coefficient of x in v is set to ceil(c/m).
4697 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4698 int *signs)
4700 int i, j;
4701 int total;
4702 isl_vec *v = NULL;
4703 struct isl_upoly *s;
4705 qp = isl_qpolynomial_cow(qp);
4706 if (!qp)
4707 return NULL;
4708 qp->div = isl_mat_cow(qp->div);
4709 if (!qp->div)
4710 goto error;
4712 total = isl_space_dim(qp->dim, isl_dim_all);
4713 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4715 for (i = 0; i < qp->div->n_row; ++i) {
4716 isl_int *row = qp->div->row[i];
4717 v = isl_vec_clr(v);
4718 if (!v)
4719 goto error;
4720 if (isl_int_lt(row[1], row[0])) {
4721 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4722 isl_int_sub_ui(v->el[0], v->el[0], 1);
4723 isl_int_submul(row[1], row[0], v->el[0]);
4725 for (j = 0; j < total; ++j) {
4726 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4727 continue;
4728 if (signs[j] < 0)
4729 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4730 else
4731 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4732 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4734 for (j = 0; j < i; ++j) {
4735 if (isl_int_sgn(row[2 + total + j]) >= 0)
4736 continue;
4737 isl_int_fdiv_q(v->el[1 + total + j],
4738 row[2 + total + j], row[0]);
4739 isl_int_submul(row[2 + total + j],
4740 row[0], v->el[1 + total + j]);
4742 for (j = i + 1; j < qp->div->n_row; ++j) {
4743 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4744 continue;
4745 isl_seq_combine(qp->div->row[j] + 1,
4746 qp->div->ctx->one, qp->div->row[j] + 1,
4747 qp->div->row[j][2 + total + i], v->el, v->size);
4749 isl_int_set_si(v->el[1 + total + i], 1);
4750 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4751 qp->div->ctx->one, v->size);
4752 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4753 isl_upoly_free(s);
4754 if (!qp->upoly)
4755 goto error;
4758 isl_vec_free(v);
4759 return qp;
4760 error:
4761 isl_vec_free(v);
4762 isl_qpolynomial_free(qp);
4763 return NULL;
4766 struct isl_to_poly_data {
4767 int sign;
4768 isl_pw_qpolynomial *res;
4769 isl_qpolynomial *qp;
4772 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4773 * We first make all integer divisions positive and then split the
4774 * quasipolynomials into terms with sign data->sign (the direction
4775 * of the requested approximation) and terms with the opposite sign.
4776 * In the first set of terms, each integer division [a/m] is
4777 * overapproximated by a/m, while in the second it is underapproximated
4778 * by (a-(m-1))/m.
4780 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4781 void *user)
4783 struct isl_to_poly_data *data = user;
4784 isl_pw_qpolynomial *t;
4785 isl_qpolynomial *qp, *up, *down;
4787 qp = isl_qpolynomial_copy(data->qp);
4788 qp = make_divs_pos(qp, signs);
4790 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4791 up = qp_drop_floors(up, 0);
4792 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4793 down = qp_drop_floors(down, 1);
4795 isl_qpolynomial_free(qp);
4796 qp = isl_qpolynomial_add(up, down);
4798 t = isl_pw_qpolynomial_alloc(orthant, qp);
4799 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4801 return 0;
4804 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4805 * the polynomial will be an overapproximation. If "sign" is negative,
4806 * it will be an underapproximation. If "sign" is zero, the approximation
4807 * will lie somewhere in between.
4809 * In particular, is sign == 0, we simply drop the floors, turning
4810 * the integer divisions into rational divisions.
4811 * Otherwise, we split the domains into orthants, make all integer divisions
4812 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4813 * depending on the requested sign and the sign of the term in which
4814 * the integer division appears.
4816 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4817 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4819 int i;
4820 struct isl_to_poly_data data;
4822 if (sign == 0)
4823 return pwqp_drop_floors(pwqp);
4825 if (!pwqp)
4826 return NULL;
4828 data.sign = sign;
4829 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4831 for (i = 0; i < pwqp->n; ++i) {
4832 if (pwqp->p[i].qp->div->n_row == 0) {
4833 isl_pw_qpolynomial *t;
4834 t = isl_pw_qpolynomial_alloc(
4835 isl_set_copy(pwqp->p[i].set),
4836 isl_qpolynomial_copy(pwqp->p[i].qp));
4837 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4838 continue;
4840 data.qp = pwqp->p[i].qp;
4841 if (isl_set_foreach_orthant(pwqp->p[i].set,
4842 &to_polynomial_on_orthant, &data) < 0)
4843 goto error;
4846 isl_pw_qpolynomial_free(pwqp);
4848 return data.res;
4849 error:
4850 isl_pw_qpolynomial_free(pwqp);
4851 isl_pw_qpolynomial_free(data.res);
4852 return NULL;
4855 static __isl_give isl_pw_qpolynomial *poly_entry(
4856 __isl_take isl_pw_qpolynomial *pwqp, void *user)
4858 int *sign = user;
4860 return isl_pw_qpolynomial_to_polynomial(pwqp, *sign);
4863 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4864 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4866 return isl_union_pw_qpolynomial_transform_inplace(upwqp,
4867 &poly_entry, &sign);
4870 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4871 __isl_take isl_qpolynomial *qp)
4873 int i, k;
4874 isl_space *dim;
4875 isl_vec *aff = NULL;
4876 isl_basic_map *bmap = NULL;
4877 unsigned pos;
4878 unsigned n_div;
4880 if (!qp)
4881 return NULL;
4882 if (!isl_upoly_is_affine(qp->upoly))
4883 isl_die(qp->dim->ctx, isl_error_invalid,
4884 "input quasi-polynomial not affine", goto error);
4885 aff = isl_qpolynomial_extract_affine(qp);
4886 if (!aff)
4887 goto error;
4888 dim = isl_qpolynomial_get_space(qp);
4889 pos = 1 + isl_space_offset(dim, isl_dim_out);
4890 n_div = qp->div->n_row;
4891 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4893 for (i = 0; i < n_div; ++i) {
4894 k = isl_basic_map_alloc_div(bmap);
4895 if (k < 0)
4896 goto error;
4897 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4898 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4899 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4900 goto error;
4902 k = isl_basic_map_alloc_equality(bmap);
4903 if (k < 0)
4904 goto error;
4905 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4906 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4907 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4909 isl_vec_free(aff);
4910 isl_qpolynomial_free(qp);
4911 bmap = isl_basic_map_finalize(bmap);
4912 return bmap;
4913 error:
4914 isl_vec_free(aff);
4915 isl_qpolynomial_free(qp);
4916 isl_basic_map_free(bmap);
4917 return NULL;