isl_qpolynomial_reset_domain_space: fix memory management annotation
[isl.git] / isl_polynomial.c
blobff6e9c4faf6373bc7002dee0d0074b7ddb7c4b3e
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 #include <isl_ctx_private.h>
13 #include <isl_map_private.h>
14 #include <isl_factorization.h>
15 #include <isl_lp_private.h>
16 #include <isl_seq.h>
17 #include <isl_union_map_private.h>
18 #include <isl_constraint_private.h>
19 #include <isl_polynomial_private.h>
20 #include <isl_point_private.h>
21 #include <isl_space_private.h>
22 #include <isl_mat_private.h>
23 #include <isl_vec_private.h>
24 #include <isl_range.h>
25 #include <isl_local.h>
26 #include <isl_local_space_private.h>
27 #include <isl_aff_private.h>
28 #include <isl_val_private.h>
29 #include <isl_config.h>
31 #undef EL_BASE
32 #define EL_BASE qpolynomial
34 #include <isl_list_templ.c>
36 #undef EL_BASE
37 #define EL_BASE pw_qpolynomial
39 #include <isl_list_templ.c>
41 static unsigned pos(__isl_keep isl_space *space, enum isl_dim_type type)
43 switch (type) {
44 case isl_dim_param: return 0;
45 case isl_dim_in: return space->nparam;
46 case isl_dim_out: return space->nparam + space->n_in;
47 default: return 0;
51 isl_bool isl_poly_is_cst(__isl_keep isl_poly *poly)
53 if (!poly)
54 return isl_bool_error;
56 return isl_bool_ok(poly->var < 0);
59 __isl_keep isl_poly_cst *isl_poly_as_cst(__isl_keep isl_poly *poly)
61 if (!poly)
62 return NULL;
64 isl_assert(poly->ctx, poly->var < 0, return NULL);
66 return (isl_poly_cst *) poly;
69 __isl_keep isl_poly_rec *isl_poly_as_rec(__isl_keep isl_poly *poly)
71 if (!poly)
72 return NULL;
74 isl_assert(poly->ctx, poly->var >= 0, return NULL);
76 return (isl_poly_rec *) poly;
79 /* Compare two polynomials.
81 * Return -1 if "poly1" is "smaller" than "poly2", 1 if "poly1" is "greater"
82 * than "poly2" and 0 if they are equal.
84 static int isl_poly_plain_cmp(__isl_keep isl_poly *poly1,
85 __isl_keep isl_poly *poly2)
87 int i;
88 isl_bool is_cst1;
89 isl_poly_rec *rec1, *rec2;
91 if (poly1 == poly2)
92 return 0;
93 is_cst1 = isl_poly_is_cst(poly1);
94 if (is_cst1 < 0)
95 return -1;
96 if (!poly2)
97 return 1;
98 if (poly1->var != poly2->var)
99 return poly1->var - poly2->var;
101 if (is_cst1) {
102 isl_poly_cst *cst1, *cst2;
103 int cmp;
105 cst1 = isl_poly_as_cst(poly1);
106 cst2 = isl_poly_as_cst(poly2);
107 if (!cst1 || !cst2)
108 return 0;
109 cmp = isl_int_cmp(cst1->n, cst2->n);
110 if (cmp != 0)
111 return cmp;
112 return isl_int_cmp(cst1->d, cst2->d);
115 rec1 = isl_poly_as_rec(poly1);
116 rec2 = isl_poly_as_rec(poly2);
117 if (!rec1 || !rec2)
118 return 0;
120 if (rec1->n != rec2->n)
121 return rec1->n - rec2->n;
123 for (i = 0; i < rec1->n; ++i) {
124 int cmp = isl_poly_plain_cmp(rec1->p[i], rec2->p[i]);
125 if (cmp != 0)
126 return cmp;
129 return 0;
132 isl_bool isl_poly_is_equal(__isl_keep isl_poly *poly1,
133 __isl_keep isl_poly *poly2)
135 int i;
136 isl_bool is_cst1;
137 isl_poly_rec *rec1, *rec2;
139 is_cst1 = isl_poly_is_cst(poly1);
140 if (is_cst1 < 0 || !poly2)
141 return isl_bool_error;
142 if (poly1 == poly2)
143 return isl_bool_true;
144 if (poly1->var != poly2->var)
145 return isl_bool_false;
146 if (is_cst1) {
147 isl_poly_cst *cst1, *cst2;
148 int r;
149 cst1 = isl_poly_as_cst(poly1);
150 cst2 = isl_poly_as_cst(poly2);
151 if (!cst1 || !cst2)
152 return isl_bool_error;
153 r = isl_int_eq(cst1->n, cst2->n) &&
154 isl_int_eq(cst1->d, cst2->d);
155 return isl_bool_ok(r);
158 rec1 = isl_poly_as_rec(poly1);
159 rec2 = isl_poly_as_rec(poly2);
160 if (!rec1 || !rec2)
161 return isl_bool_error;
163 if (rec1->n != rec2->n)
164 return isl_bool_false;
166 for (i = 0; i < rec1->n; ++i) {
167 isl_bool eq = isl_poly_is_equal(rec1->p[i], rec2->p[i]);
168 if (eq < 0 || !eq)
169 return eq;
172 return isl_bool_true;
175 isl_bool isl_poly_is_zero(__isl_keep isl_poly *poly)
177 isl_bool is_cst;
178 isl_poly_cst *cst;
180 is_cst = isl_poly_is_cst(poly);
181 if (is_cst < 0 || !is_cst)
182 return is_cst;
184 cst = isl_poly_as_cst(poly);
185 if (!cst)
186 return isl_bool_error;
188 return isl_bool_ok(isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d));
191 int isl_poly_sgn(__isl_keep isl_poly *poly)
193 isl_bool is_cst;
194 isl_poly_cst *cst;
196 is_cst = isl_poly_is_cst(poly);
197 if (is_cst < 0 || !is_cst)
198 return 0;
200 cst = isl_poly_as_cst(poly);
201 if (!cst)
202 return 0;
204 return isl_int_sgn(cst->n);
207 isl_bool isl_poly_is_nan(__isl_keep isl_poly *poly)
209 isl_bool is_cst;
210 isl_poly_cst *cst;
212 is_cst = isl_poly_is_cst(poly);
213 if (is_cst < 0 || !is_cst)
214 return is_cst;
216 cst = isl_poly_as_cst(poly);
217 if (!cst)
218 return isl_bool_error;
220 return isl_bool_ok(isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d));
223 isl_bool isl_poly_is_infty(__isl_keep isl_poly *poly)
225 isl_bool is_cst;
226 isl_poly_cst *cst;
228 is_cst = isl_poly_is_cst(poly);
229 if (is_cst < 0 || !is_cst)
230 return is_cst;
232 cst = isl_poly_as_cst(poly);
233 if (!cst)
234 return isl_bool_error;
236 return isl_bool_ok(isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d));
239 isl_bool isl_poly_is_neginfty(__isl_keep isl_poly *poly)
241 isl_bool is_cst;
242 isl_poly_cst *cst;
244 is_cst = isl_poly_is_cst(poly);
245 if (is_cst < 0 || !is_cst)
246 return is_cst;
248 cst = isl_poly_as_cst(poly);
249 if (!cst)
250 return isl_bool_error;
252 return isl_bool_ok(isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d));
255 isl_bool isl_poly_is_one(__isl_keep isl_poly *poly)
257 isl_bool is_cst;
258 isl_poly_cst *cst;
259 int r;
261 is_cst = isl_poly_is_cst(poly);
262 if (is_cst < 0 || !is_cst)
263 return is_cst;
265 cst = isl_poly_as_cst(poly);
266 if (!cst)
267 return isl_bool_error;
269 r = isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
270 return isl_bool_ok(r);
273 isl_bool isl_poly_is_negone(__isl_keep isl_poly *poly)
275 isl_bool is_cst;
276 isl_poly_cst *cst;
278 is_cst = isl_poly_is_cst(poly);
279 if (is_cst < 0 || !is_cst)
280 return is_cst;
282 cst = isl_poly_as_cst(poly);
283 if (!cst)
284 return isl_bool_error;
286 return isl_bool_ok(isl_int_is_negone(cst->n) && isl_int_is_one(cst->d));
289 __isl_give isl_poly_cst *isl_poly_cst_alloc(isl_ctx *ctx)
291 isl_poly_cst *cst;
293 cst = isl_alloc_type(ctx, struct isl_poly_cst);
294 if (!cst)
295 return NULL;
297 cst->poly.ref = 1;
298 cst->poly.ctx = ctx;
299 isl_ctx_ref(ctx);
300 cst->poly.var = -1;
302 isl_int_init(cst->n);
303 isl_int_init(cst->d);
305 return cst;
308 __isl_give isl_poly *isl_poly_zero(isl_ctx *ctx)
310 isl_poly_cst *cst;
312 cst = isl_poly_cst_alloc(ctx);
313 if (!cst)
314 return NULL;
316 isl_int_set_si(cst->n, 0);
317 isl_int_set_si(cst->d, 1);
319 return &cst->poly;
322 __isl_give isl_poly *isl_poly_one(isl_ctx *ctx)
324 isl_poly_cst *cst;
326 cst = isl_poly_cst_alloc(ctx);
327 if (!cst)
328 return NULL;
330 isl_int_set_si(cst->n, 1);
331 isl_int_set_si(cst->d, 1);
333 return &cst->poly;
336 __isl_give isl_poly *isl_poly_infty(isl_ctx *ctx)
338 isl_poly_cst *cst;
340 cst = isl_poly_cst_alloc(ctx);
341 if (!cst)
342 return NULL;
344 isl_int_set_si(cst->n, 1);
345 isl_int_set_si(cst->d, 0);
347 return &cst->poly;
350 __isl_give isl_poly *isl_poly_neginfty(isl_ctx *ctx)
352 isl_poly_cst *cst;
354 cst = isl_poly_cst_alloc(ctx);
355 if (!cst)
356 return NULL;
358 isl_int_set_si(cst->n, -1);
359 isl_int_set_si(cst->d, 0);
361 return &cst->poly;
364 __isl_give isl_poly *isl_poly_nan(isl_ctx *ctx)
366 isl_poly_cst *cst;
368 cst = isl_poly_cst_alloc(ctx);
369 if (!cst)
370 return NULL;
372 isl_int_set_si(cst->n, 0);
373 isl_int_set_si(cst->d, 0);
375 return &cst->poly;
378 __isl_give isl_poly *isl_poly_rat_cst(isl_ctx *ctx, isl_int n, isl_int d)
380 isl_poly_cst *cst;
382 cst = isl_poly_cst_alloc(ctx);
383 if (!cst)
384 return NULL;
386 isl_int_set(cst->n, n);
387 isl_int_set(cst->d, d);
389 return &cst->poly;
392 __isl_give isl_poly_rec *isl_poly_alloc_rec(isl_ctx *ctx, int var, int size)
394 isl_poly_rec *rec;
396 isl_assert(ctx, var >= 0, return NULL);
397 isl_assert(ctx, size >= 0, return NULL);
398 rec = isl_calloc(ctx, struct isl_poly_rec,
399 sizeof(struct isl_poly_rec) +
400 size * sizeof(struct isl_poly *));
401 if (!rec)
402 return NULL;
404 rec->poly.ref = 1;
405 rec->poly.ctx = ctx;
406 isl_ctx_ref(ctx);
407 rec->poly.var = var;
409 rec->n = 0;
410 rec->size = size;
412 return rec;
415 /* Return the domain space of "qp".
416 * This may be either a copy or the space itself
417 * if there is only one reference to "qp".
418 * This allows the space to be modified inplace
419 * if both the quasi-polynomial and its domain space
420 * have only a single reference.
421 * The caller is not allowed to modify "qp" between this call and
422 * a subsequent call to isl_qpolynomial_restore_domain_space.
423 * The only exception is that isl_qpolynomial_free can be called instead.
425 static __isl_give isl_space *isl_qpolynomial_take_domain_space(
426 __isl_keep isl_qpolynomial *qp)
428 isl_space *space;
430 if (!qp)
431 return NULL;
432 if (qp->ref != 1)
433 return isl_qpolynomial_get_domain_space(qp);
434 space = qp->dim;
435 qp->dim = NULL;
436 return space;
439 /* Set the domain space of "qp" to "space",
440 * where the domain space of "qp" may be missing
441 * due to a preceding call to isl_qpolynomial_take_domain_space.
442 * However, in this case, "qp" only has a single reference and
443 * then the call to isl_qpolynomial_cow has no effect.
445 static __isl_give isl_qpolynomial *isl_qpolynomial_restore_domain_space(
446 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space)
448 if (!qp || !space)
449 goto error;
451 if (qp->dim == space) {
452 isl_space_free(space);
453 return qp;
456 qp = isl_qpolynomial_cow(qp);
457 if (!qp)
458 goto error;
459 isl_space_free(qp->dim);
460 qp->dim = space;
462 return qp;
463 error:
464 isl_qpolynomial_free(qp);
465 isl_space_free(space);
466 return NULL;
469 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
470 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space)
472 return isl_qpolynomial_restore_domain_space(qp, space);
475 /* Reset the space of "qp". This function is called from isl_pw_templ.c
476 * and doesn't know if the space of an element object is represented
477 * directly or through its domain. It therefore passes along both.
479 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
480 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
481 __isl_take isl_space *domain)
483 isl_space_free(space);
484 return isl_qpolynomial_reset_domain_space(qp, domain);
487 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
489 return qp ? qp->dim->ctx : NULL;
492 /* Return the domain space of "qp".
494 static __isl_keep isl_space *isl_qpolynomial_peek_domain_space(
495 __isl_keep isl_qpolynomial *qp)
497 return qp ? qp->dim : NULL;
500 /* Return a copy of the domain space of "qp".
502 __isl_give isl_space *isl_qpolynomial_get_domain_space(
503 __isl_keep isl_qpolynomial *qp)
505 return isl_space_copy(isl_qpolynomial_peek_domain_space(qp));
508 #undef TYPE
509 #define TYPE isl_qpolynomial
510 #undef PEEK_SPACE
511 #define PEEK_SPACE peek_domain_space
513 static
514 #include "isl_type_has_equal_space_bin_templ.c"
515 static
516 #include "isl_type_check_equal_space_templ.c"
518 #undef PEEK_SPACE
520 /* Return a copy of the local variables of "qp".
522 __isl_keep isl_local *isl_qpolynomial_get_local(
523 __isl_keep isl_qpolynomial *qp)
525 return qp ? isl_local_copy(qp->div) : NULL;
528 /* Return the local variables of "qp".
529 * This may be either a copy or the local variables themselves
530 * if there is only one reference to "qp".
531 * This allows the local variables to be modified in-place
532 * if both the quasi-polynomial and its local variables
533 * have only a single reference.
534 * The caller is not allowed to modify "qp" between this call and
535 * the subsequent call to isl_qpolynomial_restore_local.
536 * The only exception is that isl_qpolynomial_free can be called instead.
538 static __isl_give isl_local *isl_qpolynomial_take_local(
539 __isl_keep isl_qpolynomial *qp)
541 isl_local *local;
543 if (!qp)
544 return NULL;
545 if (qp->ref != 1)
546 return isl_qpolynomial_get_local(qp);
547 local = qp->div;
548 qp->div = NULL;
549 return local;
552 /* Set the local variables of "qp" to "local",
553 * where the local variables of "qp" may be missing
554 * due to a preceding call to isl_qpolynomial_take_local.
555 * However, in this case, "qp" only has a single reference and
556 * then the call to isl_qpolynomial_cow has no effect.
558 static __isl_give isl_qpolynomial *isl_qpolynomial_restore_local(
559 __isl_keep isl_qpolynomial *qp, __isl_take isl_local *local)
561 if (!qp || !local)
562 goto error;
564 if (qp->div == local) {
565 isl_local_free(local);
566 return qp;
569 qp = isl_qpolynomial_cow(qp);
570 if (!qp)
571 goto error;
572 isl_local_free(qp->div);
573 qp->div = local;
575 return qp;
576 error:
577 isl_qpolynomial_free(qp);
578 isl_local_free(local);
579 return NULL;
582 /* Return a copy of the local space on which "qp" is defined.
584 static __isl_give isl_local_space *isl_qpolynomial_get_domain_local_space(
585 __isl_keep isl_qpolynomial *qp)
587 isl_space *space;
588 isl_local *local;
590 if (!qp)
591 return NULL;
593 space = isl_qpolynomial_get_domain_space(qp);
594 local = isl_qpolynomial_get_local(qp);
595 return isl_local_space_alloc_div(space, local);
598 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
600 isl_space *space;
601 if (!qp)
602 return NULL;
603 space = isl_space_copy(qp->dim);
604 space = isl_space_from_domain(space);
605 space = isl_space_add_dims(space, isl_dim_out, 1);
606 return space;
609 /* Return the number of variables of the given type in the domain of "qp".
611 isl_size isl_qpolynomial_domain_dim(__isl_keep isl_qpolynomial *qp,
612 enum isl_dim_type type)
614 isl_space *space;
615 isl_size dim;
617 space = isl_qpolynomial_peek_domain_space(qp);
619 if (!space)
620 return isl_size_error;
621 if (type == isl_dim_div)
622 return qp->div->n_row;
623 dim = isl_space_dim(space, type);
624 if (dim < 0)
625 return isl_size_error;
626 if (type == isl_dim_all) {
627 isl_size n_div;
629 n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div);
630 if (n_div < 0)
631 return isl_size_error;
632 dim += n_div;
634 return dim;
637 /* Given the type of a dimension of an isl_qpolynomial,
638 * return the type of the corresponding dimension in its domain.
639 * This function is only called for "type" equal to isl_dim_in or
640 * isl_dim_param.
642 static enum isl_dim_type domain_type(enum isl_dim_type type)
644 return type == isl_dim_in ? isl_dim_set : type;
647 /* Externally, an isl_qpolynomial has a map space, but internally, the
648 * ls field corresponds to the domain of that space.
650 isl_size isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
651 enum isl_dim_type type)
653 if (!qp)
654 return isl_size_error;
655 if (type == isl_dim_out)
656 return 1;
657 type = domain_type(type);
658 return isl_qpolynomial_domain_dim(qp, type);
661 /* Return the offset of the first variable of type "type" within
662 * the variables of the domain of "qp".
664 static isl_size isl_qpolynomial_domain_var_offset(
665 __isl_keep isl_qpolynomial *qp, enum isl_dim_type type)
667 isl_space *space;
669 space = isl_qpolynomial_peek_domain_space(qp);
671 switch (type) {
672 case isl_dim_param:
673 case isl_dim_set: return isl_space_offset(space, type);
674 case isl_dim_div: return isl_space_dim(space, isl_dim_all);
675 case isl_dim_cst:
676 default:
677 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
678 "invalid dimension type", return isl_size_error);
682 /* Return the offset of the first coefficient of type "type" in
683 * the domain of "qp".
685 unsigned isl_qpolynomial_domain_offset(__isl_keep isl_qpolynomial *qp,
686 enum isl_dim_type type)
688 switch (type) {
689 case isl_dim_cst:
690 return 0;
691 case isl_dim_param:
692 case isl_dim_set:
693 case isl_dim_div:
694 return 1 + isl_qpolynomial_domain_var_offset(qp, type);
695 default:
696 return 0;
700 isl_bool isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
702 return qp ? isl_poly_is_zero(qp->poly) : isl_bool_error;
705 isl_bool isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
707 return qp ? isl_poly_is_one(qp->poly) : isl_bool_error;
710 isl_bool isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
712 return qp ? isl_poly_is_nan(qp->poly) : isl_bool_error;
715 isl_bool isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
717 return qp ? isl_poly_is_infty(qp->poly) : isl_bool_error;
720 isl_bool isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
722 return qp ? isl_poly_is_neginfty(qp->poly) : isl_bool_error;
725 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
727 return qp ? isl_poly_sgn(qp->poly) : 0;
730 static void poly_free_cst(__isl_take isl_poly_cst *cst)
732 isl_int_clear(cst->n);
733 isl_int_clear(cst->d);
736 static void poly_free_rec(__isl_take isl_poly_rec *rec)
738 int i;
740 for (i = 0; i < rec->n; ++i)
741 isl_poly_free(rec->p[i]);
744 __isl_give isl_poly *isl_poly_copy(__isl_keep isl_poly *poly)
746 if (!poly)
747 return NULL;
749 poly->ref++;
750 return poly;
753 __isl_give isl_poly *isl_poly_dup_cst(__isl_keep isl_poly *poly)
755 isl_poly_cst *cst;
756 isl_poly_cst *dup;
758 cst = isl_poly_as_cst(poly);
759 if (!cst)
760 return NULL;
762 dup = isl_poly_as_cst(isl_poly_zero(poly->ctx));
763 if (!dup)
764 return NULL;
765 isl_int_set(dup->n, cst->n);
766 isl_int_set(dup->d, cst->d);
768 return &dup->poly;
771 __isl_give isl_poly *isl_poly_dup_rec(__isl_keep isl_poly *poly)
773 int i;
774 isl_poly_rec *rec;
775 isl_poly_rec *dup;
777 rec = isl_poly_as_rec(poly);
778 if (!rec)
779 return NULL;
781 dup = isl_poly_alloc_rec(poly->ctx, poly->var, rec->n);
782 if (!dup)
783 return NULL;
785 for (i = 0; i < rec->n; ++i) {
786 dup->p[i] = isl_poly_copy(rec->p[i]);
787 if (!dup->p[i])
788 goto error;
789 dup->n++;
792 return &dup->poly;
793 error:
794 isl_poly_free(&dup->poly);
795 return NULL;
798 __isl_give isl_poly *isl_poly_dup(__isl_keep isl_poly *poly)
800 isl_bool is_cst;
802 is_cst = isl_poly_is_cst(poly);
803 if (is_cst < 0)
804 return NULL;
805 if (is_cst)
806 return isl_poly_dup_cst(poly);
807 else
808 return isl_poly_dup_rec(poly);
811 __isl_give isl_poly *isl_poly_cow(__isl_take isl_poly *poly)
813 if (!poly)
814 return NULL;
816 if (poly->ref == 1)
817 return poly;
818 poly->ref--;
819 return isl_poly_dup(poly);
822 __isl_null isl_poly *isl_poly_free(__isl_take isl_poly *poly)
824 if (!poly)
825 return NULL;
827 if (--poly->ref > 0)
828 return NULL;
830 if (poly->var < 0)
831 poly_free_cst((isl_poly_cst *) poly);
832 else
833 poly_free_rec((isl_poly_rec *) poly);
835 isl_ctx_deref(poly->ctx);
836 free(poly);
837 return NULL;
840 static void isl_poly_cst_reduce(__isl_keep isl_poly_cst *cst)
842 isl_int gcd;
844 isl_int_init(gcd);
845 isl_int_gcd(gcd, cst->n, cst->d);
846 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
847 isl_int_divexact(cst->n, cst->n, gcd);
848 isl_int_divexact(cst->d, cst->d, gcd);
850 isl_int_clear(gcd);
853 __isl_give isl_poly *isl_poly_sum_cst(__isl_take isl_poly *poly1,
854 __isl_take isl_poly *poly2)
856 isl_poly_cst *cst1;
857 isl_poly_cst *cst2;
859 poly1 = isl_poly_cow(poly1);
860 if (!poly1 || !poly2)
861 goto error;
863 cst1 = isl_poly_as_cst(poly1);
864 cst2 = isl_poly_as_cst(poly2);
866 if (isl_int_eq(cst1->d, cst2->d))
867 isl_int_add(cst1->n, cst1->n, cst2->n);
868 else {
869 isl_int_mul(cst1->n, cst1->n, cst2->d);
870 isl_int_addmul(cst1->n, cst2->n, cst1->d);
871 isl_int_mul(cst1->d, cst1->d, cst2->d);
874 isl_poly_cst_reduce(cst1);
876 isl_poly_free(poly2);
877 return poly1;
878 error:
879 isl_poly_free(poly1);
880 isl_poly_free(poly2);
881 return NULL;
884 static __isl_give isl_poly *replace_by_zero(__isl_take isl_poly *poly)
886 struct isl_ctx *ctx;
888 if (!poly)
889 return NULL;
890 ctx = poly->ctx;
891 isl_poly_free(poly);
892 return isl_poly_zero(ctx);
895 static __isl_give isl_poly *replace_by_constant_term(__isl_take isl_poly *poly)
897 isl_poly_rec *rec;
898 isl_poly *cst;
900 if (!poly)
901 return NULL;
903 rec = isl_poly_as_rec(poly);
904 if (!rec)
905 goto error;
906 cst = isl_poly_copy(rec->p[0]);
907 isl_poly_free(poly);
908 return cst;
909 error:
910 isl_poly_free(poly);
911 return NULL;
914 __isl_give isl_poly *isl_poly_sum(__isl_take isl_poly *poly1,
915 __isl_take isl_poly *poly2)
917 int i;
918 isl_bool is_zero, is_nan, is_cst;
919 isl_poly_rec *rec1, *rec2;
921 if (!poly1 || !poly2)
922 goto error;
924 is_nan = isl_poly_is_nan(poly1);
925 if (is_nan < 0)
926 goto error;
927 if (is_nan) {
928 isl_poly_free(poly2);
929 return poly1;
932 is_nan = isl_poly_is_nan(poly2);
933 if (is_nan < 0)
934 goto error;
935 if (is_nan) {
936 isl_poly_free(poly1);
937 return poly2;
940 is_zero = isl_poly_is_zero(poly1);
941 if (is_zero < 0)
942 goto error;
943 if (is_zero) {
944 isl_poly_free(poly1);
945 return poly2;
948 is_zero = isl_poly_is_zero(poly2);
949 if (is_zero < 0)
950 goto error;
951 if (is_zero) {
952 isl_poly_free(poly2);
953 return poly1;
956 if (poly1->var < poly2->var)
957 return isl_poly_sum(poly2, poly1);
959 if (poly2->var < poly1->var) {
960 isl_poly_rec *rec;
961 isl_bool is_infty;
963 is_infty = isl_poly_is_infty(poly2);
964 if (is_infty >= 0 && !is_infty)
965 is_infty = isl_poly_is_neginfty(poly2);
966 if (is_infty < 0)
967 goto error;
968 if (is_infty) {
969 isl_poly_free(poly1);
970 return poly2;
972 poly1 = isl_poly_cow(poly1);
973 rec = isl_poly_as_rec(poly1);
974 if (!rec)
975 goto error;
976 rec->p[0] = isl_poly_sum(rec->p[0], poly2);
977 if (rec->n == 1)
978 poly1 = replace_by_constant_term(poly1);
979 return poly1;
982 is_cst = isl_poly_is_cst(poly1);
983 if (is_cst < 0)
984 goto error;
985 if (is_cst)
986 return isl_poly_sum_cst(poly1, poly2);
988 rec1 = isl_poly_as_rec(poly1);
989 rec2 = isl_poly_as_rec(poly2);
990 if (!rec1 || !rec2)
991 goto error;
993 if (rec1->n < rec2->n)
994 return isl_poly_sum(poly2, poly1);
996 poly1 = isl_poly_cow(poly1);
997 rec1 = isl_poly_as_rec(poly1);
998 if (!rec1)
999 goto error;
1001 for (i = rec2->n - 1; i >= 0; --i) {
1002 isl_bool is_zero;
1004 rec1->p[i] = isl_poly_sum(rec1->p[i],
1005 isl_poly_copy(rec2->p[i]));
1006 if (!rec1->p[i])
1007 goto error;
1008 if (i != rec1->n - 1)
1009 continue;
1010 is_zero = isl_poly_is_zero(rec1->p[i]);
1011 if (is_zero < 0)
1012 goto error;
1013 if (is_zero) {
1014 isl_poly_free(rec1->p[i]);
1015 rec1->n--;
1019 if (rec1->n == 0)
1020 poly1 = replace_by_zero(poly1);
1021 else if (rec1->n == 1)
1022 poly1 = replace_by_constant_term(poly1);
1024 isl_poly_free(poly2);
1026 return poly1;
1027 error:
1028 isl_poly_free(poly1);
1029 isl_poly_free(poly2);
1030 return NULL;
1033 __isl_give isl_poly *isl_poly_cst_add_isl_int(__isl_take isl_poly *poly,
1034 isl_int v)
1036 isl_poly_cst *cst;
1038 poly = isl_poly_cow(poly);
1039 if (!poly)
1040 return NULL;
1042 cst = isl_poly_as_cst(poly);
1044 isl_int_addmul(cst->n, cst->d, v);
1046 return poly;
1049 __isl_give isl_poly *isl_poly_add_isl_int(__isl_take isl_poly *poly, isl_int v)
1051 isl_bool is_cst;
1052 isl_poly_rec *rec;
1054 is_cst = isl_poly_is_cst(poly);
1055 if (is_cst < 0)
1056 return isl_poly_free(poly);
1057 if (is_cst)
1058 return isl_poly_cst_add_isl_int(poly, v);
1060 poly = isl_poly_cow(poly);
1061 rec = isl_poly_as_rec(poly);
1062 if (!rec)
1063 goto error;
1065 rec->p[0] = isl_poly_add_isl_int(rec->p[0], v);
1066 if (!rec->p[0])
1067 goto error;
1069 return poly;
1070 error:
1071 isl_poly_free(poly);
1072 return NULL;
1075 __isl_give isl_poly *isl_poly_cst_mul_isl_int(__isl_take isl_poly *poly,
1076 isl_int v)
1078 isl_bool is_zero;
1079 isl_poly_cst *cst;
1081 is_zero = isl_poly_is_zero(poly);
1082 if (is_zero < 0)
1083 return isl_poly_free(poly);
1084 if (is_zero)
1085 return poly;
1087 poly = isl_poly_cow(poly);
1088 if (!poly)
1089 return NULL;
1091 cst = isl_poly_as_cst(poly);
1093 isl_int_mul(cst->n, cst->n, v);
1095 return poly;
1098 __isl_give isl_poly *isl_poly_mul_isl_int(__isl_take isl_poly *poly, isl_int v)
1100 int i;
1101 isl_bool is_cst;
1102 isl_poly_rec *rec;
1104 is_cst = isl_poly_is_cst(poly);
1105 if (is_cst < 0)
1106 return isl_poly_free(poly);
1107 if (is_cst)
1108 return isl_poly_cst_mul_isl_int(poly, v);
1110 poly = isl_poly_cow(poly);
1111 rec = isl_poly_as_rec(poly);
1112 if (!rec)
1113 goto error;
1115 for (i = 0; i < rec->n; ++i) {
1116 rec->p[i] = isl_poly_mul_isl_int(rec->p[i], v);
1117 if (!rec->p[i])
1118 goto error;
1121 return poly;
1122 error:
1123 isl_poly_free(poly);
1124 return NULL;
1127 /* Multiply the constant polynomial "poly" by "v".
1129 static __isl_give isl_poly *isl_poly_cst_scale_val(__isl_take isl_poly *poly,
1130 __isl_keep isl_val *v)
1132 isl_bool is_zero;
1133 isl_poly_cst *cst;
1135 is_zero = isl_poly_is_zero(poly);
1136 if (is_zero < 0)
1137 return isl_poly_free(poly);
1138 if (is_zero)
1139 return poly;
1141 poly = isl_poly_cow(poly);
1142 if (!poly)
1143 return NULL;
1145 cst = isl_poly_as_cst(poly);
1147 isl_int_mul(cst->n, cst->n, v->n);
1148 isl_int_mul(cst->d, cst->d, v->d);
1149 isl_poly_cst_reduce(cst);
1151 return poly;
1154 /* Multiply the polynomial "poly" by "v".
1156 static __isl_give isl_poly *isl_poly_scale_val(__isl_take isl_poly *poly,
1157 __isl_keep isl_val *v)
1159 int i;
1160 isl_bool is_cst;
1161 isl_poly_rec *rec;
1163 is_cst = isl_poly_is_cst(poly);
1164 if (is_cst < 0)
1165 return isl_poly_free(poly);
1166 if (is_cst)
1167 return isl_poly_cst_scale_val(poly, v);
1169 poly = isl_poly_cow(poly);
1170 rec = isl_poly_as_rec(poly);
1171 if (!rec)
1172 goto error;
1174 for (i = 0; i < rec->n; ++i) {
1175 rec->p[i] = isl_poly_scale_val(rec->p[i], v);
1176 if (!rec->p[i])
1177 goto error;
1180 return poly;
1181 error:
1182 isl_poly_free(poly);
1183 return NULL;
1186 __isl_give isl_poly *isl_poly_mul_cst(__isl_take isl_poly *poly1,
1187 __isl_take isl_poly *poly2)
1189 isl_poly_cst *cst1;
1190 isl_poly_cst *cst2;
1192 poly1 = isl_poly_cow(poly1);
1193 if (!poly1 || !poly2)
1194 goto error;
1196 cst1 = isl_poly_as_cst(poly1);
1197 cst2 = isl_poly_as_cst(poly2);
1199 isl_int_mul(cst1->n, cst1->n, cst2->n);
1200 isl_int_mul(cst1->d, cst1->d, cst2->d);
1202 isl_poly_cst_reduce(cst1);
1204 isl_poly_free(poly2);
1205 return poly1;
1206 error:
1207 isl_poly_free(poly1);
1208 isl_poly_free(poly2);
1209 return NULL;
1212 __isl_give isl_poly *isl_poly_mul_rec(__isl_take isl_poly *poly1,
1213 __isl_take isl_poly *poly2)
1215 isl_poly_rec *rec1;
1216 isl_poly_rec *rec2;
1217 isl_poly_rec *res = NULL;
1218 int i, j;
1219 int size;
1221 rec1 = isl_poly_as_rec(poly1);
1222 rec2 = isl_poly_as_rec(poly2);
1223 if (!rec1 || !rec2)
1224 goto error;
1225 size = rec1->n + rec2->n - 1;
1226 res = isl_poly_alloc_rec(poly1->ctx, poly1->var, size);
1227 if (!res)
1228 goto error;
1230 for (i = 0; i < rec1->n; ++i) {
1231 res->p[i] = isl_poly_mul(isl_poly_copy(rec2->p[0]),
1232 isl_poly_copy(rec1->p[i]));
1233 if (!res->p[i])
1234 goto error;
1235 res->n++;
1237 for (; i < size; ++i) {
1238 res->p[i] = isl_poly_zero(poly1->ctx);
1239 if (!res->p[i])
1240 goto error;
1241 res->n++;
1243 for (i = 0; i < rec1->n; ++i) {
1244 for (j = 1; j < rec2->n; ++j) {
1245 isl_poly *poly;
1246 poly = isl_poly_mul(isl_poly_copy(rec2->p[j]),
1247 isl_poly_copy(rec1->p[i]));
1248 res->p[i + j] = isl_poly_sum(res->p[i + j], poly);
1249 if (!res->p[i + j])
1250 goto error;
1254 isl_poly_free(poly1);
1255 isl_poly_free(poly2);
1257 return &res->poly;
1258 error:
1259 isl_poly_free(poly1);
1260 isl_poly_free(poly2);
1261 isl_poly_free(&res->poly);
1262 return NULL;
1265 __isl_give isl_poly *isl_poly_mul(__isl_take isl_poly *poly1,
1266 __isl_take isl_poly *poly2)
1268 isl_bool is_zero, is_nan, is_one, is_cst;
1270 if (!poly1 || !poly2)
1271 goto error;
1273 is_nan = isl_poly_is_nan(poly1);
1274 if (is_nan < 0)
1275 goto error;
1276 if (is_nan) {
1277 isl_poly_free(poly2);
1278 return poly1;
1281 is_nan = isl_poly_is_nan(poly2);
1282 if (is_nan < 0)
1283 goto error;
1284 if (is_nan) {
1285 isl_poly_free(poly1);
1286 return poly2;
1289 is_zero = isl_poly_is_zero(poly1);
1290 if (is_zero < 0)
1291 goto error;
1292 if (is_zero) {
1293 isl_poly_free(poly2);
1294 return poly1;
1297 is_zero = isl_poly_is_zero(poly2);
1298 if (is_zero < 0)
1299 goto error;
1300 if (is_zero) {
1301 isl_poly_free(poly1);
1302 return poly2;
1305 is_one = isl_poly_is_one(poly1);
1306 if (is_one < 0)
1307 goto error;
1308 if (is_one) {
1309 isl_poly_free(poly1);
1310 return poly2;
1313 is_one = isl_poly_is_one(poly2);
1314 if (is_one < 0)
1315 goto error;
1316 if (is_one) {
1317 isl_poly_free(poly2);
1318 return poly1;
1321 if (poly1->var < poly2->var)
1322 return isl_poly_mul(poly2, poly1);
1324 if (poly2->var < poly1->var) {
1325 int i;
1326 isl_poly_rec *rec;
1327 isl_bool is_infty;
1329 is_infty = isl_poly_is_infty(poly2);
1330 if (is_infty >= 0 && !is_infty)
1331 is_infty = isl_poly_is_neginfty(poly2);
1332 if (is_infty < 0)
1333 goto error;
1334 if (is_infty) {
1335 isl_ctx *ctx = poly1->ctx;
1336 isl_poly_free(poly1);
1337 isl_poly_free(poly2);
1338 return isl_poly_nan(ctx);
1340 poly1 = isl_poly_cow(poly1);
1341 rec = isl_poly_as_rec(poly1);
1342 if (!rec)
1343 goto error;
1345 for (i = 0; i < rec->n; ++i) {
1346 rec->p[i] = isl_poly_mul(rec->p[i],
1347 isl_poly_copy(poly2));
1348 if (!rec->p[i])
1349 goto error;
1351 isl_poly_free(poly2);
1352 return poly1;
1355 is_cst = isl_poly_is_cst(poly1);
1356 if (is_cst < 0)
1357 goto error;
1358 if (is_cst)
1359 return isl_poly_mul_cst(poly1, poly2);
1361 return isl_poly_mul_rec(poly1, poly2);
1362 error:
1363 isl_poly_free(poly1);
1364 isl_poly_free(poly2);
1365 return NULL;
1368 __isl_give isl_poly *isl_poly_pow(__isl_take isl_poly *poly, unsigned power)
1370 isl_poly *res;
1372 if (!poly)
1373 return NULL;
1374 if (power == 1)
1375 return poly;
1377 if (power % 2)
1378 res = isl_poly_copy(poly);
1379 else
1380 res = isl_poly_one(poly->ctx);
1382 while (power >>= 1) {
1383 poly = isl_poly_mul(poly, isl_poly_copy(poly));
1384 if (power % 2)
1385 res = isl_poly_mul(res, isl_poly_copy(poly));
1388 isl_poly_free(poly);
1389 return res;
1392 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *space,
1393 unsigned n_div, __isl_take isl_poly *poly)
1395 struct isl_qpolynomial *qp = NULL;
1396 isl_size total;
1398 total = isl_space_dim(space, isl_dim_all);
1399 if (total < 0 || !poly)
1400 goto error;
1402 if (!isl_space_is_set(space))
1403 isl_die(isl_space_get_ctx(space), isl_error_invalid,
1404 "domain of polynomial should be a set", goto error);
1406 qp = isl_calloc_type(space->ctx, struct isl_qpolynomial);
1407 if (!qp)
1408 goto error;
1410 qp->ref = 1;
1411 qp->div = isl_mat_alloc(space->ctx, n_div, 1 + 1 + total + n_div);
1412 if (!qp->div)
1413 goto error;
1415 qp->dim = space;
1416 qp->poly = poly;
1418 return qp;
1419 error:
1420 isl_space_free(space);
1421 isl_poly_free(poly);
1422 isl_qpolynomial_free(qp);
1423 return NULL;
1426 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1428 if (!qp)
1429 return NULL;
1431 qp->ref++;
1432 return qp;
1435 /* Return a copy of the polynomial expression of "qp".
1437 __isl_give isl_poly *isl_qpolynomial_get_poly(__isl_keep isl_qpolynomial *qp)
1439 return qp ? isl_poly_copy(qp->poly) : NULL;
1442 /* Return the polynomial expression of "qp".
1443 * This may be either a copy or the polynomial expression itself
1444 * if there is only one reference to "qp".
1445 * This allows the polynomial expression to be modified inplace
1446 * if both the quasi-polynomial and its polynomial expression
1447 * have only a single reference.
1448 * The caller is not allowed to modify "qp" between this call and
1449 * a subsequent call to isl_qpolynomial_restore_poly.
1450 * The only exception is that isl_qpolynomial_free can be called instead.
1452 static __isl_give isl_poly *isl_qpolynomial_take_poly(
1453 __isl_keep isl_qpolynomial *qp)
1455 isl_poly *poly;
1457 if (!qp)
1458 return NULL;
1459 if (qp->ref != 1)
1460 return isl_qpolynomial_get_poly(qp);
1461 poly = qp->poly;
1462 qp->poly = NULL;
1463 return poly;
1466 /* Set the polynomial expression of "qp" to "space",
1467 * where the polynomial expression of "qp" may be missing
1468 * due to a preceding call to isl_qpolynomial_take_poly.
1469 * However, in this case, "qp" only has a single reference and
1470 * then the call to isl_qpolynomial_cow has no effect.
1472 static __isl_give isl_qpolynomial *isl_qpolynomial_restore_poly(
1473 __isl_keep isl_qpolynomial *qp, __isl_take isl_poly *poly)
1475 if (!qp || !poly)
1476 goto error;
1478 if (qp->poly == poly) {
1479 isl_poly_free(poly);
1480 return qp;
1483 qp = isl_qpolynomial_cow(qp);
1484 if (!qp)
1485 goto error;
1486 isl_poly_free(qp->poly);
1487 qp->poly = poly;
1489 return qp;
1490 error:
1491 isl_qpolynomial_free(qp);
1492 isl_poly_free(poly);
1493 return NULL;
1496 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1498 isl_poly *poly;
1499 struct isl_qpolynomial *dup;
1501 if (!qp)
1502 return NULL;
1504 poly = isl_qpolynomial_get_poly(qp);
1505 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1506 poly);
1507 if (!dup)
1508 return NULL;
1509 isl_mat_free(dup->div);
1510 dup->div = isl_qpolynomial_get_local(qp);
1511 if (!dup->div)
1512 goto error;
1514 return dup;
1515 error:
1516 isl_qpolynomial_free(dup);
1517 return NULL;
1520 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1522 if (!qp)
1523 return NULL;
1525 if (qp->ref == 1)
1526 return qp;
1527 qp->ref--;
1528 return isl_qpolynomial_dup(qp);
1531 __isl_null isl_qpolynomial *isl_qpolynomial_free(
1532 __isl_take isl_qpolynomial *qp)
1534 if (!qp)
1535 return NULL;
1537 if (--qp->ref > 0)
1538 return NULL;
1540 isl_space_free(qp->dim);
1541 isl_mat_free(qp->div);
1542 isl_poly_free(qp->poly);
1544 free(qp);
1545 return NULL;
1548 __isl_give isl_poly *isl_poly_var_pow(isl_ctx *ctx, int pos, int power)
1550 int i;
1551 isl_poly_rec *rec;
1552 isl_poly_cst *cst;
1554 rec = isl_poly_alloc_rec(ctx, pos, 1 + power);
1555 if (!rec)
1556 return NULL;
1557 for (i = 0; i < 1 + power; ++i) {
1558 rec->p[i] = isl_poly_zero(ctx);
1559 if (!rec->p[i])
1560 goto error;
1561 rec->n++;
1563 cst = isl_poly_as_cst(rec->p[power]);
1564 isl_int_set_si(cst->n, 1);
1566 return &rec->poly;
1567 error:
1568 isl_poly_free(&rec->poly);
1569 return NULL;
1572 /* r array maps original positions to new positions.
1574 static __isl_give isl_poly *reorder(__isl_take isl_poly *poly, int *r)
1576 int i;
1577 isl_bool is_cst;
1578 isl_poly_rec *rec;
1579 isl_poly *base;
1580 isl_poly *res;
1582 is_cst = isl_poly_is_cst(poly);
1583 if (is_cst < 0)
1584 return isl_poly_free(poly);
1585 if (is_cst)
1586 return poly;
1588 rec = isl_poly_as_rec(poly);
1589 if (!rec)
1590 goto error;
1592 isl_assert(poly->ctx, rec->n >= 1, goto error);
1594 base = isl_poly_var_pow(poly->ctx, r[poly->var], 1);
1595 res = reorder(isl_poly_copy(rec->p[rec->n - 1]), r);
1597 for (i = rec->n - 2; i >= 0; --i) {
1598 res = isl_poly_mul(res, isl_poly_copy(base));
1599 res = isl_poly_sum(res, reorder(isl_poly_copy(rec->p[i]), r));
1602 isl_poly_free(base);
1603 isl_poly_free(poly);
1605 return res;
1606 error:
1607 isl_poly_free(poly);
1608 return NULL;
1611 static isl_bool compatible_divs(__isl_keep isl_mat *div1,
1612 __isl_keep isl_mat *div2)
1614 int n_row, n_col;
1615 isl_bool equal;
1617 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1618 div1->n_col >= div2->n_col,
1619 return isl_bool_error);
1621 if (div1->n_row == div2->n_row)
1622 return isl_mat_is_equal(div1, div2);
1624 n_row = div1->n_row;
1625 n_col = div1->n_col;
1626 div1->n_row = div2->n_row;
1627 div1->n_col = div2->n_col;
1629 equal = isl_mat_is_equal(div1, div2);
1631 div1->n_row = n_row;
1632 div1->n_col = n_col;
1634 return equal;
1637 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1639 int li, lj;
1641 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1642 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1644 if (li != lj)
1645 return li - lj;
1647 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1650 struct isl_div_sort_info {
1651 isl_mat *div;
1652 int row;
1655 static int div_sort_cmp(const void *p1, const void *p2)
1657 const struct isl_div_sort_info *i1, *i2;
1658 i1 = (const struct isl_div_sort_info *) p1;
1659 i2 = (const struct isl_div_sort_info *) p2;
1661 return cmp_row(i1->div, i1->row, i2->row);
1664 /* Sort divs and remove duplicates.
1666 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1668 int i;
1669 int skip;
1670 int len;
1671 struct isl_div_sort_info *array = NULL;
1672 int *pos = NULL, *at = NULL;
1673 int *reordering = NULL;
1674 isl_size div_pos;
1676 if (!qp)
1677 return NULL;
1678 if (qp->div->n_row <= 1)
1679 return qp;
1681 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
1682 if (div_pos < 0)
1683 return isl_qpolynomial_free(qp);
1685 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1686 qp->div->n_row);
1687 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1688 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1689 len = qp->div->n_col - 2;
1690 reordering = isl_alloc_array(qp->div->ctx, int, len);
1691 if (!array || !pos || !at || !reordering)
1692 goto error;
1694 for (i = 0; i < qp->div->n_row; ++i) {
1695 array[i].div = qp->div;
1696 array[i].row = i;
1697 pos[i] = i;
1698 at[i] = i;
1701 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1702 div_sort_cmp);
1704 for (i = 0; i < div_pos; ++i)
1705 reordering[i] = i;
1707 for (i = 0; i < qp->div->n_row; ++i) {
1708 if (pos[array[i].row] == i)
1709 continue;
1710 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1711 pos[at[i]] = pos[array[i].row];
1712 at[pos[array[i].row]] = at[i];
1713 at[i] = array[i].row;
1714 pos[array[i].row] = i;
1717 skip = 0;
1718 for (i = 0; i < len - div_pos; ++i) {
1719 if (i > 0 &&
1720 isl_seq_eq(qp->div->row[i - skip - 1],
1721 qp->div->row[i - skip], qp->div->n_col)) {
1722 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1723 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1724 2 + div_pos + i - skip);
1725 qp->div = isl_mat_drop_cols(qp->div,
1726 2 + div_pos + i - skip, 1);
1727 skip++;
1729 reordering[div_pos + array[i].row] = div_pos + i - skip;
1732 qp->poly = reorder(qp->poly, reordering);
1734 if (!qp->poly || !qp->div)
1735 goto error;
1737 free(at);
1738 free(pos);
1739 free(array);
1740 free(reordering);
1742 return qp;
1743 error:
1744 free(at);
1745 free(pos);
1746 free(array);
1747 free(reordering);
1748 isl_qpolynomial_free(qp);
1749 return NULL;
1752 static __isl_give isl_poly *expand(__isl_take isl_poly *poly, int *exp,
1753 int first)
1755 int i;
1756 isl_bool is_cst;
1757 isl_poly_rec *rec;
1759 is_cst = isl_poly_is_cst(poly);
1760 if (is_cst < 0)
1761 return isl_poly_free(poly);
1762 if (is_cst)
1763 return poly;
1765 if (poly->var < first)
1766 return poly;
1768 if (exp[poly->var - first] == poly->var - first)
1769 return poly;
1771 poly = isl_poly_cow(poly);
1772 if (!poly)
1773 goto error;
1775 poly->var = exp[poly->var - first] + first;
1777 rec = isl_poly_as_rec(poly);
1778 if (!rec)
1779 goto error;
1781 for (i = 0; i < rec->n; ++i) {
1782 rec->p[i] = expand(rec->p[i], exp, first);
1783 if (!rec->p[i])
1784 goto error;
1787 return poly;
1788 error:
1789 isl_poly_free(poly);
1790 return NULL;
1793 static __isl_give isl_qpolynomial *with_merged_divs(
1794 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1795 __isl_take isl_qpolynomial *qp2),
1796 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1798 int *exp1 = NULL;
1799 int *exp2 = NULL;
1800 isl_mat *div = NULL;
1801 int n_div1, n_div2;
1803 qp1 = isl_qpolynomial_cow(qp1);
1804 qp2 = isl_qpolynomial_cow(qp2);
1806 if (!qp1 || !qp2)
1807 goto error;
1809 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1810 qp1->div->n_col >= qp2->div->n_col, goto error);
1812 n_div1 = qp1->div->n_row;
1813 n_div2 = qp2->div->n_row;
1814 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1);
1815 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2);
1816 if ((n_div1 && !exp1) || (n_div2 && !exp2))
1817 goto error;
1819 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1820 if (!div)
1821 goto error;
1823 isl_mat_free(qp1->div);
1824 qp1->div = isl_mat_copy(div);
1825 isl_mat_free(qp2->div);
1826 qp2->div = isl_mat_copy(div);
1828 qp1->poly = expand(qp1->poly, exp1, div->n_col - div->n_row - 2);
1829 qp2->poly = expand(qp2->poly, exp2, div->n_col - div->n_row - 2);
1831 if (!qp1->poly || !qp2->poly)
1832 goto error;
1834 isl_mat_free(div);
1835 free(exp1);
1836 free(exp2);
1838 return fn(qp1, qp2);
1839 error:
1840 isl_mat_free(div);
1841 free(exp1);
1842 free(exp2);
1843 isl_qpolynomial_free(qp1);
1844 isl_qpolynomial_free(qp2);
1845 return NULL;
1848 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1849 __isl_take isl_qpolynomial *qp2)
1851 isl_bool compatible;
1852 isl_poly *poly;
1854 if (isl_qpolynomial_check_equal_space(qp1, qp2) < 0)
1855 goto error;
1857 if (qp1->div->n_row < qp2->div->n_row)
1858 return isl_qpolynomial_add(qp2, qp1);
1860 compatible = compatible_divs(qp1->div, qp2->div);
1861 if (compatible < 0)
1862 goto error;
1863 if (!compatible)
1864 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1866 poly = isl_qpolynomial_take_poly(qp1);
1867 poly = isl_poly_sum(poly, isl_qpolynomial_get_poly(qp2));
1868 qp1 = isl_qpolynomial_restore_poly(qp1, poly);
1870 isl_qpolynomial_free(qp2);
1872 return qp1;
1873 error:
1874 isl_qpolynomial_free(qp1);
1875 isl_qpolynomial_free(qp2);
1876 return NULL;
1879 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1880 __isl_keep isl_set *dom,
1881 __isl_take isl_qpolynomial *qp1,
1882 __isl_take isl_qpolynomial *qp2)
1884 qp1 = isl_qpolynomial_add(qp1, qp2);
1885 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1886 return qp1;
1889 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1890 __isl_take isl_qpolynomial *qp2)
1892 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1895 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1896 __isl_take isl_qpolynomial *qp, isl_int v)
1898 isl_poly *poly;
1900 if (isl_int_is_zero(v))
1901 return qp;
1903 poly = isl_qpolynomial_take_poly(qp);
1904 poly = isl_poly_add_isl_int(poly, v);
1905 qp = isl_qpolynomial_restore_poly(qp, poly);
1907 return qp;
1910 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1912 if (!qp)
1913 return NULL;
1915 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1918 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1919 __isl_take isl_qpolynomial *qp, isl_int v)
1921 isl_poly *poly;
1923 if (isl_int_is_one(v))
1924 return qp;
1926 if (qp && isl_int_is_zero(v)) {
1927 isl_qpolynomial *zero;
1928 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1929 isl_qpolynomial_free(qp);
1930 return zero;
1933 poly = isl_qpolynomial_take_poly(qp);
1934 poly = isl_poly_mul_isl_int(poly, v);
1935 qp = isl_qpolynomial_restore_poly(qp, poly);
1937 return qp;
1940 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1941 __isl_take isl_qpolynomial *qp, isl_int v)
1943 return isl_qpolynomial_mul_isl_int(qp, v);
1946 /* Multiply "qp" by "v".
1948 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1949 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1951 isl_poly *poly;
1953 if (!qp || !v)
1954 goto error;
1956 if (!isl_val_is_rat(v))
1957 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1958 "expecting rational factor", goto error);
1960 if (isl_val_is_one(v)) {
1961 isl_val_free(v);
1962 return qp;
1965 if (isl_val_is_zero(v)) {
1966 isl_space *space;
1968 space = isl_qpolynomial_get_domain_space(qp);
1969 isl_qpolynomial_free(qp);
1970 isl_val_free(v);
1971 return isl_qpolynomial_zero_on_domain(space);
1974 poly = isl_qpolynomial_take_poly(qp);
1975 poly = isl_poly_scale_val(poly, v);
1976 qp = isl_qpolynomial_restore_poly(qp, poly);
1978 isl_val_free(v);
1979 return qp;
1980 error:
1981 isl_val_free(v);
1982 isl_qpolynomial_free(qp);
1983 return NULL;
1986 /* Divide "qp" by "v".
1988 __isl_give isl_qpolynomial *isl_qpolynomial_scale_down_val(
1989 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1991 if (!qp || !v)
1992 goto error;
1994 if (!isl_val_is_rat(v))
1995 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1996 "expecting rational factor", goto error);
1997 if (isl_val_is_zero(v))
1998 isl_die(isl_val_get_ctx(v), isl_error_invalid,
1999 "cannot scale down by zero", goto error);
2001 return isl_qpolynomial_scale_val(qp, isl_val_inv(v));
2002 error:
2003 isl_val_free(v);
2004 isl_qpolynomial_free(qp);
2005 return NULL;
2008 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
2009 __isl_take isl_qpolynomial *qp2)
2011 isl_bool compatible;
2012 isl_poly *poly;
2014 if (isl_qpolynomial_check_equal_space(qp1, qp2) < 0)
2015 goto error;
2017 if (qp1->div->n_row < qp2->div->n_row)
2018 return isl_qpolynomial_mul(qp2, qp1);
2020 compatible = compatible_divs(qp1->div, qp2->div);
2021 if (compatible < 0)
2022 goto error;
2023 if (!compatible)
2024 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
2026 poly = isl_qpolynomial_take_poly(qp1);
2027 poly = isl_poly_mul(poly, isl_qpolynomial_get_poly(qp2));
2028 qp1 = isl_qpolynomial_restore_poly(qp1, poly);
2030 isl_qpolynomial_free(qp2);
2032 return qp1;
2033 error:
2034 isl_qpolynomial_free(qp1);
2035 isl_qpolynomial_free(qp2);
2036 return NULL;
2039 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
2040 unsigned power)
2042 isl_poly *poly;
2044 poly = isl_qpolynomial_take_poly(qp);
2045 poly = isl_poly_pow(poly, power);
2046 qp = isl_qpolynomial_restore_poly(qp, poly);
2048 return qp;
2051 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
2052 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
2054 int i;
2056 if (power == 1)
2057 return pwqp;
2059 pwqp = isl_pw_qpolynomial_cow(pwqp);
2060 if (!pwqp)
2061 return NULL;
2063 for (i = 0; i < pwqp->n; ++i) {
2064 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
2065 if (!pwqp->p[i].qp)
2066 return isl_pw_qpolynomial_free(pwqp);
2069 return pwqp;
2072 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
2073 __isl_take isl_space *domain)
2075 if (!domain)
2076 return NULL;
2077 return isl_qpolynomial_alloc(domain, 0, isl_poly_zero(domain->ctx));
2080 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
2081 __isl_take isl_space *domain)
2083 if (!domain)
2084 return NULL;
2085 return isl_qpolynomial_alloc(domain, 0, isl_poly_one(domain->ctx));
2088 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
2089 __isl_take isl_space *domain)
2091 if (!domain)
2092 return NULL;
2093 return isl_qpolynomial_alloc(domain, 0, isl_poly_infty(domain->ctx));
2096 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
2097 __isl_take isl_space *domain)
2099 if (!domain)
2100 return NULL;
2101 return isl_qpolynomial_alloc(domain, 0, isl_poly_neginfty(domain->ctx));
2104 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
2105 __isl_take isl_space *domain)
2107 if (!domain)
2108 return NULL;
2109 return isl_qpolynomial_alloc(domain, 0, isl_poly_nan(domain->ctx));
2112 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
2113 __isl_take isl_space *domain,
2114 isl_int v)
2116 struct isl_qpolynomial *qp;
2117 isl_poly_cst *cst;
2119 qp = isl_qpolynomial_zero_on_domain(domain);
2120 if (!qp)
2121 return NULL;
2123 cst = isl_poly_as_cst(qp->poly);
2124 isl_int_set(cst->n, v);
2126 return qp;
2129 isl_bool isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
2130 isl_int *n, isl_int *d)
2132 isl_bool is_cst;
2133 isl_poly_cst *cst;
2135 if (!qp)
2136 return isl_bool_error;
2138 is_cst = isl_poly_is_cst(qp->poly);
2139 if (is_cst < 0 || !is_cst)
2140 return is_cst;
2142 cst = isl_poly_as_cst(qp->poly);
2143 if (!cst)
2144 return isl_bool_error;
2146 if (n)
2147 isl_int_set(*n, cst->n);
2148 if (d)
2149 isl_int_set(*d, cst->d);
2151 return isl_bool_true;
2154 /* Return the constant term of "poly".
2156 static __isl_give isl_val *isl_poly_get_constant_val(__isl_keep isl_poly *poly)
2158 isl_bool is_cst;
2159 isl_poly_cst *cst;
2161 if (!poly)
2162 return NULL;
2164 while ((is_cst = isl_poly_is_cst(poly)) == isl_bool_false) {
2165 isl_poly_rec *rec;
2167 rec = isl_poly_as_rec(poly);
2168 if (!rec)
2169 return NULL;
2170 poly = rec->p[0];
2172 if (is_cst < 0)
2173 return NULL;
2175 cst = isl_poly_as_cst(poly);
2176 if (!cst)
2177 return NULL;
2178 return isl_val_rat_from_isl_int(cst->poly.ctx, cst->n, cst->d);
2181 /* Return the constant term of "qp".
2183 __isl_give isl_val *isl_qpolynomial_get_constant_val(
2184 __isl_keep isl_qpolynomial *qp)
2186 if (!qp)
2187 return NULL;
2189 return isl_poly_get_constant_val(qp->poly);
2192 isl_bool isl_poly_is_affine(__isl_keep isl_poly *poly)
2194 isl_bool is_cst;
2195 isl_poly_rec *rec;
2197 if (!poly)
2198 return isl_bool_error;
2200 if (poly->var < 0)
2201 return isl_bool_true;
2203 rec = isl_poly_as_rec(poly);
2204 if (!rec)
2205 return isl_bool_error;
2207 if (rec->n > 2)
2208 return isl_bool_false;
2210 isl_assert(poly->ctx, rec->n > 1, return isl_bool_error);
2212 is_cst = isl_poly_is_cst(rec->p[1]);
2213 if (is_cst < 0 || !is_cst)
2214 return is_cst;
2216 return isl_poly_is_affine(rec->p[0]);
2219 isl_bool isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
2221 if (!qp)
2222 return isl_bool_error;
2224 if (qp->div->n_row > 0)
2225 return isl_bool_false;
2227 return isl_poly_is_affine(qp->poly);
2230 static void update_coeff(__isl_keep isl_vec *aff,
2231 __isl_keep isl_poly_cst *cst, int pos)
2233 isl_int gcd;
2234 isl_int f;
2236 if (isl_int_is_zero(cst->n))
2237 return;
2239 isl_int_init(gcd);
2240 isl_int_init(f);
2241 isl_int_gcd(gcd, cst->d, aff->el[0]);
2242 isl_int_divexact(f, cst->d, gcd);
2243 isl_int_divexact(gcd, aff->el[0], gcd);
2244 isl_seq_scale(aff->el, aff->el, f, aff->size);
2245 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
2246 isl_int_clear(gcd);
2247 isl_int_clear(f);
2250 int isl_poly_update_affine(__isl_keep isl_poly *poly, __isl_keep isl_vec *aff)
2252 isl_poly_cst *cst;
2253 isl_poly_rec *rec;
2255 if (!poly || !aff)
2256 return -1;
2258 if (poly->var < 0) {
2259 isl_poly_cst *cst;
2261 cst = isl_poly_as_cst(poly);
2262 if (!cst)
2263 return -1;
2264 update_coeff(aff, cst, 0);
2265 return 0;
2268 rec = isl_poly_as_rec(poly);
2269 if (!rec)
2270 return -1;
2271 isl_assert(poly->ctx, rec->n == 2, return -1);
2273 cst = isl_poly_as_cst(rec->p[1]);
2274 if (!cst)
2275 return -1;
2276 update_coeff(aff, cst, 1 + poly->var);
2278 return isl_poly_update_affine(rec->p[0], aff);
2281 __isl_give isl_vec *isl_qpolynomial_extract_affine(
2282 __isl_keep isl_qpolynomial *qp)
2284 isl_vec *aff;
2285 isl_size d;
2287 d = isl_qpolynomial_domain_dim(qp, isl_dim_all);
2288 if (d < 0)
2289 return NULL;
2291 aff = isl_vec_alloc(qp->div->ctx, 2 + d);
2292 if (!aff)
2293 return NULL;
2295 isl_seq_clr(aff->el + 1, 1 + d);
2296 isl_int_set_si(aff->el[0], 1);
2298 if (isl_poly_update_affine(qp->poly, aff) < 0)
2299 goto error;
2301 return aff;
2302 error:
2303 isl_vec_free(aff);
2304 return NULL;
2307 /* Compare two quasi-polynomials.
2309 * Return -1 if "qp1" is "smaller" than "qp2", 1 if "qp1" is "greater"
2310 * than "qp2" and 0 if they are equal.
2312 int isl_qpolynomial_plain_cmp(__isl_keep isl_qpolynomial *qp1,
2313 __isl_keep isl_qpolynomial *qp2)
2315 int cmp;
2317 if (qp1 == qp2)
2318 return 0;
2319 if (!qp1)
2320 return -1;
2321 if (!qp2)
2322 return 1;
2324 cmp = isl_space_cmp(qp1->dim, qp2->dim);
2325 if (cmp != 0)
2326 return cmp;
2328 cmp = isl_local_cmp(qp1->div, qp2->div);
2329 if (cmp != 0)
2330 return cmp;
2332 return isl_poly_plain_cmp(qp1->poly, qp2->poly);
2335 /* Is "qp1" obviously equal to "qp2"?
2337 * NaN is not equal to anything, not even to another NaN.
2339 isl_bool isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
2340 __isl_keep isl_qpolynomial *qp2)
2342 isl_bool equal;
2344 if (!qp1 || !qp2)
2345 return isl_bool_error;
2347 if (isl_qpolynomial_is_nan(qp1) || isl_qpolynomial_is_nan(qp2))
2348 return isl_bool_false;
2350 equal = isl_space_is_equal(qp1->dim, qp2->dim);
2351 if (equal < 0 || !equal)
2352 return equal;
2354 equal = isl_mat_is_equal(qp1->div, qp2->div);
2355 if (equal < 0 || !equal)
2356 return equal;
2358 return isl_poly_is_equal(qp1->poly, qp2->poly);
2361 static isl_stat poly_update_den(__isl_keep isl_poly *poly, isl_int *d)
2363 int i;
2364 isl_bool is_cst;
2365 isl_poly_rec *rec;
2367 is_cst = isl_poly_is_cst(poly);
2368 if (is_cst < 0)
2369 return isl_stat_error;
2370 if (is_cst) {
2371 isl_poly_cst *cst;
2372 cst = isl_poly_as_cst(poly);
2373 if (!cst)
2374 return isl_stat_error;
2375 isl_int_lcm(*d, *d, cst->d);
2376 return isl_stat_ok;
2379 rec = isl_poly_as_rec(poly);
2380 if (!rec)
2381 return isl_stat_error;
2383 for (i = 0; i < rec->n; ++i)
2384 poly_update_den(rec->p[i], d);
2386 return isl_stat_ok;
2389 __isl_give isl_val *isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp)
2391 isl_val *d;
2393 if (!qp)
2394 return NULL;
2395 d = isl_val_one(isl_qpolynomial_get_ctx(qp));
2396 if (!d)
2397 return NULL;
2398 if (poly_update_den(qp->poly, &d->n) < 0)
2399 return isl_val_free(d);
2400 return d;
2403 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
2404 __isl_take isl_space *domain, int pos, int power)
2406 struct isl_ctx *ctx;
2408 if (!domain)
2409 return NULL;
2411 ctx = domain->ctx;
2413 return isl_qpolynomial_alloc(domain, 0,
2414 isl_poly_var_pow(ctx, pos, power));
2417 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(
2418 __isl_take isl_space *domain, enum isl_dim_type type, unsigned pos)
2420 isl_size off;
2422 if (isl_space_check_is_set(domain ) < 0)
2423 goto error;
2424 if (isl_space_check_range(domain, type, pos, 1) < 0)
2425 goto error;
2427 off = isl_space_offset(domain, type);
2428 if (off < 0)
2429 goto error;
2431 return isl_qpolynomial_var_pow_on_domain(domain, off + pos, 1);
2432 error:
2433 isl_space_free(domain);
2434 return NULL;
2437 __isl_give isl_poly *isl_poly_subs(__isl_take isl_poly *poly,
2438 unsigned first, unsigned n, __isl_keep isl_poly **subs)
2440 int i;
2441 isl_bool is_cst;
2442 isl_poly_rec *rec;
2443 isl_poly *base, *res;
2445 is_cst = isl_poly_is_cst(poly);
2446 if (is_cst < 0)
2447 return isl_poly_free(poly);
2448 if (is_cst)
2449 return poly;
2451 if (poly->var < first)
2452 return poly;
2454 rec = isl_poly_as_rec(poly);
2455 if (!rec)
2456 goto error;
2458 isl_assert(poly->ctx, rec->n >= 1, goto error);
2460 if (poly->var >= first + n)
2461 base = isl_poly_var_pow(poly->ctx, poly->var, 1);
2462 else
2463 base = isl_poly_copy(subs[poly->var - first]);
2465 res = isl_poly_subs(isl_poly_copy(rec->p[rec->n - 1]), first, n, subs);
2466 for (i = rec->n - 2; i >= 0; --i) {
2467 isl_poly *t;
2468 t = isl_poly_subs(isl_poly_copy(rec->p[i]), first, n, subs);
2469 res = isl_poly_mul(res, isl_poly_copy(base));
2470 res = isl_poly_sum(res, t);
2473 isl_poly_free(base);
2474 isl_poly_free(poly);
2476 return res;
2477 error:
2478 isl_poly_free(poly);
2479 return NULL;
2482 __isl_give isl_poly *isl_poly_from_affine(isl_ctx *ctx, isl_int *f,
2483 isl_int denom, unsigned len)
2485 int i;
2486 isl_poly *poly;
2488 isl_assert(ctx, len >= 1, return NULL);
2490 poly = isl_poly_rat_cst(ctx, f[0], denom);
2491 for (i = 0; i < len - 1; ++i) {
2492 isl_poly *t;
2493 isl_poly *c;
2495 if (isl_int_is_zero(f[1 + i]))
2496 continue;
2498 c = isl_poly_rat_cst(ctx, f[1 + i], denom);
2499 t = isl_poly_var_pow(ctx, i, 1);
2500 t = isl_poly_mul(c, t);
2501 poly = isl_poly_sum(poly, t);
2504 return poly;
2507 /* Remove common factor of non-constant terms and denominator.
2509 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2511 isl_ctx *ctx = qp->div->ctx;
2512 unsigned total = qp->div->n_col - 2;
2514 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2515 isl_int_gcd(ctx->normalize_gcd,
2516 ctx->normalize_gcd, qp->div->row[div][0]);
2517 if (isl_int_is_one(ctx->normalize_gcd))
2518 return;
2520 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2521 ctx->normalize_gcd, total);
2522 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2523 ctx->normalize_gcd);
2524 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2525 ctx->normalize_gcd);
2528 /* Replace the integer division identified by "div" by the polynomial "s".
2529 * The integer division is assumed not to appear in the definition
2530 * of any other integer divisions.
2532 static __isl_give isl_qpolynomial *substitute_div(
2533 __isl_take isl_qpolynomial *qp, int div, __isl_take isl_poly *s)
2535 int i;
2536 isl_size div_pos;
2537 int *reordering;
2538 isl_ctx *ctx;
2540 if (!qp || !s)
2541 goto error;
2543 qp = isl_qpolynomial_cow(qp);
2544 if (!qp)
2545 goto error;
2547 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
2548 if (div_pos < 0)
2549 goto error;
2550 qp->poly = isl_poly_subs(qp->poly, div_pos + div, 1, &s);
2551 if (!qp->poly)
2552 goto error;
2554 ctx = isl_qpolynomial_get_ctx(qp);
2555 reordering = isl_alloc_array(ctx, int, div_pos + qp->div->n_row);
2556 if (!reordering)
2557 goto error;
2558 for (i = 0; i < div_pos + div; ++i)
2559 reordering[i] = i;
2560 for (i = div_pos + div + 1; i < div_pos + qp->div->n_row; ++i)
2561 reordering[i] = i - 1;
2562 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2563 qp->div = isl_mat_drop_cols(qp->div, 2 + div_pos + div, 1);
2564 qp->poly = reorder(qp->poly, reordering);
2565 free(reordering);
2567 if (!qp->poly || !qp->div)
2568 goto error;
2570 isl_poly_free(s);
2571 return qp;
2572 error:
2573 isl_qpolynomial_free(qp);
2574 isl_poly_free(s);
2575 return NULL;
2578 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2579 * divisions because d is equal to 1 by their definition, i.e., e.
2581 static __isl_give isl_qpolynomial *substitute_non_divs(
2582 __isl_take isl_qpolynomial *qp)
2584 int i, j;
2585 isl_size div_pos;
2586 isl_poly *s;
2588 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
2589 if (div_pos < 0)
2590 return isl_qpolynomial_free(qp);
2592 for (i = 0; qp && i < qp->div->n_row; ++i) {
2593 if (!isl_int_is_one(qp->div->row[i][0]))
2594 continue;
2595 for (j = i + 1; j < qp->div->n_row; ++j) {
2596 if (isl_int_is_zero(qp->div->row[j][2 + div_pos + i]))
2597 continue;
2598 isl_seq_combine(qp->div->row[j] + 1,
2599 qp->div->ctx->one, qp->div->row[j] + 1,
2600 qp->div->row[j][2 + div_pos + i],
2601 qp->div->row[i] + 1, 1 + div_pos + i);
2602 isl_int_set_si(qp->div->row[j][2 + div_pos + i], 0);
2603 normalize_div(qp, j);
2605 s = isl_poly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2606 qp->div->row[i][0], qp->div->n_col - 1);
2607 qp = substitute_div(qp, i, s);
2608 --i;
2611 return qp;
2614 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2615 * with d the denominator. When replacing the coefficient e of x by
2616 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2617 * inside the division, so we need to add floor(e/d) * x outside.
2618 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2619 * to adjust the coefficient of x in each later div that depends on the
2620 * current div "div" and also in the affine expressions in the rows of "mat"
2621 * (if they too depend on "div").
2623 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2624 __isl_keep isl_mat **mat)
2626 int i, j;
2627 isl_int v;
2628 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2630 isl_int_init(v);
2631 for (i = 0; i < 1 + total + div; ++i) {
2632 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2633 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2634 continue;
2635 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2636 isl_int_fdiv_r(qp->div->row[div][1 + i],
2637 qp->div->row[div][1 + i], qp->div->row[div][0]);
2638 *mat = isl_mat_col_addmul(*mat, i, v, 1 + total + div);
2639 for (j = div + 1; j < qp->div->n_row; ++j) {
2640 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2641 continue;
2642 isl_int_addmul(qp->div->row[j][1 + i],
2643 v, qp->div->row[j][2 + total + div]);
2646 isl_int_clear(v);
2649 /* Check if the last non-zero coefficient is bigger that half of the
2650 * denominator. If so, we will invert the div to further reduce the number
2651 * of distinct divs that may appear.
2652 * If the last non-zero coefficient is exactly half the denominator,
2653 * then we continue looking for earlier coefficients that are bigger
2654 * than half the denominator.
2656 static int needs_invert(__isl_keep isl_mat *div, int row)
2658 int i;
2659 int cmp;
2661 for (i = div->n_col - 1; i >= 1; --i) {
2662 if (isl_int_is_zero(div->row[row][i]))
2663 continue;
2664 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2665 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2666 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2667 if (cmp)
2668 return cmp > 0;
2669 if (i == 1)
2670 return 1;
2673 return 0;
2676 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2677 * We only invert the coefficients of e (and the coefficient of q in
2678 * later divs and in the rows of "mat"). After calling this function, the
2679 * coefficients of e should be reduced again.
2681 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2682 __isl_keep isl_mat **mat)
2684 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2686 isl_seq_neg(qp->div->row[div] + 1,
2687 qp->div->row[div] + 1, qp->div->n_col - 1);
2688 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2689 isl_int_add(qp->div->row[div][1],
2690 qp->div->row[div][1], qp->div->row[div][0]);
2691 *mat = isl_mat_col_neg(*mat, 1 + total + div);
2692 isl_mat_col_mul(qp->div, 2 + total + div,
2693 qp->div->ctx->negone, 2 + total + div);
2696 /* Reduce all divs of "qp" to have coefficients
2697 * in the interval [0, d-1], with d the denominator and such that the
2698 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2699 * The modifications to the integer divisions need to be reflected
2700 * in the factors of the polynomial that refer to the original
2701 * integer divisions. To this end, the modifications are collected
2702 * as a set of affine expressions and then plugged into the polynomial.
2704 * After the reduction, some divs may have become redundant or identical,
2705 * so we call substitute_non_divs and sort_divs. If these functions
2706 * eliminate divs or merge two or more divs into one, the coefficients
2707 * of the enclosing divs may have to be reduced again, so we call
2708 * ourselves recursively if the number of divs decreases.
2710 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2712 int i;
2713 isl_ctx *ctx;
2714 isl_mat *mat;
2715 isl_poly **s;
2716 unsigned o_div;
2717 isl_size n_div, total, new_n_div;
2719 total = isl_qpolynomial_domain_dim(qp, isl_dim_all);
2720 n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div);
2721 o_div = isl_qpolynomial_domain_offset(qp, isl_dim_div);
2722 if (total < 0 || n_div < 0)
2723 return isl_qpolynomial_free(qp);
2724 ctx = isl_qpolynomial_get_ctx(qp);
2725 mat = isl_mat_zero(ctx, n_div, 1 + total);
2727 for (i = 0; i < n_div; ++i)
2728 mat = isl_mat_set_element_si(mat, i, o_div + i, 1);
2730 for (i = 0; i < qp->div->n_row; ++i) {
2731 normalize_div(qp, i);
2732 reduce_div(qp, i, &mat);
2733 if (needs_invert(qp->div, i)) {
2734 invert_div(qp, i, &mat);
2735 reduce_div(qp, i, &mat);
2738 if (!mat)
2739 goto error;
2741 s = isl_alloc_array(ctx, struct isl_poly *, n_div);
2742 if (n_div && !s)
2743 goto error;
2744 for (i = 0; i < n_div; ++i)
2745 s[i] = isl_poly_from_affine(ctx, mat->row[i], ctx->one,
2746 1 + total);
2747 qp->poly = isl_poly_subs(qp->poly, o_div - 1, n_div, s);
2748 for (i = 0; i < n_div; ++i)
2749 isl_poly_free(s[i]);
2750 free(s);
2751 if (!qp->poly)
2752 goto error;
2754 isl_mat_free(mat);
2756 qp = substitute_non_divs(qp);
2757 qp = sort_divs(qp);
2758 new_n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div);
2759 if (new_n_div < 0)
2760 return isl_qpolynomial_free(qp);
2761 if (new_n_div < n_div)
2762 return reduce_divs(qp);
2764 return qp;
2765 error:
2766 isl_qpolynomial_free(qp);
2767 isl_mat_free(mat);
2768 return NULL;
2771 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2772 __isl_take isl_space *domain, const isl_int n, const isl_int d)
2774 struct isl_qpolynomial *qp;
2775 isl_poly_cst *cst;
2777 qp = isl_qpolynomial_zero_on_domain(domain);
2778 if (!qp)
2779 return NULL;
2781 cst = isl_poly_as_cst(qp->poly);
2782 isl_int_set(cst->n, n);
2783 isl_int_set(cst->d, d);
2785 return qp;
2788 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2790 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2791 __isl_take isl_space *domain, __isl_take isl_val *val)
2793 isl_qpolynomial *qp;
2794 isl_poly_cst *cst;
2796 qp = isl_qpolynomial_zero_on_domain(domain);
2797 if (!qp || !val)
2798 goto error;
2800 cst = isl_poly_as_cst(qp->poly);
2801 isl_int_set(cst->n, val->n);
2802 isl_int_set(cst->d, val->d);
2804 isl_val_free(val);
2805 return qp;
2806 error:
2807 isl_val_free(val);
2808 isl_qpolynomial_free(qp);
2809 return NULL;
2812 static isl_stat poly_set_active(__isl_keep isl_poly *poly, int *active, int d)
2814 isl_bool is_cst;
2815 isl_poly_rec *rec;
2816 int i;
2818 is_cst = isl_poly_is_cst(poly);
2819 if (is_cst < 0)
2820 return isl_stat_error;
2821 if (is_cst)
2822 return isl_stat_ok;
2824 if (poly->var < d)
2825 active[poly->var] = 1;
2827 rec = isl_poly_as_rec(poly);
2828 for (i = 0; i < rec->n; ++i)
2829 if (poly_set_active(rec->p[i], active, d) < 0)
2830 return isl_stat_error;
2832 return isl_stat_ok;
2835 static isl_stat set_active(__isl_keep isl_qpolynomial *qp, int *active)
2837 int i, j;
2838 isl_size d;
2839 isl_space *space;
2841 space = isl_qpolynomial_peek_domain_space(qp);
2842 d = isl_space_dim(space, isl_dim_all);
2843 if (d < 0 || !active)
2844 return isl_stat_error;
2846 for (i = 0; i < d; ++i)
2847 for (j = 0; j < qp->div->n_row; ++j) {
2848 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2849 continue;
2850 active[i] = 1;
2851 break;
2854 return poly_set_active(qp->poly, active, d);
2857 #undef TYPE
2858 #define TYPE isl_qpolynomial
2859 static
2860 #include "check_type_range_templ.c"
2862 isl_bool isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2863 enum isl_dim_type type, unsigned first, unsigned n)
2865 int i;
2866 int *active = NULL;
2867 isl_bool involves = isl_bool_false;
2868 isl_size offset;
2869 isl_size d;
2870 isl_space *space;
2872 if (!qp)
2873 return isl_bool_error;
2874 if (n == 0)
2875 return isl_bool_false;
2877 if (isl_qpolynomial_check_range(qp, type, first, n) < 0)
2878 return isl_bool_error;
2879 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2880 type == isl_dim_in, return isl_bool_error);
2882 space = isl_qpolynomial_peek_domain_space(qp);
2883 d = isl_space_dim(space, isl_dim_all);
2884 if (d < 0)
2885 return isl_bool_error;
2886 active = isl_calloc_array(qp->dim->ctx, int, d);
2887 if (set_active(qp, active) < 0)
2888 goto error;
2890 offset = isl_qpolynomial_domain_var_offset(qp, domain_type(type));
2891 if (offset < 0)
2892 goto error;
2893 first += offset;
2894 for (i = 0; i < n; ++i)
2895 if (active[first + i]) {
2896 involves = isl_bool_true;
2897 break;
2900 free(active);
2902 return involves;
2903 error:
2904 free(active);
2905 return isl_bool_error;
2908 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2909 * of the divs that do appear in the quasi-polynomial.
2911 static __isl_give isl_qpolynomial *remove_redundant_divs(
2912 __isl_take isl_qpolynomial *qp)
2914 int i, j;
2915 isl_size div_pos;
2916 int len;
2917 int skip;
2918 int *active = NULL;
2919 int *reordering = NULL;
2920 int redundant = 0;
2921 int n_div;
2922 isl_ctx *ctx;
2924 if (!qp)
2925 return NULL;
2926 if (qp->div->n_row == 0)
2927 return qp;
2929 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
2930 if (div_pos < 0)
2931 return isl_qpolynomial_free(qp);
2932 len = qp->div->n_col - 2;
2933 ctx = isl_qpolynomial_get_ctx(qp);
2934 active = isl_calloc_array(ctx, int, len);
2935 if (!active)
2936 goto error;
2938 if (poly_set_active(qp->poly, active, len) < 0)
2939 goto error;
2941 for (i = qp->div->n_row - 1; i >= 0; --i) {
2942 if (!active[div_pos + i]) {
2943 redundant = 1;
2944 continue;
2946 for (j = 0; j < i; ++j) {
2947 if (isl_int_is_zero(qp->div->row[i][2 + div_pos + j]))
2948 continue;
2949 active[div_pos + j] = 1;
2950 break;
2954 if (!redundant) {
2955 free(active);
2956 return qp;
2959 reordering = isl_alloc_array(qp->div->ctx, int, len);
2960 if (!reordering)
2961 goto error;
2963 for (i = 0; i < div_pos; ++i)
2964 reordering[i] = i;
2966 skip = 0;
2967 n_div = qp->div->n_row;
2968 for (i = 0; i < n_div; ++i) {
2969 if (!active[div_pos + i]) {
2970 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2971 qp->div = isl_mat_drop_cols(qp->div,
2972 2 + div_pos + i - skip, 1);
2973 skip++;
2975 reordering[div_pos + i] = div_pos + i - skip;
2978 qp->poly = reorder(qp->poly, reordering);
2980 if (!qp->poly || !qp->div)
2981 goto error;
2983 free(active);
2984 free(reordering);
2986 return qp;
2987 error:
2988 free(active);
2989 free(reordering);
2990 isl_qpolynomial_free(qp);
2991 return NULL;
2994 __isl_give isl_poly *isl_poly_drop(__isl_take isl_poly *poly,
2995 unsigned first, unsigned n)
2997 int i;
2998 isl_poly_rec *rec;
3000 if (!poly)
3001 return NULL;
3002 if (n == 0 || poly->var < 0 || poly->var < first)
3003 return poly;
3004 if (poly->var < first + n) {
3005 poly = replace_by_constant_term(poly);
3006 return isl_poly_drop(poly, first, n);
3008 poly = isl_poly_cow(poly);
3009 if (!poly)
3010 return NULL;
3011 poly->var -= n;
3012 rec = isl_poly_as_rec(poly);
3013 if (!rec)
3014 goto error;
3016 for (i = 0; i < rec->n; ++i) {
3017 rec->p[i] = isl_poly_drop(rec->p[i], first, n);
3018 if (!rec->p[i])
3019 goto error;
3022 return poly;
3023 error:
3024 isl_poly_free(poly);
3025 return NULL;
3028 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
3029 __isl_take isl_qpolynomial *qp,
3030 enum isl_dim_type type, unsigned pos, const char *s)
3032 isl_space *space;
3034 if (!qp)
3035 return NULL;
3036 if (type == isl_dim_out)
3037 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
3038 "cannot set name of output/set dimension",
3039 return isl_qpolynomial_free(qp));
3040 type = domain_type(type);
3041 space = isl_qpolynomial_take_domain_space(qp);
3042 space = isl_space_set_dim_name(space, type, pos, s);
3043 qp = isl_qpolynomial_restore_domain_space(qp, space);
3044 return qp;
3047 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
3048 __isl_take isl_qpolynomial *qp,
3049 enum isl_dim_type type, unsigned first, unsigned n)
3051 isl_space *space;
3052 isl_size offset;
3054 if (!qp)
3055 return NULL;
3056 if (type == isl_dim_out)
3057 isl_die(qp->dim->ctx, isl_error_invalid,
3058 "cannot drop output/set dimension",
3059 goto error);
3060 if (isl_qpolynomial_check_range(qp, type, first, n) < 0)
3061 return isl_qpolynomial_free(qp);
3062 type = domain_type(type);
3063 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3064 return qp;
3067 isl_assert(qp->dim->ctx, type == isl_dim_param ||
3068 type == isl_dim_set, goto error);
3070 space = isl_qpolynomial_take_domain_space(qp);
3071 space = isl_space_drop_dims(space, type, first, n);
3072 qp = isl_qpolynomial_restore_domain_space(qp, space);
3074 qp = isl_qpolynomial_cow(qp);
3075 if (!qp)
3076 return NULL;
3078 offset = isl_qpolynomial_domain_var_offset(qp, type);
3079 if (offset < 0)
3080 goto error;
3081 first += offset;
3083 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
3084 if (!qp->div)
3085 goto error;
3087 qp->poly = isl_poly_drop(qp->poly, first, n);
3088 if (!qp->poly)
3089 goto error;
3091 return qp;
3092 error:
3093 isl_qpolynomial_free(qp);
3094 return NULL;
3097 /* Project the domain of the quasi-polynomial onto its parameter space.
3098 * The quasi-polynomial may not involve any of the domain dimensions.
3100 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
3101 __isl_take isl_qpolynomial *qp)
3103 isl_space *space;
3104 isl_size n;
3105 isl_bool involves;
3107 n = isl_qpolynomial_dim(qp, isl_dim_in);
3108 if (n < 0)
3109 return isl_qpolynomial_free(qp);
3110 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
3111 if (involves < 0)
3112 return isl_qpolynomial_free(qp);
3113 if (involves)
3114 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
3115 "polynomial involves some of the domain dimensions",
3116 return isl_qpolynomial_free(qp));
3117 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
3118 space = isl_qpolynomial_get_domain_space(qp);
3119 space = isl_space_params(space);
3120 qp = isl_qpolynomial_reset_domain_space(qp, space);
3121 return qp;
3124 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
3125 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
3127 int i, j, k;
3128 isl_int denom;
3129 unsigned total;
3130 unsigned n_div;
3131 isl_poly *poly;
3133 if (!eq)
3134 goto error;
3135 if (eq->n_eq == 0) {
3136 isl_basic_set_free(eq);
3137 return qp;
3140 qp = isl_qpolynomial_cow(qp);
3141 if (!qp)
3142 goto error;
3143 qp->div = isl_mat_cow(qp->div);
3144 if (!qp->div)
3145 goto error;
3147 total = isl_basic_set_offset(eq, isl_dim_div);
3148 n_div = eq->n_div;
3149 isl_int_init(denom);
3150 for (i = 0; i < eq->n_eq; ++i) {
3151 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
3152 if (j < 0 || j == 0 || j >= total)
3153 continue;
3155 for (k = 0; k < qp->div->n_row; ++k) {
3156 if (isl_int_is_zero(qp->div->row[k][1 + j]))
3157 continue;
3158 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
3159 &qp->div->row[k][0]);
3160 normalize_div(qp, k);
3163 if (isl_int_is_pos(eq->eq[i][j]))
3164 isl_seq_neg(eq->eq[i], eq->eq[i], total);
3165 isl_int_abs(denom, eq->eq[i][j]);
3166 isl_int_set_si(eq->eq[i][j], 0);
3168 poly = isl_poly_from_affine(qp->dim->ctx,
3169 eq->eq[i], denom, total);
3170 qp->poly = isl_poly_subs(qp->poly, j - 1, 1, &poly);
3171 isl_poly_free(poly);
3173 isl_int_clear(denom);
3175 if (!qp->poly)
3176 goto error;
3178 isl_basic_set_free(eq);
3180 qp = substitute_non_divs(qp);
3181 qp = sort_divs(qp);
3183 return qp;
3184 error:
3185 isl_basic_set_free(eq);
3186 isl_qpolynomial_free(qp);
3187 return NULL;
3190 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
3192 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
3193 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
3195 if (!qp || !eq)
3196 goto error;
3197 if (qp->div->n_row > 0)
3198 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
3199 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
3200 error:
3201 isl_basic_set_free(eq);
3202 isl_qpolynomial_free(qp);
3203 return NULL;
3206 /* Look for equalities among the variables shared by context and qp
3207 * and the integer divisions of qp, if any.
3208 * The equalities are then used to eliminate variables and/or integer
3209 * divisions from qp.
3211 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
3212 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
3214 isl_local_space *ls;
3215 isl_basic_set *aff;
3217 ls = isl_qpolynomial_get_domain_local_space(qp);
3218 context = isl_local_space_lift_set(ls, context);
3220 aff = isl_set_affine_hull(context);
3221 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
3224 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
3225 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
3227 isl_space *space = isl_qpolynomial_get_domain_space(qp);
3228 isl_set *dom_context = isl_set_universe(space);
3229 dom_context = isl_set_intersect_params(dom_context, context);
3230 return isl_qpolynomial_gist(qp, dom_context);
3233 /* Return a zero isl_qpolynomial in the given space.
3235 * This is a helper function for isl_pw_*_as_* that ensures a uniform
3236 * interface over all piecewise types.
3238 static __isl_give isl_qpolynomial *isl_qpolynomial_zero_in_space(
3239 __isl_take isl_space *space)
3241 return isl_qpolynomial_zero_on_domain(isl_space_domain(space));
3244 #define isl_qpolynomial_involves_nan isl_qpolynomial_is_nan
3246 #undef PW
3247 #define PW isl_pw_qpolynomial
3248 #undef BASE
3249 #define BASE qpolynomial
3250 #undef EL_IS_ZERO
3251 #define EL_IS_ZERO is_zero
3252 #undef ZERO
3253 #define ZERO zero
3254 #undef IS_ZERO
3255 #define IS_ZERO is_zero
3256 #undef FIELD
3257 #define FIELD qp
3258 #undef DEFAULT_IS_ZERO
3259 #define DEFAULT_IS_ZERO 1
3261 #include <isl_pw_templ.c>
3262 #include <isl_pw_un_op_templ.c>
3263 #include <isl_pw_add_disjoint_templ.c>
3264 #include <isl_pw_domain_reverse_templ.c>
3265 #include <isl_pw_eval.c>
3266 #include <isl_pw_fix_templ.c>
3267 #include <isl_pw_from_range_templ.c>
3268 #include <isl_pw_insert_dims_templ.c>
3269 #include <isl_pw_lift_templ.c>
3270 #include <isl_pw_morph_templ.c>
3271 #include <isl_pw_move_dims_templ.c>
3272 #include <isl_pw_neg_templ.c>
3273 #include <isl_pw_opt_templ.c>
3274 #include <isl_pw_split_dims_templ.c>
3275 #include <isl_pw_sub_templ.c>
3277 #undef BASE
3278 #define BASE pw_qpolynomial
3280 #include <isl_union_single.c>
3281 #include <isl_union_domain_reverse_templ.c>
3282 #include <isl_union_eval.c>
3283 #include <isl_union_neg.c>
3284 #include <isl_union_sub_templ.c>
3286 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
3288 if (!pwqp)
3289 return -1;
3291 if (pwqp->n != -1)
3292 return 0;
3294 if (!isl_set_plain_is_universe(pwqp->p[0].set))
3295 return 0;
3297 return isl_qpolynomial_is_one(pwqp->p[0].qp);
3300 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
3301 __isl_take isl_pw_qpolynomial *pwqp1,
3302 __isl_take isl_pw_qpolynomial *pwqp2)
3304 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
3307 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
3308 __isl_take isl_pw_qpolynomial *pwqp1,
3309 __isl_take isl_pw_qpolynomial *pwqp2)
3311 int i, j, n;
3312 struct isl_pw_qpolynomial *res;
3314 if (!pwqp1 || !pwqp2)
3315 goto error;
3317 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
3318 goto error);
3320 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
3321 isl_pw_qpolynomial_free(pwqp2);
3322 return pwqp1;
3325 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
3326 isl_pw_qpolynomial_free(pwqp1);
3327 return pwqp2;
3330 if (isl_pw_qpolynomial_is_one(pwqp1)) {
3331 isl_pw_qpolynomial_free(pwqp1);
3332 return pwqp2;
3335 if (isl_pw_qpolynomial_is_one(pwqp2)) {
3336 isl_pw_qpolynomial_free(pwqp2);
3337 return pwqp1;
3340 n = pwqp1->n * pwqp2->n;
3341 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
3343 for (i = 0; i < pwqp1->n; ++i) {
3344 for (j = 0; j < pwqp2->n; ++j) {
3345 struct isl_set *common;
3346 struct isl_qpolynomial *prod;
3347 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
3348 isl_set_copy(pwqp2->p[j].set));
3349 if (isl_set_plain_is_empty(common)) {
3350 isl_set_free(common);
3351 continue;
3354 prod = isl_qpolynomial_mul(
3355 isl_qpolynomial_copy(pwqp1->p[i].qp),
3356 isl_qpolynomial_copy(pwqp2->p[j].qp));
3358 res = isl_pw_qpolynomial_add_piece(res, common, prod);
3362 isl_pw_qpolynomial_free(pwqp1);
3363 isl_pw_qpolynomial_free(pwqp2);
3365 return res;
3366 error:
3367 isl_pw_qpolynomial_free(pwqp1);
3368 isl_pw_qpolynomial_free(pwqp2);
3369 return NULL;
3372 __isl_give isl_val *isl_poly_eval(__isl_take isl_poly *poly,
3373 __isl_take isl_vec *vec)
3375 int i;
3376 isl_bool is_cst;
3377 isl_poly_rec *rec;
3378 isl_val *res;
3379 isl_val *base;
3381 is_cst = isl_poly_is_cst(poly);
3382 if (is_cst < 0)
3383 goto error;
3384 if (is_cst) {
3385 isl_vec_free(vec);
3386 res = isl_poly_get_constant_val(poly);
3387 isl_poly_free(poly);
3388 return res;
3391 rec = isl_poly_as_rec(poly);
3392 if (!rec || !vec)
3393 goto error;
3395 isl_assert(poly->ctx, rec->n >= 1, goto error);
3397 base = isl_val_rat_from_isl_int(poly->ctx,
3398 vec->el[1 + poly->var], vec->el[0]);
3400 res = isl_poly_eval(isl_poly_copy(rec->p[rec->n - 1]),
3401 isl_vec_copy(vec));
3403 for (i = rec->n - 2; i >= 0; --i) {
3404 res = isl_val_mul(res, isl_val_copy(base));
3405 res = isl_val_add(res, isl_poly_eval(isl_poly_copy(rec->p[i]),
3406 isl_vec_copy(vec)));
3409 isl_val_free(base);
3410 isl_poly_free(poly);
3411 isl_vec_free(vec);
3412 return res;
3413 error:
3414 isl_poly_free(poly);
3415 isl_vec_free(vec);
3416 return NULL;
3419 /* Evaluate "qp" in the void point "pnt".
3420 * In particular, return the value NaN.
3422 static __isl_give isl_val *eval_void(__isl_take isl_qpolynomial *qp,
3423 __isl_take isl_point *pnt)
3425 isl_ctx *ctx;
3427 ctx = isl_point_get_ctx(pnt);
3428 isl_qpolynomial_free(qp);
3429 isl_point_free(pnt);
3430 return isl_val_nan(ctx);
3433 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
3434 __isl_take isl_point *pnt)
3436 isl_bool is_void;
3437 isl_vec *ext;
3438 isl_val *v;
3440 if (!qp || !pnt)
3441 goto error;
3442 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
3443 is_void = isl_point_is_void(pnt);
3444 if (is_void < 0)
3445 goto error;
3446 if (is_void)
3447 return eval_void(qp, pnt);
3449 ext = isl_local_extend_point_vec(qp->div, isl_vec_copy(pnt->vec));
3451 v = isl_poly_eval(isl_qpolynomial_get_poly(qp), ext);
3453 isl_qpolynomial_free(qp);
3454 isl_point_free(pnt);
3456 return v;
3457 error:
3458 isl_qpolynomial_free(qp);
3459 isl_point_free(pnt);
3460 return NULL;
3463 int isl_poly_cmp(__isl_keep isl_poly_cst *cst1, __isl_keep isl_poly_cst *cst2)
3465 int cmp;
3466 isl_int t;
3467 isl_int_init(t);
3468 isl_int_mul(t, cst1->n, cst2->d);
3469 isl_int_submul(t, cst2->n, cst1->d);
3470 cmp = isl_int_sgn(t);
3471 isl_int_clear(t);
3472 return cmp;
3475 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
3476 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
3477 unsigned first, unsigned n)
3479 unsigned total;
3480 unsigned g_pos;
3481 int *exp;
3482 isl_space *space;
3484 if (!qp)
3485 return NULL;
3486 if (type == isl_dim_out)
3487 isl_die(qp->div->ctx, isl_error_invalid,
3488 "cannot insert output/set dimensions",
3489 goto error);
3490 if (isl_qpolynomial_check_range(qp, type, first, 0) < 0)
3491 return isl_qpolynomial_free(qp);
3492 type = domain_type(type);
3493 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3494 return qp;
3496 qp = isl_qpolynomial_cow(qp);
3497 if (!qp)
3498 return NULL;
3500 g_pos = pos(qp->dim, type) + first;
3502 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3503 if (!qp->div)
3504 goto error;
3506 total = qp->div->n_col - 2;
3507 if (total > g_pos) {
3508 int i;
3509 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3510 if (!exp)
3511 goto error;
3512 for (i = 0; i < total - g_pos; ++i)
3513 exp[i] = i + n;
3514 qp->poly = expand(qp->poly, exp, g_pos);
3515 free(exp);
3516 if (!qp->poly)
3517 goto error;
3520 space = isl_qpolynomial_take_domain_space(qp);
3521 space = isl_space_insert_dims(space, type, first, n);
3522 qp = isl_qpolynomial_restore_domain_space(qp, space);
3524 return qp;
3525 error:
3526 isl_qpolynomial_free(qp);
3527 return NULL;
3530 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3531 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3533 isl_size pos;
3535 pos = isl_qpolynomial_dim(qp, type);
3536 if (pos < 0)
3537 return isl_qpolynomial_free(qp);
3539 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3542 static int *reordering_move(isl_ctx *ctx,
3543 unsigned len, unsigned dst, unsigned src, unsigned n)
3545 int i;
3546 int *reordering;
3548 reordering = isl_alloc_array(ctx, int, len);
3549 if (!reordering)
3550 return NULL;
3552 if (dst <= src) {
3553 for (i = 0; i < dst; ++i)
3554 reordering[i] = i;
3555 for (i = 0; i < n; ++i)
3556 reordering[src + i] = dst + i;
3557 for (i = 0; i < src - dst; ++i)
3558 reordering[dst + i] = dst + n + i;
3559 for (i = 0; i < len - src - n; ++i)
3560 reordering[src + n + i] = src + n + i;
3561 } else {
3562 for (i = 0; i < src; ++i)
3563 reordering[i] = i;
3564 for (i = 0; i < n; ++i)
3565 reordering[src + i] = dst + i;
3566 for (i = 0; i < dst - src; ++i)
3567 reordering[src + n + i] = src + i;
3568 for (i = 0; i < len - dst - n; ++i)
3569 reordering[dst + n + i] = dst + n + i;
3572 return reordering;
3575 /* Move the "n" variables starting at "src_pos" of "local" to "dst_pos".
3576 * Only modify the polynomial expression and the local variables of "qp".
3577 * The caller is responsible for modifying the space accordingly.
3579 static __isl_give isl_qpolynomial *local_poly_move_dims(
3580 __isl_take isl_qpolynomial *qp,
3581 unsigned dst_pos, unsigned src_pos, unsigned n)
3583 isl_ctx *ctx;
3584 isl_size total;
3585 int *reordering;
3586 isl_local *local;
3587 isl_poly *poly;
3589 local = isl_qpolynomial_take_local(qp);
3590 local = isl_local_move_vars(local, dst_pos, src_pos, n);
3591 qp = isl_qpolynomial_restore_local(qp, local);
3592 qp = sort_divs(qp);
3594 total = isl_qpolynomial_domain_dim(qp, isl_dim_all);
3595 if (total < 0)
3596 return isl_qpolynomial_free(qp);
3597 ctx = isl_qpolynomial_get_ctx(qp);
3598 reordering = reordering_move(ctx, total, dst_pos, src_pos, n);
3599 if (!reordering)
3600 return isl_qpolynomial_free(qp);
3602 poly = isl_qpolynomial_take_poly(qp);
3603 poly = reorder(poly, reordering);
3604 qp = isl_qpolynomial_restore_poly(qp, poly);
3605 free(reordering);
3607 return qp;
3610 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3611 __isl_take isl_qpolynomial *qp,
3612 enum isl_dim_type dst_type, unsigned dst_pos,
3613 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3615 isl_ctx *ctx;
3616 unsigned g_dst_pos;
3617 unsigned g_src_pos;
3618 isl_size src_off, dst_off;
3619 isl_space *space;
3621 if (!qp)
3622 return NULL;
3624 ctx = isl_qpolynomial_get_ctx(qp);
3625 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3626 isl_die(ctx, isl_error_invalid,
3627 "cannot move output/set dimension",
3628 return isl_qpolynomial_free(qp));
3629 if (src_type == isl_dim_div || dst_type == isl_dim_div)
3630 isl_die(ctx, isl_error_invalid, "cannot move local variables",
3631 return isl_qpolynomial_free(qp));
3632 if (isl_qpolynomial_check_range(qp, src_type, src_pos, n) < 0)
3633 return isl_qpolynomial_free(qp);
3634 if (dst_type == isl_dim_in)
3635 dst_type = isl_dim_set;
3636 if (src_type == isl_dim_in)
3637 src_type = isl_dim_set;
3639 if (n == 0 &&
3640 !isl_space_is_named_or_nested(qp->dim, src_type) &&
3641 !isl_space_is_named_or_nested(qp->dim, dst_type))
3642 return qp;
3644 src_off = isl_qpolynomial_domain_var_offset(qp, src_type);
3645 dst_off = isl_qpolynomial_domain_var_offset(qp, dst_type);
3646 if (src_off < 0 || dst_off < 0)
3647 return isl_qpolynomial_free(qp);
3649 g_dst_pos = dst_off + dst_pos;
3650 g_src_pos = src_off + src_pos;
3651 if (dst_type > src_type)
3652 g_dst_pos -= n;
3654 qp = local_poly_move_dims(qp, g_dst_pos, g_src_pos, n);
3656 space = isl_qpolynomial_take_domain_space(qp);
3657 space = isl_space_move_dims(space, dst_type, dst_pos,
3658 src_type, src_pos, n);
3659 qp = isl_qpolynomial_restore_domain_space(qp, space);
3661 return qp;
3664 /* Given a quasi-polynomial on a domain (A -> B),
3665 * interchange A and B in the wrapped domain
3666 * to obtain a quasi-polynomial on the domain (B -> A).
3668 __isl_give isl_qpolynomial *isl_qpolynomial_domain_reverse(
3669 __isl_take isl_qpolynomial *qp)
3671 isl_space *space;
3672 isl_size n_in, n_out, offset;
3674 space = isl_qpolynomial_peek_domain_space(qp);
3675 offset = isl_space_offset(space, isl_dim_set);
3676 n_in = isl_space_wrapped_dim(space, isl_dim_set, isl_dim_in);
3677 n_out = isl_space_wrapped_dim(space, isl_dim_set, isl_dim_out);
3678 if (offset < 0 || n_in < 0 || n_out < 0)
3679 return isl_qpolynomial_free(qp);
3681 qp = local_poly_move_dims(qp, offset, offset + n_in, n_out);
3683 space = isl_qpolynomial_take_domain_space(qp);
3684 space = isl_space_wrapped_reverse(space);
3685 qp = isl_qpolynomial_restore_domain_space(qp, space);
3687 return qp;
3690 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(
3691 __isl_take isl_space *space, isl_int *f, isl_int denom)
3693 isl_size d;
3694 isl_poly *poly;
3696 space = isl_space_domain(space);
3697 if (!space)
3698 return NULL;
3700 d = isl_space_dim(space, isl_dim_all);
3701 poly = d < 0 ? NULL : isl_poly_from_affine(space->ctx, f, denom, 1 + d);
3703 return isl_qpolynomial_alloc(space, 0, poly);
3706 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3708 isl_ctx *ctx;
3709 isl_poly *poly;
3710 isl_qpolynomial *qp;
3712 if (!aff)
3713 return NULL;
3715 ctx = isl_aff_get_ctx(aff);
3716 poly = isl_poly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3717 aff->v->size - 1);
3719 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3720 aff->ls->div->n_row, poly);
3721 if (!qp)
3722 goto error;
3724 isl_mat_free(qp->div);
3725 qp->div = isl_mat_copy(aff->ls->div);
3726 qp->div = isl_mat_cow(qp->div);
3727 if (!qp->div)
3728 goto error;
3730 isl_aff_free(aff);
3731 qp = reduce_divs(qp);
3732 qp = remove_redundant_divs(qp);
3733 return qp;
3734 error:
3735 isl_aff_free(aff);
3736 return isl_qpolynomial_free(qp);
3739 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3740 __isl_take isl_pw_aff *pwaff)
3742 int i;
3743 isl_pw_qpolynomial *pwqp;
3745 if (!pwaff)
3746 return NULL;
3748 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3749 pwaff->n);
3751 for (i = 0; i < pwaff->n; ++i) {
3752 isl_set *dom;
3753 isl_qpolynomial *qp;
3755 dom = isl_set_copy(pwaff->p[i].set);
3756 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3757 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3760 isl_pw_aff_free(pwaff);
3761 return pwqp;
3764 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3765 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3767 isl_aff *aff;
3769 aff = isl_constraint_get_bound(c, type, pos);
3770 isl_constraint_free(c);
3771 return isl_qpolynomial_from_aff(aff);
3774 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3775 * in "qp" by subs[i].
3777 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3778 __isl_take isl_qpolynomial *qp,
3779 enum isl_dim_type type, unsigned first, unsigned n,
3780 __isl_keep isl_qpolynomial **subs)
3782 int i;
3783 isl_poly *poly;
3784 isl_poly **polys;
3786 if (n == 0)
3787 return qp;
3789 if (!qp)
3790 return NULL;
3792 if (type == isl_dim_out)
3793 isl_die(qp->dim->ctx, isl_error_invalid,
3794 "cannot substitute output/set dimension",
3795 goto error);
3796 if (isl_qpolynomial_check_range(qp, type, first, n) < 0)
3797 return isl_qpolynomial_free(qp);
3798 type = domain_type(type);
3800 for (i = 0; i < n; ++i)
3801 if (!subs[i])
3802 goto error;
3804 for (i = 0; i < n; ++i)
3805 if (isl_qpolynomial_check_equal_space(qp, subs[i]) < 0)
3806 goto error;
3808 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3809 for (i = 0; i < n; ++i)
3810 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3812 first += pos(qp->dim, type);
3814 polys = isl_alloc_array(qp->dim->ctx, struct isl_poly *, n);
3815 if (!polys)
3816 goto error;
3817 for (i = 0; i < n; ++i)
3818 polys[i] = subs[i]->poly;
3820 poly = isl_qpolynomial_take_poly(qp);
3821 poly = isl_poly_subs(poly, first, n, polys);
3822 qp = isl_qpolynomial_restore_poly(qp, poly);
3824 free(polys);
3826 return qp;
3827 error:
3828 isl_qpolynomial_free(qp);
3829 return NULL;
3832 /* Extend "bset" with extra set dimensions for each integer division
3833 * in "qp" and then call "fn" with the extended bset and the polynomial
3834 * that results from replacing each of the integer divisions by the
3835 * corresponding extra set dimension.
3837 isl_stat isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3838 __isl_keep isl_basic_set *bset,
3839 isl_stat (*fn)(__isl_take isl_basic_set *bset,
3840 __isl_take isl_qpolynomial *poly, void *user), void *user)
3842 isl_space *space;
3843 isl_local_space *ls;
3844 isl_poly *poly;
3845 isl_qpolynomial *polynomial;
3847 if (!qp || !bset)
3848 return isl_stat_error;
3849 if (qp->div->n_row == 0)
3850 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3851 user);
3853 space = isl_space_copy(qp->dim);
3854 space = isl_space_add_dims(space, isl_dim_set, qp->div->n_row);
3855 poly = isl_qpolynomial_get_poly(qp);
3856 polynomial = isl_qpolynomial_alloc(space, 0, poly);
3857 bset = isl_basic_set_copy(bset);
3858 ls = isl_qpolynomial_get_domain_local_space(qp);
3859 bset = isl_local_space_lift_basic_set(ls, bset);
3861 return fn(bset, polynomial, user);
3864 /* Return total degree in variables first (inclusive) up to last (exclusive).
3866 int isl_poly_degree(__isl_keep isl_poly *poly, int first, int last)
3868 int deg = -1;
3869 int i;
3870 isl_bool is_zero, is_cst;
3871 isl_poly_rec *rec;
3873 is_zero = isl_poly_is_zero(poly);
3874 if (is_zero < 0)
3875 return -2;
3876 if (is_zero)
3877 return -1;
3878 is_cst = isl_poly_is_cst(poly);
3879 if (is_cst < 0)
3880 return -2;
3881 if (is_cst || poly->var < first)
3882 return 0;
3884 rec = isl_poly_as_rec(poly);
3885 if (!rec)
3886 return -2;
3888 for (i = 0; i < rec->n; ++i) {
3889 int d;
3891 is_zero = isl_poly_is_zero(rec->p[i]);
3892 if (is_zero < 0)
3893 return -2;
3894 if (is_zero)
3895 continue;
3896 d = isl_poly_degree(rec->p[i], first, last);
3897 if (poly->var < last)
3898 d += i;
3899 if (d > deg)
3900 deg = d;
3903 return deg;
3906 /* Return total degree in set variables.
3908 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3910 isl_size ovar;
3911 isl_size nvar;
3913 if (!poly)
3914 return -2;
3916 ovar = isl_space_offset(poly->dim, isl_dim_set);
3917 nvar = isl_space_dim(poly->dim, isl_dim_set);
3918 if (ovar < 0 || nvar < 0)
3919 return -2;
3920 return isl_poly_degree(poly->poly, ovar, ovar + nvar);
3923 __isl_give isl_poly *isl_poly_coeff(__isl_keep isl_poly *poly,
3924 unsigned pos, int deg)
3926 int i;
3927 isl_bool is_cst;
3928 isl_poly_rec *rec;
3930 is_cst = isl_poly_is_cst(poly);
3931 if (is_cst < 0)
3932 return NULL;
3933 if (is_cst || poly->var < pos) {
3934 if (deg == 0)
3935 return isl_poly_copy(poly);
3936 else
3937 return isl_poly_zero(poly->ctx);
3940 rec = isl_poly_as_rec(poly);
3941 if (!rec)
3942 return NULL;
3944 if (poly->var == pos) {
3945 if (deg < rec->n)
3946 return isl_poly_copy(rec->p[deg]);
3947 else
3948 return isl_poly_zero(poly->ctx);
3951 poly = isl_poly_copy(poly);
3952 poly = isl_poly_cow(poly);
3953 rec = isl_poly_as_rec(poly);
3954 if (!rec)
3955 goto error;
3957 for (i = 0; i < rec->n; ++i) {
3958 isl_poly *t;
3959 t = isl_poly_coeff(rec->p[i], pos, deg);
3960 if (!t)
3961 goto error;
3962 isl_poly_free(rec->p[i]);
3963 rec->p[i] = t;
3966 return poly;
3967 error:
3968 isl_poly_free(poly);
3969 return NULL;
3972 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3974 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3975 __isl_keep isl_qpolynomial *qp,
3976 enum isl_dim_type type, unsigned t_pos, int deg)
3978 unsigned g_pos;
3979 isl_poly *poly;
3980 isl_qpolynomial *c;
3982 if (!qp)
3983 return NULL;
3985 if (type == isl_dim_out)
3986 isl_die(qp->div->ctx, isl_error_invalid,
3987 "output/set dimension does not have a coefficient",
3988 return NULL);
3989 if (isl_qpolynomial_check_range(qp, type, t_pos, 1) < 0)
3990 return NULL;
3991 type = domain_type(type);
3993 g_pos = pos(qp->dim, type) + t_pos;
3994 poly = isl_poly_coeff(qp->poly, g_pos, deg);
3996 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim),
3997 qp->div->n_row, poly);
3998 if (!c)
3999 return NULL;
4000 isl_mat_free(c->div);
4001 c->div = isl_qpolynomial_get_local(qp);
4002 if (!c->div)
4003 goto error;
4004 return c;
4005 error:
4006 isl_qpolynomial_free(c);
4007 return NULL;
4010 /* Homogenize the polynomial in the variables first (inclusive) up to
4011 * last (exclusive) by inserting powers of variable first.
4012 * Variable first is assumed not to appear in the input.
4014 __isl_give isl_poly *isl_poly_homogenize(__isl_take isl_poly *poly, int deg,
4015 int target, int first, int last)
4017 int i;
4018 isl_bool is_zero, is_cst;
4019 isl_poly_rec *rec;
4021 is_zero = isl_poly_is_zero(poly);
4022 if (is_zero < 0)
4023 return isl_poly_free(poly);
4024 if (is_zero)
4025 return poly;
4026 if (deg == target)
4027 return poly;
4028 is_cst = isl_poly_is_cst(poly);
4029 if (is_cst < 0)
4030 return isl_poly_free(poly);
4031 if (is_cst || poly->var < first) {
4032 isl_poly *hom;
4034 hom = isl_poly_var_pow(poly->ctx, first, target - deg);
4035 if (!hom)
4036 goto error;
4037 rec = isl_poly_as_rec(hom);
4038 rec->p[target - deg] = isl_poly_mul(rec->p[target - deg], poly);
4040 return hom;
4043 poly = isl_poly_cow(poly);
4044 rec = isl_poly_as_rec(poly);
4045 if (!rec)
4046 goto error;
4048 for (i = 0; i < rec->n; ++i) {
4049 is_zero = isl_poly_is_zero(rec->p[i]);
4050 if (is_zero < 0)
4051 return isl_poly_free(poly);
4052 if (is_zero)
4053 continue;
4054 rec->p[i] = isl_poly_homogenize(rec->p[i],
4055 poly->var < last ? deg + i : i, target,
4056 first, last);
4057 if (!rec->p[i])
4058 goto error;
4061 return poly;
4062 error:
4063 isl_poly_free(poly);
4064 return NULL;
4067 /* Homogenize the polynomial in the set variables by introducing
4068 * powers of an extra set variable at position 0.
4070 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
4071 __isl_take isl_qpolynomial *poly)
4073 isl_size ovar;
4074 isl_size nvar;
4075 int deg = isl_qpolynomial_degree(poly);
4077 if (deg < -1)
4078 goto error;
4080 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
4081 poly = isl_qpolynomial_cow(poly);
4082 if (!poly)
4083 goto error;
4085 ovar = isl_space_offset(poly->dim, isl_dim_set);
4086 nvar = isl_space_dim(poly->dim, isl_dim_set);
4087 if (ovar < 0 || nvar < 0)
4088 return isl_qpolynomial_free(poly);
4089 poly->poly = isl_poly_homogenize(poly->poly, 0, deg, ovar, ovar + nvar);
4090 if (!poly->poly)
4091 goto error;
4093 return poly;
4094 error:
4095 isl_qpolynomial_free(poly);
4096 return NULL;
4099 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *space,
4100 __isl_take isl_mat *div)
4102 isl_term *term;
4103 isl_size d;
4104 int n;
4106 d = isl_space_dim(space, isl_dim_all);
4107 if (d < 0 || !div)
4108 goto error;
4110 n = d + div->n_row;
4112 term = isl_calloc(space->ctx, struct isl_term,
4113 sizeof(struct isl_term) + (n - 1) * sizeof(int));
4114 if (!term)
4115 goto error;
4117 term->ref = 1;
4118 term->dim = space;
4119 term->div = div;
4120 isl_int_init(term->n);
4121 isl_int_init(term->d);
4123 return term;
4124 error:
4125 isl_space_free(space);
4126 isl_mat_free(div);
4127 return NULL;
4130 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
4132 if (!term)
4133 return NULL;
4135 term->ref++;
4136 return term;
4139 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
4141 int i;
4142 isl_term *dup;
4143 isl_size total;
4145 total = isl_term_dim(term, isl_dim_all);
4146 if (total < 0)
4147 return NULL;
4149 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
4150 if (!dup)
4151 return NULL;
4153 isl_int_set(dup->n, term->n);
4154 isl_int_set(dup->d, term->d);
4156 for (i = 0; i < total; ++i)
4157 dup->pow[i] = term->pow[i];
4159 return dup;
4162 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
4164 if (!term)
4165 return NULL;
4167 if (term->ref == 1)
4168 return term;
4169 term->ref--;
4170 return isl_term_dup(term);
4173 __isl_null isl_term *isl_term_free(__isl_take isl_term *term)
4175 if (!term)
4176 return NULL;
4178 if (--term->ref > 0)
4179 return NULL;
4181 isl_space_free(term->dim);
4182 isl_mat_free(term->div);
4183 isl_int_clear(term->n);
4184 isl_int_clear(term->d);
4185 free(term);
4187 return NULL;
4190 isl_size isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
4192 isl_size dim;
4194 if (!term)
4195 return isl_size_error;
4197 switch (type) {
4198 case isl_dim_param:
4199 case isl_dim_in:
4200 case isl_dim_out: return isl_space_dim(term->dim, type);
4201 case isl_dim_div: return term->div->n_row;
4202 case isl_dim_all: dim = isl_space_dim(term->dim, isl_dim_all);
4203 if (dim < 0)
4204 return isl_size_error;
4205 return dim + term->div->n_row;
4206 default: return isl_size_error;
4210 /* Return the space of "term".
4212 static __isl_keep isl_space *isl_term_peek_space(__isl_keep isl_term *term)
4214 return term ? term->dim : NULL;
4217 /* Return the offset of the first variable of type "type" within
4218 * the variables of "term".
4220 static isl_size isl_term_offset(__isl_keep isl_term *term,
4221 enum isl_dim_type type)
4223 isl_space *space;
4225 space = isl_term_peek_space(term);
4226 if (!space)
4227 return isl_size_error;
4229 switch (type) {
4230 case isl_dim_param:
4231 case isl_dim_set: return isl_space_offset(space, type);
4232 case isl_dim_div: return isl_space_dim(space, isl_dim_all);
4233 default:
4234 isl_die(isl_term_get_ctx(term), isl_error_invalid,
4235 "invalid dimension type", return isl_size_error);
4239 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
4241 return term ? term->dim->ctx : NULL;
4244 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
4246 if (!term)
4247 return;
4248 isl_int_set(*n, term->n);
4251 /* Return the coefficient of the term "term".
4253 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
4255 if (!term)
4256 return NULL;
4258 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
4259 term->n, term->d);
4262 #undef TYPE
4263 #define TYPE isl_term
4264 static
4265 #include "check_type_range_templ.c"
4267 isl_size isl_term_get_exp(__isl_keep isl_term *term,
4268 enum isl_dim_type type, unsigned pos)
4270 isl_size offset;
4272 if (isl_term_check_range(term, type, pos, 1) < 0)
4273 return isl_size_error;
4274 offset = isl_term_offset(term, type);
4275 if (offset < 0)
4276 return isl_size_error;
4278 return term->pow[offset + pos];
4281 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
4283 isl_local_space *ls;
4284 isl_aff *aff;
4286 if (isl_term_check_range(term, isl_dim_div, pos, 1) < 0)
4287 return NULL;
4289 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
4290 isl_mat_copy(term->div));
4291 aff = isl_aff_alloc(ls);
4292 if (!aff)
4293 return NULL;
4295 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
4297 aff = isl_aff_normalize(aff);
4299 return aff;
4302 __isl_give isl_term *isl_poly_foreach_term(__isl_keep isl_poly *poly,
4303 isl_stat (*fn)(__isl_take isl_term *term, void *user),
4304 __isl_take isl_term *term, void *user)
4306 int i;
4307 isl_bool is_zero, is_bad, is_cst;
4308 isl_poly_rec *rec;
4310 is_zero = isl_poly_is_zero(poly);
4311 if (is_zero < 0 || !term)
4312 goto error;
4314 if (is_zero)
4315 return term;
4317 is_cst = isl_poly_is_cst(poly);
4318 is_bad = isl_poly_is_nan(poly);
4319 if (is_bad >= 0 && !is_bad)
4320 is_bad = isl_poly_is_infty(poly);
4321 if (is_bad >= 0 && !is_bad)
4322 is_bad = isl_poly_is_neginfty(poly);
4323 if (is_cst < 0 || is_bad < 0)
4324 return isl_term_free(term);
4325 if (is_bad)
4326 isl_die(isl_term_get_ctx(term), isl_error_invalid,
4327 "cannot handle NaN/infty polynomial",
4328 return isl_term_free(term));
4330 if (is_cst) {
4331 isl_poly_cst *cst;
4332 cst = isl_poly_as_cst(poly);
4333 if (!cst)
4334 goto error;
4335 term = isl_term_cow(term);
4336 if (!term)
4337 goto error;
4338 isl_int_set(term->n, cst->n);
4339 isl_int_set(term->d, cst->d);
4340 if (fn(isl_term_copy(term), user) < 0)
4341 goto error;
4342 return term;
4345 rec = isl_poly_as_rec(poly);
4346 if (!rec)
4347 goto error;
4349 for (i = 0; i < rec->n; ++i) {
4350 term = isl_term_cow(term);
4351 if (!term)
4352 goto error;
4353 term->pow[poly->var] = i;
4354 term = isl_poly_foreach_term(rec->p[i], fn, term, user);
4355 if (!term)
4356 goto error;
4358 term = isl_term_cow(term);
4359 if (!term)
4360 return NULL;
4361 term->pow[poly->var] = 0;
4363 return term;
4364 error:
4365 isl_term_free(term);
4366 return NULL;
4369 isl_stat isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
4370 isl_stat (*fn)(__isl_take isl_term *term, void *user), void *user)
4372 isl_local *local;
4373 isl_term *term;
4375 if (!qp)
4376 return isl_stat_error;
4378 local = isl_qpolynomial_get_local(qp);
4379 term = isl_term_alloc(isl_space_copy(qp->dim), local);
4380 if (!term)
4381 return isl_stat_error;
4383 term = isl_poly_foreach_term(qp->poly, fn, term, user);
4385 isl_term_free(term);
4387 return term ? isl_stat_ok : isl_stat_error;
4390 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
4392 isl_poly *poly;
4393 isl_qpolynomial *qp;
4394 int i;
4395 isl_size n;
4397 n = isl_term_dim(term, isl_dim_all);
4398 if (n < 0)
4399 term = isl_term_free(term);
4400 if (!term)
4401 return NULL;
4403 poly = isl_poly_rat_cst(term->dim->ctx, term->n, term->d);
4404 for (i = 0; i < n; ++i) {
4405 if (!term->pow[i])
4406 continue;
4407 poly = isl_poly_mul(poly,
4408 isl_poly_var_pow(term->dim->ctx, i, term->pow[i]));
4411 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim),
4412 term->div->n_row, poly);
4413 if (!qp)
4414 goto error;
4415 isl_mat_free(qp->div);
4416 qp->div = isl_mat_copy(term->div);
4417 if (!qp->div)
4418 goto error;
4420 isl_term_free(term);
4421 return qp;
4422 error:
4423 isl_qpolynomial_free(qp);
4424 isl_term_free(term);
4425 return NULL;
4428 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
4429 __isl_take isl_space *space)
4431 int i;
4432 int extra;
4433 isl_size total, d_set, d_qp;
4435 if (!qp || !space)
4436 goto error;
4438 if (isl_space_is_equal(qp->dim, space)) {
4439 isl_space_free(space);
4440 return qp;
4443 qp = isl_qpolynomial_cow(qp);
4444 if (!qp)
4445 goto error;
4447 d_set = isl_space_dim(space, isl_dim_set);
4448 d_qp = isl_qpolynomial_domain_dim(qp, isl_dim_set);
4449 extra = d_set - d_qp;
4450 total = isl_space_dim(qp->dim, isl_dim_all);
4451 if (d_set < 0 || d_qp < 0 || total < 0)
4452 goto error;
4453 if (qp->div->n_row) {
4454 int *exp;
4456 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
4457 if (!exp)
4458 goto error;
4459 for (i = 0; i < qp->div->n_row; ++i)
4460 exp[i] = extra + i;
4461 qp->poly = expand(qp->poly, exp, total);
4462 free(exp);
4463 if (!qp->poly)
4464 goto error;
4466 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
4467 if (!qp->div)
4468 goto error;
4469 for (i = 0; i < qp->div->n_row; ++i)
4470 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
4472 isl_space_free(isl_qpolynomial_take_domain_space(qp));
4473 qp = isl_qpolynomial_restore_domain_space(qp, space);
4475 return qp;
4476 error:
4477 isl_space_free(space);
4478 isl_qpolynomial_free(qp);
4479 return NULL;
4482 /* For each parameter or variable that does not appear in qp,
4483 * first eliminate the variable from all constraints and then set it to zero.
4485 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
4486 __isl_keep isl_qpolynomial *qp)
4488 int *active = NULL;
4489 int i;
4490 isl_size d;
4491 isl_size nparam;
4492 isl_size nvar;
4494 d = isl_set_dim(set, isl_dim_all);
4495 if (d < 0 || !qp)
4496 goto error;
4498 active = isl_calloc_array(set->ctx, int, d);
4499 if (set_active(qp, active) < 0)
4500 goto error;
4502 for (i = 0; i < d; ++i)
4503 if (!active[i])
4504 break;
4506 if (i == d) {
4507 free(active);
4508 return set;
4511 nparam = isl_set_dim(set, isl_dim_param);
4512 nvar = isl_set_dim(set, isl_dim_set);
4513 if (nparam < 0 || nvar < 0)
4514 goto error;
4515 for (i = 0; i < nparam; ++i) {
4516 if (active[i])
4517 continue;
4518 set = isl_set_eliminate(set, isl_dim_param, i, 1);
4519 set = isl_set_fix_si(set, isl_dim_param, i, 0);
4521 for (i = 0; i < nvar; ++i) {
4522 if (active[nparam + i])
4523 continue;
4524 set = isl_set_eliminate(set, isl_dim_set, i, 1);
4525 set = isl_set_fix_si(set, isl_dim_set, i, 0);
4528 free(active);
4530 return set;
4531 error:
4532 free(active);
4533 isl_set_free(set);
4534 return NULL;
4537 struct isl_opt_data {
4538 isl_qpolynomial *qp;
4539 int first;
4540 isl_val *opt;
4541 int max;
4544 static isl_stat opt_fn(__isl_take isl_point *pnt, void *user)
4546 struct isl_opt_data *data = (struct isl_opt_data *)user;
4547 isl_val *val;
4549 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
4550 if (data->first) {
4551 data->first = 0;
4552 data->opt = val;
4553 } else if (data->max) {
4554 data->opt = isl_val_max(data->opt, val);
4555 } else {
4556 data->opt = isl_val_min(data->opt, val);
4559 return isl_stat_ok;
4562 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
4563 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
4565 struct isl_opt_data data = { NULL, 1, NULL, max };
4566 isl_bool is_cst;
4568 if (!set || !qp)
4569 goto error;
4571 is_cst = isl_poly_is_cst(qp->poly);
4572 if (is_cst < 0)
4573 goto error;
4574 if (is_cst) {
4575 isl_set_free(set);
4576 data.opt = isl_qpolynomial_get_constant_val(qp);
4577 isl_qpolynomial_free(qp);
4578 return data.opt;
4581 set = fix_inactive(set, qp);
4583 data.qp = qp;
4584 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
4585 goto error;
4587 if (data.first)
4588 data.opt = isl_val_zero(isl_set_get_ctx(set));
4590 isl_set_free(set);
4591 isl_qpolynomial_free(qp);
4592 return data.opt;
4593 error:
4594 isl_set_free(set);
4595 isl_qpolynomial_free(qp);
4596 isl_val_free(data.opt);
4597 return NULL;
4600 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
4601 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
4603 int i;
4604 int n_sub;
4605 isl_ctx *ctx;
4606 isl_space *space;
4607 isl_poly **subs;
4608 isl_mat *mat, *diag;
4610 qp = isl_qpolynomial_cow(qp);
4612 space = isl_qpolynomial_peek_domain_space(qp);
4613 if (isl_morph_check_applies(morph, space) < 0)
4614 goto error;
4616 ctx = isl_qpolynomial_get_ctx(qp);
4617 n_sub = morph->inv->n_row - 1;
4618 if (morph->inv->n_row != morph->inv->n_col)
4619 n_sub += qp->div->n_row;
4620 subs = isl_calloc_array(ctx, struct isl_poly *, n_sub);
4621 if (n_sub && !subs)
4622 goto error;
4624 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4625 subs[i] = isl_poly_from_affine(ctx, morph->inv->row[1 + i],
4626 morph->inv->row[0][0], morph->inv->n_col);
4627 if (morph->inv->n_row != morph->inv->n_col)
4628 for (i = 0; i < qp->div->n_row; ++i)
4629 subs[morph->inv->n_row - 1 + i] =
4630 isl_poly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4632 qp->poly = isl_poly_subs(qp->poly, 0, n_sub, subs);
4634 for (i = 0; i < n_sub; ++i)
4635 isl_poly_free(subs[i]);
4636 free(subs);
4638 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4639 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4640 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4641 mat = isl_mat_diagonal(mat, diag);
4642 qp->div = isl_mat_product(qp->div, mat);
4644 if (!qp->poly || !qp->div)
4645 goto error;
4647 isl_space_free(isl_qpolynomial_take_domain_space(qp));
4648 space = isl_space_copy(morph->ran->dim);
4649 qp = isl_qpolynomial_restore_domain_space(qp, space);
4651 isl_morph_free(morph);
4653 return qp;
4654 error:
4655 isl_qpolynomial_free(qp);
4656 isl_morph_free(morph);
4657 return NULL;
4660 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4661 __isl_take isl_union_pw_qpolynomial *upwqp1,
4662 __isl_take isl_union_pw_qpolynomial *upwqp2)
4664 return isl_union_pw_qpolynomial_match_bin_op(upwqp1, upwqp2,
4665 &isl_pw_qpolynomial_mul);
4668 /* Reorder the dimension of "qp" according to the given reordering.
4670 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4671 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4673 isl_space *space;
4674 isl_poly *poly;
4675 isl_local *local;
4677 if (!qp)
4678 goto error;
4680 r = isl_reordering_extend(r, qp->div->n_row);
4681 if (!r)
4682 goto error;
4684 local = isl_qpolynomial_take_local(qp);
4685 local = isl_local_reorder(local, isl_reordering_copy(r));
4686 qp = isl_qpolynomial_restore_local(qp, local);
4688 poly = isl_qpolynomial_take_poly(qp);
4689 poly = reorder(poly, r->pos);
4690 qp = isl_qpolynomial_restore_poly(qp, poly);
4692 space = isl_reordering_get_space(r);
4693 qp = isl_qpolynomial_reset_domain_space(qp, space);
4695 isl_reordering_free(r);
4696 return qp;
4697 error:
4698 isl_qpolynomial_free(qp);
4699 isl_reordering_free(r);
4700 return NULL;
4703 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4704 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4706 isl_space *domain_space;
4707 isl_bool equal_params;
4709 domain_space = isl_qpolynomial_peek_domain_space(qp);
4710 equal_params = isl_space_has_equal_params(domain_space, model);
4711 if (equal_params < 0)
4712 goto error;
4713 if (!equal_params) {
4714 isl_reordering *exp;
4716 exp = isl_parameter_alignment_reordering(domain_space, model);
4717 qp = isl_qpolynomial_realign_domain(qp, exp);
4720 isl_space_free(model);
4721 return qp;
4722 error:
4723 isl_space_free(model);
4724 isl_qpolynomial_free(qp);
4725 return NULL;
4728 struct isl_split_periods_data {
4729 int max_periods;
4730 isl_pw_qpolynomial *res;
4733 /* Create a slice where the integer division "div" has the fixed value "v".
4734 * In particular, if "div" refers to floor(f/m), then create a slice
4736 * m v <= f <= m v + (m - 1)
4738 * or
4740 * f - m v >= 0
4741 * -f + m v + (m - 1) >= 0
4743 static __isl_give isl_set *set_div_slice(__isl_take isl_space *space,
4744 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4746 isl_size total;
4747 isl_basic_set *bset = NULL;
4748 int k;
4750 total = isl_space_dim(space, isl_dim_all);
4751 if (total < 0 || !qp)
4752 goto error;
4754 bset = isl_basic_set_alloc_space(isl_space_copy(space), 0, 0, 2);
4756 k = isl_basic_set_alloc_inequality(bset);
4757 if (k < 0)
4758 goto error;
4759 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4760 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4762 k = isl_basic_set_alloc_inequality(bset);
4763 if (k < 0)
4764 goto error;
4765 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4766 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4767 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4768 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4770 isl_space_free(space);
4771 return isl_set_from_basic_set(bset);
4772 error:
4773 isl_basic_set_free(bset);
4774 isl_space_free(space);
4775 return NULL;
4778 static isl_stat split_periods(__isl_take isl_set *set,
4779 __isl_take isl_qpolynomial *qp, void *user);
4781 /* Create a slice of the domain "set" such that integer division "div"
4782 * has the fixed value "v" and add the results to data->res,
4783 * replacing the integer division by "v" in "qp".
4785 static isl_stat set_div(__isl_take isl_set *set,
4786 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4787 struct isl_split_periods_data *data)
4789 int i;
4790 isl_size div_pos;
4791 isl_set *slice;
4792 isl_poly *cst;
4794 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4795 set = isl_set_intersect(set, slice);
4797 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
4798 if (div_pos < 0)
4799 goto error;
4801 for (i = div + 1; i < qp->div->n_row; ++i) {
4802 if (isl_int_is_zero(qp->div->row[i][2 + div_pos + div]))
4803 continue;
4804 isl_int_addmul(qp->div->row[i][1],
4805 qp->div->row[i][2 + div_pos + div], v);
4806 isl_int_set_si(qp->div->row[i][2 + div_pos + div], 0);
4809 cst = isl_poly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4810 qp = substitute_div(qp, div, cst);
4812 return split_periods(set, qp, data);
4813 error:
4814 isl_set_free(set);
4815 isl_qpolynomial_free(qp);
4816 return isl_stat_error;
4819 /* Split the domain "set" such that integer division "div"
4820 * has a fixed value (ranging from "min" to "max") on each slice
4821 * and add the results to data->res.
4823 static isl_stat split_div(__isl_take isl_set *set,
4824 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4825 struct isl_split_periods_data *data)
4827 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4828 isl_set *set_i = isl_set_copy(set);
4829 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4831 if (set_div(set_i, qp_i, div, min, data) < 0)
4832 goto error;
4834 isl_set_free(set);
4835 isl_qpolynomial_free(qp);
4836 return isl_stat_ok;
4837 error:
4838 isl_set_free(set);
4839 isl_qpolynomial_free(qp);
4840 return isl_stat_error;
4843 /* If "qp" refers to any integer division
4844 * that can only attain "max_periods" distinct values on "set"
4845 * then split the domain along those distinct values.
4846 * Add the results (or the original if no splitting occurs)
4847 * to data->res.
4849 static isl_stat split_periods(__isl_take isl_set *set,
4850 __isl_take isl_qpolynomial *qp, void *user)
4852 int i;
4853 isl_pw_qpolynomial *pwqp;
4854 struct isl_split_periods_data *data;
4855 isl_int min, max;
4856 isl_size div_pos;
4857 isl_stat r = isl_stat_ok;
4859 data = (struct isl_split_periods_data *)user;
4861 if (!set || !qp)
4862 goto error;
4864 if (qp->div->n_row == 0) {
4865 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4866 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4867 return isl_stat_ok;
4870 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
4871 if (div_pos < 0)
4872 goto error;
4874 isl_int_init(min);
4875 isl_int_init(max);
4876 for (i = 0; i < qp->div->n_row; ++i) {
4877 enum isl_lp_result lp_res;
4879 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + div_pos,
4880 qp->div->n_row) != -1)
4881 continue;
4883 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4884 set->ctx->one, &min, NULL, NULL);
4885 if (lp_res == isl_lp_error)
4886 goto error2;
4887 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4888 continue;
4889 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4891 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4892 set->ctx->one, &max, NULL, NULL);
4893 if (lp_res == isl_lp_error)
4894 goto error2;
4895 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4896 continue;
4897 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4899 isl_int_sub(max, max, min);
4900 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4901 isl_int_add(max, max, min);
4902 break;
4906 if (i < qp->div->n_row) {
4907 r = split_div(set, qp, i, min, max, data);
4908 } else {
4909 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4910 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4913 isl_int_clear(max);
4914 isl_int_clear(min);
4916 return r;
4917 error2:
4918 isl_int_clear(max);
4919 isl_int_clear(min);
4920 error:
4921 isl_set_free(set);
4922 isl_qpolynomial_free(qp);
4923 return isl_stat_error;
4926 /* If any quasi-polynomial in pwqp refers to any integer division
4927 * that can only attain "max_periods" distinct values on its domain
4928 * then split the domain along those distinct values.
4930 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4931 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4933 struct isl_split_periods_data data;
4935 data.max_periods = max_periods;
4936 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4938 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4939 goto error;
4941 isl_pw_qpolynomial_free(pwqp);
4943 return data.res;
4944 error:
4945 isl_pw_qpolynomial_free(data.res);
4946 isl_pw_qpolynomial_free(pwqp);
4947 return NULL;
4950 /* Construct a piecewise quasipolynomial that is constant on the given
4951 * domain. In particular, it is
4952 * 0 if cst == 0
4953 * 1 if cst == 1
4954 * infinity if cst == -1
4956 * If cst == -1, then explicitly check whether the domain is empty and,
4957 * if so, return 0 instead.
4959 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4960 __isl_take isl_basic_set *bset, int cst)
4962 isl_space *space;
4963 isl_qpolynomial *qp;
4965 if (cst < 0 && isl_basic_set_is_empty(bset) == isl_bool_true)
4966 cst = 0;
4967 if (!bset)
4968 return NULL;
4970 bset = isl_basic_set_params(bset);
4971 space = isl_basic_set_get_space(bset);
4972 if (cst < 0)
4973 qp = isl_qpolynomial_infty_on_domain(space);
4974 else if (cst == 0)
4975 qp = isl_qpolynomial_zero_on_domain(space);
4976 else
4977 qp = isl_qpolynomial_one_on_domain(space);
4978 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4981 /* Internal data structure for multiplicative_call_factor_pw_qpolynomial.
4982 * "fn" is the function that is called on each factor.
4983 * "pwpq" collects the results.
4985 struct isl_multiplicative_call_data_pw_qpolynomial {
4986 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset);
4987 isl_pw_qpolynomial *pwqp;
4990 /* Call "fn" on "bset" and return the result,
4991 * but first check if "bset" has any redundant constraints or
4992 * implicit equality constraints.
4993 * If so, there may be further opportunities for detecting factors or
4994 * removing equality constraints, so recursively call
4995 * the top-level isl_basic_set_multiplicative_call.
4997 static __isl_give isl_pw_qpolynomial *multiplicative_call_base(
4998 __isl_take isl_basic_set *bset,
4999 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
5001 isl_size n1, n2, n_eq;
5003 n1 = isl_basic_set_n_constraint(bset);
5004 if (n1 < 0)
5005 bset = isl_basic_set_free(bset);
5006 bset = isl_basic_set_remove_redundancies(bset);
5007 bset = isl_basic_set_detect_equalities(bset);
5008 n2 = isl_basic_set_n_constraint(bset);
5009 n_eq = isl_basic_set_n_equality(bset);
5010 if (n2 < 0 || n_eq < 0)
5011 bset = isl_basic_set_free(bset);
5012 else if (n2 < n1 || n_eq > 0)
5013 return isl_basic_set_multiplicative_call(bset, fn);
5014 return fn(bset);
5017 /* isl_factorizer_every_factor_basic_set callback that applies
5018 * data->fn to the factor "bset" and multiplies in the result
5019 * in data->pwqp.
5021 static isl_bool multiplicative_call_factor_pw_qpolynomial(
5022 __isl_keep isl_basic_set *bset, void *user)
5024 struct isl_multiplicative_call_data_pw_qpolynomial *data = user;
5025 isl_pw_qpolynomial *res;
5027 bset = isl_basic_set_copy(bset);
5028 res = multiplicative_call_base(bset, data->fn);
5029 data->pwqp = isl_pw_qpolynomial_mul(data->pwqp, res);
5030 if (!data->pwqp)
5031 return isl_bool_error;
5033 return isl_bool_true;
5036 /* Factor bset, call fn on each of the factors and return the product.
5038 * If no factors can be found, simply call fn on the input.
5039 * Otherwise, construct the factors based on the factorizer,
5040 * call fn on each factor and compute the product.
5042 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
5043 __isl_take isl_basic_set *bset,
5044 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
5046 struct isl_multiplicative_call_data_pw_qpolynomial data = { fn };
5047 isl_space *space;
5048 isl_set *set;
5049 isl_factorizer *f;
5050 isl_qpolynomial *qp;
5051 isl_bool every;
5053 f = isl_basic_set_factorizer(bset);
5054 if (!f)
5055 goto error;
5056 if (f->n_group == 0) {
5057 isl_factorizer_free(f);
5058 return multiplicative_call_base(bset, fn);
5061 space = isl_basic_set_get_space(bset);
5062 space = isl_space_params(space);
5063 set = isl_set_universe(isl_space_copy(space));
5064 qp = isl_qpolynomial_one_on_domain(space);
5065 data.pwqp = isl_pw_qpolynomial_alloc(set, qp);
5067 every = isl_factorizer_every_factor_basic_set(f,
5068 &multiplicative_call_factor_pw_qpolynomial, &data);
5069 if (every < 0)
5070 data.pwqp = isl_pw_qpolynomial_free(data.pwqp);
5072 isl_basic_set_free(bset);
5073 isl_factorizer_free(f);
5075 return data.pwqp;
5076 error:
5077 isl_basic_set_free(bset);
5078 return NULL;
5081 /* Factor bset, call fn on each of the factors and return the product.
5082 * The function is assumed to evaluate to zero on empty domains,
5083 * to one on zero-dimensional domains and to infinity on unbounded domains
5084 * and will not be called explicitly on zero-dimensional or unbounded domains.
5086 * We first check for some special cases and remove all equalities.
5087 * Then we hand over control to compressed_multiplicative_call.
5089 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
5090 __isl_take isl_basic_set *bset,
5091 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
5093 isl_bool bounded;
5094 isl_size dim;
5095 isl_morph *morph;
5096 isl_pw_qpolynomial *pwqp;
5098 if (!bset)
5099 return NULL;
5101 if (isl_basic_set_plain_is_empty(bset))
5102 return constant_on_domain(bset, 0);
5104 dim = isl_basic_set_dim(bset, isl_dim_set);
5105 if (dim < 0)
5106 goto error;
5107 if (dim == 0)
5108 return constant_on_domain(bset, 1);
5110 bounded = isl_basic_set_is_bounded(bset);
5111 if (bounded < 0)
5112 goto error;
5113 if (!bounded)
5114 return constant_on_domain(bset, -1);
5116 if (bset->n_eq == 0)
5117 return compressed_multiplicative_call(bset, fn);
5119 morph = isl_basic_set_full_compression(bset);
5120 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
5122 pwqp = compressed_multiplicative_call(bset, fn);
5124 morph = isl_morph_dom_params(morph);
5125 morph = isl_morph_ran_params(morph);
5126 morph = isl_morph_inverse(morph);
5128 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
5130 return pwqp;
5131 error:
5132 isl_basic_set_free(bset);
5133 return NULL;
5136 /* Drop all floors in "qp", turning each integer division [a/m] into
5137 * a rational division a/m. If "down" is set, then the integer division
5138 * is replaced by (a-(m-1))/m instead.
5140 static __isl_give isl_qpolynomial *qp_drop_floors(
5141 __isl_take isl_qpolynomial *qp, int down)
5143 int i;
5144 isl_poly *s;
5146 if (!qp)
5147 return NULL;
5148 if (qp->div->n_row == 0)
5149 return qp;
5151 qp = isl_qpolynomial_cow(qp);
5152 if (!qp)
5153 return NULL;
5155 for (i = qp->div->n_row - 1; i >= 0; --i) {
5156 if (down) {
5157 isl_int_sub(qp->div->row[i][1],
5158 qp->div->row[i][1], qp->div->row[i][0]);
5159 isl_int_add_ui(qp->div->row[i][1],
5160 qp->div->row[i][1], 1);
5162 s = isl_poly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
5163 qp->div->row[i][0], qp->div->n_col - 1);
5164 qp = substitute_div(qp, i, s);
5165 if (!qp)
5166 return NULL;
5169 return qp;
5172 /* Drop all floors in "pwqp", turning each integer division [a/m] into
5173 * a rational division a/m.
5175 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
5176 __isl_take isl_pw_qpolynomial *pwqp)
5178 int i;
5180 if (!pwqp)
5181 return NULL;
5183 if (isl_pw_qpolynomial_is_zero(pwqp))
5184 return pwqp;
5186 pwqp = isl_pw_qpolynomial_cow(pwqp);
5187 if (!pwqp)
5188 return NULL;
5190 for (i = 0; i < pwqp->n; ++i) {
5191 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
5192 if (!pwqp->p[i].qp)
5193 goto error;
5196 return pwqp;
5197 error:
5198 isl_pw_qpolynomial_free(pwqp);
5199 return NULL;
5202 /* Adjust all the integer divisions in "qp" such that they are at least
5203 * one over the given orthant (identified by "signs"). This ensures
5204 * that they will still be non-negative even after subtracting (m-1)/m.
5206 * In particular, f is replaced by f' + v, changing f = [a/m]
5207 * to f' = [(a - m v)/m].
5208 * If the constant term k in a is smaller than m,
5209 * the constant term of v is set to floor(k/m) - 1.
5210 * For any other term, if the coefficient c and the variable x have
5211 * the same sign, then no changes are needed.
5212 * Otherwise, if the variable is positive (and c is negative),
5213 * then the coefficient of x in v is set to floor(c/m).
5214 * If the variable is negative (and c is positive),
5215 * then the coefficient of x in v is set to ceil(c/m).
5217 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
5218 int *signs)
5220 int i, j;
5221 isl_size div_pos;
5222 isl_vec *v = NULL;
5223 isl_poly *s;
5225 qp = isl_qpolynomial_cow(qp);
5226 div_pos = isl_qpolynomial_domain_var_offset(qp, isl_dim_div);
5227 if (div_pos < 0)
5228 return isl_qpolynomial_free(qp);
5229 qp->div = isl_mat_cow(qp->div);
5230 if (!qp->div)
5231 goto error;
5233 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
5235 for (i = 0; i < qp->div->n_row; ++i) {
5236 isl_int *row = qp->div->row[i];
5237 v = isl_vec_clr(v);
5238 if (!v)
5239 goto error;
5240 if (isl_int_lt(row[1], row[0])) {
5241 isl_int_fdiv_q(v->el[0], row[1], row[0]);
5242 isl_int_sub_ui(v->el[0], v->el[0], 1);
5243 isl_int_submul(row[1], row[0], v->el[0]);
5245 for (j = 0; j < div_pos; ++j) {
5246 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
5247 continue;
5248 if (signs[j] < 0)
5249 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
5250 else
5251 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
5252 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
5254 for (j = 0; j < i; ++j) {
5255 if (isl_int_sgn(row[2 + div_pos + j]) >= 0)
5256 continue;
5257 isl_int_fdiv_q(v->el[1 + div_pos + j],
5258 row[2 + div_pos + j], row[0]);
5259 isl_int_submul(row[2 + div_pos + j],
5260 row[0], v->el[1 + div_pos + j]);
5262 for (j = i + 1; j < qp->div->n_row; ++j) {
5263 if (isl_int_is_zero(qp->div->row[j][2 + div_pos + i]))
5264 continue;
5265 isl_seq_combine(qp->div->row[j] + 1,
5266 qp->div->ctx->one, qp->div->row[j] + 1,
5267 qp->div->row[j][2 + div_pos + i], v->el,
5268 v->size);
5270 isl_int_set_si(v->el[1 + div_pos + i], 1);
5271 s = isl_poly_from_affine(qp->dim->ctx, v->el,
5272 qp->div->ctx->one, v->size);
5273 qp->poly = isl_poly_subs(qp->poly, div_pos + i, 1, &s);
5274 isl_poly_free(s);
5275 if (!qp->poly)
5276 goto error;
5279 isl_vec_free(v);
5280 return qp;
5281 error:
5282 isl_vec_free(v);
5283 isl_qpolynomial_free(qp);
5284 return NULL;
5287 struct isl_to_poly_data {
5288 int sign;
5289 isl_pw_qpolynomial *res;
5290 isl_qpolynomial *qp;
5293 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
5294 * We first make all integer divisions positive and then split the
5295 * quasipolynomials into terms with sign data->sign (the direction
5296 * of the requested approximation) and terms with the opposite sign.
5297 * In the first set of terms, each integer division [a/m] is
5298 * overapproximated by a/m, while in the second it is underapproximated
5299 * by (a-(m-1))/m.
5301 static isl_stat to_polynomial_on_orthant(__isl_take isl_set *orthant,
5302 int *signs, void *user)
5304 struct isl_to_poly_data *data = user;
5305 isl_pw_qpolynomial *t;
5306 isl_qpolynomial *qp, *up, *down;
5308 qp = isl_qpolynomial_copy(data->qp);
5309 qp = make_divs_pos(qp, signs);
5311 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
5312 up = qp_drop_floors(up, 0);
5313 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
5314 down = qp_drop_floors(down, 1);
5316 isl_qpolynomial_free(qp);
5317 qp = isl_qpolynomial_add(up, down);
5319 t = isl_pw_qpolynomial_alloc(orthant, qp);
5320 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
5322 return isl_stat_ok;
5325 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
5326 * the polynomial will be an overapproximation. If "sign" is negative,
5327 * it will be an underapproximation. If "sign" is zero, the approximation
5328 * will lie somewhere in between.
5330 * In particular, is sign == 0, we simply drop the floors, turning
5331 * the integer divisions into rational divisions.
5332 * Otherwise, we split the domains into orthants, make all integer divisions
5333 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
5334 * depending on the requested sign and the sign of the term in which
5335 * the integer division appears.
5337 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
5338 __isl_take isl_pw_qpolynomial *pwqp, int sign)
5340 int i;
5341 struct isl_to_poly_data data;
5343 if (sign == 0)
5344 return pwqp_drop_floors(pwqp);
5346 if (!pwqp)
5347 return NULL;
5349 data.sign = sign;
5350 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
5352 for (i = 0; i < pwqp->n; ++i) {
5353 if (pwqp->p[i].qp->div->n_row == 0) {
5354 isl_pw_qpolynomial *t;
5355 t = isl_pw_qpolynomial_alloc(
5356 isl_set_copy(pwqp->p[i].set),
5357 isl_qpolynomial_copy(pwqp->p[i].qp));
5358 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
5359 continue;
5361 data.qp = pwqp->p[i].qp;
5362 if (isl_set_foreach_orthant(pwqp->p[i].set,
5363 &to_polynomial_on_orthant, &data) < 0)
5364 goto error;
5367 isl_pw_qpolynomial_free(pwqp);
5369 return data.res;
5370 error:
5371 isl_pw_qpolynomial_free(pwqp);
5372 isl_pw_qpolynomial_free(data.res);
5373 return NULL;
5376 static __isl_give isl_pw_qpolynomial *poly_entry(
5377 __isl_take isl_pw_qpolynomial *pwqp, void *user)
5379 int *sign = user;
5381 return isl_pw_qpolynomial_to_polynomial(pwqp, *sign);
5384 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
5385 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
5387 return isl_union_pw_qpolynomial_transform_inplace(upwqp,
5388 &poly_entry, &sign);
5391 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
5392 __isl_take isl_qpolynomial *qp)
5394 isl_local_space *ls;
5395 isl_vec *vec;
5396 isl_aff *aff;
5397 isl_basic_map *bmap;
5398 isl_bool is_affine;
5400 if (!qp)
5401 return NULL;
5402 is_affine = isl_poly_is_affine(qp->poly);
5403 if (is_affine < 0)
5404 goto error;
5405 if (!is_affine)
5406 isl_die(qp->dim->ctx, isl_error_invalid,
5407 "input quasi-polynomial not affine", goto error);
5408 ls = isl_qpolynomial_get_domain_local_space(qp);
5409 vec = isl_qpolynomial_extract_affine(qp);
5410 aff = isl_aff_alloc_vec(ls, vec);
5411 bmap = isl_basic_map_from_aff(aff);
5412 isl_qpolynomial_free(qp);
5413 return bmap;
5414 error:
5415 isl_qpolynomial_free(qp);
5416 return NULL;