isl_pw_templ.c: extract out isl_pw_eval.c
[isl.git] / isl_polynomial.c
blob4813a489f04a8072e2522216bc8f120fe89bbfe4
1 /*
2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
11 #include <stdlib.h>
12 #define ISL_DIM_H
13 #include <isl_ctx_private.h>
14 #include <isl_map_private.h>
15 #include <isl_factorization.h>
16 #include <isl_lp_private.h>
17 #include <isl_seq.h>
18 #include <isl_union_map_private.h>
19 #include <isl_constraint_private.h>
20 #include <isl_polynomial_private.h>
21 #include <isl_point_private.h>
22 #include <isl_space_private.h>
23 #include <isl_mat_private.h>
24 #include <isl_vec_private.h>
25 #include <isl_range.h>
26 #include <isl_local.h>
27 #include <isl_local_space_private.h>
28 #include <isl_aff_private.h>
29 #include <isl_val_private.h>
30 #include <isl_config.h>
32 static unsigned pos(__isl_keep isl_space *dim, enum isl_dim_type type)
34 switch (type) {
35 case isl_dim_param: return 0;
36 case isl_dim_in: return dim->nparam;
37 case isl_dim_out: return dim->nparam + dim->n_in;
38 default: return 0;
42 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
44 if (!up)
45 return -1;
47 return up->var < 0;
50 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
52 if (!up)
53 return NULL;
55 isl_assert(up->ctx, up->var < 0, return NULL);
57 return (struct isl_upoly_cst *)up;
60 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
62 if (!up)
63 return NULL;
65 isl_assert(up->ctx, up->var >= 0, return NULL);
67 return (struct isl_upoly_rec *)up;
70 /* Compare two polynomials.
72 * Return -1 if "up1" is "smaller" than "up2", 1 if "up1" is "greater"
73 * than "up2" and 0 if they are equal.
75 static int isl_upoly_plain_cmp(__isl_keep struct isl_upoly *up1,
76 __isl_keep struct isl_upoly *up2)
78 int i;
79 struct isl_upoly_rec *rec1, *rec2;
81 if (up1 == up2)
82 return 0;
83 if (!up1)
84 return -1;
85 if (!up2)
86 return 1;
87 if (up1->var != up2->var)
88 return up1->var - up2->var;
90 if (isl_upoly_is_cst(up1)) {
91 struct isl_upoly_cst *cst1, *cst2;
92 int cmp;
94 cst1 = isl_upoly_as_cst(up1);
95 cst2 = isl_upoly_as_cst(up2);
96 if (!cst1 || !cst2)
97 return 0;
98 cmp = isl_int_cmp(cst1->n, cst2->n);
99 if (cmp != 0)
100 return cmp;
101 return isl_int_cmp(cst1->d, cst2->d);
104 rec1 = isl_upoly_as_rec(up1);
105 rec2 = isl_upoly_as_rec(up2);
106 if (!rec1 || !rec2)
107 return 0;
109 if (rec1->n != rec2->n)
110 return rec1->n - rec2->n;
112 for (i = 0; i < rec1->n; ++i) {
113 int cmp = isl_upoly_plain_cmp(rec1->p[i], rec2->p[i]);
114 if (cmp != 0)
115 return cmp;
118 return 0;
121 isl_bool isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
122 __isl_keep struct isl_upoly *up2)
124 int i;
125 struct isl_upoly_rec *rec1, *rec2;
127 if (!up1 || !up2)
128 return isl_bool_error;
129 if (up1 == up2)
130 return isl_bool_true;
131 if (up1->var != up2->var)
132 return isl_bool_false;
133 if (isl_upoly_is_cst(up1)) {
134 struct isl_upoly_cst *cst1, *cst2;
135 cst1 = isl_upoly_as_cst(up1);
136 cst2 = isl_upoly_as_cst(up2);
137 if (!cst1 || !cst2)
138 return isl_bool_error;
139 return isl_int_eq(cst1->n, cst2->n) &&
140 isl_int_eq(cst1->d, cst2->d);
143 rec1 = isl_upoly_as_rec(up1);
144 rec2 = isl_upoly_as_rec(up2);
145 if (!rec1 || !rec2)
146 return isl_bool_error;
148 if (rec1->n != rec2->n)
149 return isl_bool_false;
151 for (i = 0; i < rec1->n; ++i) {
152 isl_bool eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
153 if (eq < 0 || !eq)
154 return eq;
157 return isl_bool_true;
160 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
162 struct isl_upoly_cst *cst;
164 if (!up)
165 return -1;
166 if (!isl_upoly_is_cst(up))
167 return 0;
169 cst = isl_upoly_as_cst(up);
170 if (!cst)
171 return -1;
173 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
176 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
178 struct isl_upoly_cst *cst;
180 if (!up)
181 return 0;
182 if (!isl_upoly_is_cst(up))
183 return 0;
185 cst = isl_upoly_as_cst(up);
186 if (!cst)
187 return 0;
189 return isl_int_sgn(cst->n);
192 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
194 struct isl_upoly_cst *cst;
196 if (!up)
197 return -1;
198 if (!isl_upoly_is_cst(up))
199 return 0;
201 cst = isl_upoly_as_cst(up);
202 if (!cst)
203 return -1;
205 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
208 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
210 struct isl_upoly_cst *cst;
212 if (!up)
213 return -1;
214 if (!isl_upoly_is_cst(up))
215 return 0;
217 cst = isl_upoly_as_cst(up);
218 if (!cst)
219 return -1;
221 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
224 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
226 struct isl_upoly_cst *cst;
228 if (!up)
229 return -1;
230 if (!isl_upoly_is_cst(up))
231 return 0;
233 cst = isl_upoly_as_cst(up);
234 if (!cst)
235 return -1;
237 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
240 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
242 struct isl_upoly_cst *cst;
244 if (!up)
245 return -1;
246 if (!isl_upoly_is_cst(up))
247 return 0;
249 cst = isl_upoly_as_cst(up);
250 if (!cst)
251 return -1;
253 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
256 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
258 struct isl_upoly_cst *cst;
260 if (!up)
261 return -1;
262 if (!isl_upoly_is_cst(up))
263 return 0;
265 cst = isl_upoly_as_cst(up);
266 if (!cst)
267 return -1;
269 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
272 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
274 struct isl_upoly_cst *cst;
276 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
277 if (!cst)
278 return NULL;
280 cst->up.ref = 1;
281 cst->up.ctx = ctx;
282 isl_ctx_ref(ctx);
283 cst->up.var = -1;
285 isl_int_init(cst->n);
286 isl_int_init(cst->d);
288 return cst;
291 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
293 struct isl_upoly_cst *cst;
295 cst = isl_upoly_cst_alloc(ctx);
296 if (!cst)
297 return NULL;
299 isl_int_set_si(cst->n, 0);
300 isl_int_set_si(cst->d, 1);
302 return &cst->up;
305 __isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
307 struct isl_upoly_cst *cst;
309 cst = isl_upoly_cst_alloc(ctx);
310 if (!cst)
311 return NULL;
313 isl_int_set_si(cst->n, 1);
314 isl_int_set_si(cst->d, 1);
316 return &cst->up;
319 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
321 struct isl_upoly_cst *cst;
323 cst = isl_upoly_cst_alloc(ctx);
324 if (!cst)
325 return NULL;
327 isl_int_set_si(cst->n, 1);
328 isl_int_set_si(cst->d, 0);
330 return &cst->up;
333 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
335 struct isl_upoly_cst *cst;
337 cst = isl_upoly_cst_alloc(ctx);
338 if (!cst)
339 return NULL;
341 isl_int_set_si(cst->n, -1);
342 isl_int_set_si(cst->d, 0);
344 return &cst->up;
347 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
349 struct isl_upoly_cst *cst;
351 cst = isl_upoly_cst_alloc(ctx);
352 if (!cst)
353 return NULL;
355 isl_int_set_si(cst->n, 0);
356 isl_int_set_si(cst->d, 0);
358 return &cst->up;
361 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
362 isl_int n, isl_int d)
364 struct isl_upoly_cst *cst;
366 cst = isl_upoly_cst_alloc(ctx);
367 if (!cst)
368 return NULL;
370 isl_int_set(cst->n, n);
371 isl_int_set(cst->d, d);
373 return &cst->up;
376 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
377 int var, int size)
379 struct isl_upoly_rec *rec;
381 isl_assert(ctx, var >= 0, return NULL);
382 isl_assert(ctx, size >= 0, return NULL);
383 rec = isl_calloc(ctx, struct isl_upoly_rec,
384 sizeof(struct isl_upoly_rec) +
385 size * sizeof(struct isl_upoly *));
386 if (!rec)
387 return NULL;
389 rec->up.ref = 1;
390 rec->up.ctx = ctx;
391 isl_ctx_ref(ctx);
392 rec->up.var = var;
394 rec->n = 0;
395 rec->size = size;
397 return rec;
400 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
401 __isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
403 qp = isl_qpolynomial_cow(qp);
404 if (!qp || !dim)
405 goto error;
407 isl_space_free(qp->dim);
408 qp->dim = dim;
410 return qp;
411 error:
412 isl_qpolynomial_free(qp);
413 isl_space_free(dim);
414 return NULL;
417 /* Reset the space of "qp". This function is called from isl_pw_templ.c
418 * and doesn't know if the space of an element object is represented
419 * directly or through its domain. It therefore passes along both.
421 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
422 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
423 __isl_take isl_space *domain)
425 isl_space_free(space);
426 return isl_qpolynomial_reset_domain_space(qp, domain);
429 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
431 return qp ? qp->dim->ctx : NULL;
434 __isl_give isl_space *isl_qpolynomial_get_domain_space(
435 __isl_keep isl_qpolynomial *qp)
437 return qp ? isl_space_copy(qp->dim) : NULL;
440 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
442 isl_space *space;
443 if (!qp)
444 return NULL;
445 space = isl_space_copy(qp->dim);
446 space = isl_space_from_domain(space);
447 space = isl_space_add_dims(space, isl_dim_out, 1);
448 return space;
451 /* Return the number of variables of the given type in the domain of "qp".
453 unsigned isl_qpolynomial_domain_dim(__isl_keep isl_qpolynomial *qp,
454 enum isl_dim_type type)
456 if (!qp)
457 return 0;
458 if (type == isl_dim_div)
459 return qp->div->n_row;
460 if (type == isl_dim_all)
461 return isl_space_dim(qp->dim, isl_dim_all) +
462 isl_qpolynomial_domain_dim(qp, isl_dim_div);
463 return isl_space_dim(qp->dim, type);
466 /* Externally, an isl_qpolynomial has a map space, but internally, the
467 * ls field corresponds to the domain of that space.
469 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
470 enum isl_dim_type type)
472 if (!qp)
473 return 0;
474 if (type == isl_dim_out)
475 return 1;
476 if (type == isl_dim_in)
477 type = isl_dim_set;
478 return isl_qpolynomial_domain_dim(qp, type);
481 /* Return the offset of the first coefficient of type "type" in
482 * the domain of "qp".
484 unsigned isl_qpolynomial_domain_offset(__isl_keep isl_qpolynomial *qp,
485 enum isl_dim_type type)
487 if (!qp)
488 return 0;
489 switch (type) {
490 case isl_dim_cst:
491 return 0;
492 case isl_dim_param:
493 case isl_dim_set:
494 return 1 + isl_space_offset(qp->dim, type);
495 case isl_dim_div:
496 return 1 + isl_space_dim(qp->dim, isl_dim_all);
497 default:
498 return 0;
502 isl_bool isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
504 return qp ? isl_upoly_is_zero(qp->upoly) : isl_bool_error;
507 isl_bool isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
509 return qp ? isl_upoly_is_one(qp->upoly) : isl_bool_error;
512 isl_bool isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
514 return qp ? isl_upoly_is_nan(qp->upoly) : isl_bool_error;
517 isl_bool isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
519 return qp ? isl_upoly_is_infty(qp->upoly) : isl_bool_error;
522 isl_bool isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
524 return qp ? isl_upoly_is_neginfty(qp->upoly) : isl_bool_error;
527 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
529 return qp ? isl_upoly_sgn(qp->upoly) : 0;
532 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
534 isl_int_clear(cst->n);
535 isl_int_clear(cst->d);
538 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
540 int i;
542 for (i = 0; i < rec->n; ++i)
543 isl_upoly_free(rec->p[i]);
546 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
548 if (!up)
549 return NULL;
551 up->ref++;
552 return up;
555 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
557 struct isl_upoly_cst *cst;
558 struct isl_upoly_cst *dup;
560 cst = isl_upoly_as_cst(up);
561 if (!cst)
562 return NULL;
564 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
565 if (!dup)
566 return NULL;
567 isl_int_set(dup->n, cst->n);
568 isl_int_set(dup->d, cst->d);
570 return &dup->up;
573 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
575 int i;
576 struct isl_upoly_rec *rec;
577 struct isl_upoly_rec *dup;
579 rec = isl_upoly_as_rec(up);
580 if (!rec)
581 return NULL;
583 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
584 if (!dup)
585 return NULL;
587 for (i = 0; i < rec->n; ++i) {
588 dup->p[i] = isl_upoly_copy(rec->p[i]);
589 if (!dup->p[i])
590 goto error;
591 dup->n++;
594 return &dup->up;
595 error:
596 isl_upoly_free(&dup->up);
597 return NULL;
600 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
602 if (!up)
603 return NULL;
605 if (isl_upoly_is_cst(up))
606 return isl_upoly_dup_cst(up);
607 else
608 return isl_upoly_dup_rec(up);
611 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
613 if (!up)
614 return NULL;
616 if (up->ref == 1)
617 return up;
618 up->ref--;
619 return isl_upoly_dup(up);
622 __isl_null struct isl_upoly *isl_upoly_free(__isl_take struct isl_upoly *up)
624 if (!up)
625 return NULL;
627 if (--up->ref > 0)
628 return NULL;
630 if (up->var < 0)
631 upoly_free_cst((struct isl_upoly_cst *)up);
632 else
633 upoly_free_rec((struct isl_upoly_rec *)up);
635 isl_ctx_deref(up->ctx);
636 free(up);
637 return NULL;
640 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
642 isl_int gcd;
644 isl_int_init(gcd);
645 isl_int_gcd(gcd, cst->n, cst->d);
646 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
647 isl_int_divexact(cst->n, cst->n, gcd);
648 isl_int_divexact(cst->d, cst->d, gcd);
650 isl_int_clear(gcd);
653 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
654 __isl_take struct isl_upoly *up2)
656 struct isl_upoly_cst *cst1;
657 struct isl_upoly_cst *cst2;
659 up1 = isl_upoly_cow(up1);
660 if (!up1 || !up2)
661 goto error;
663 cst1 = isl_upoly_as_cst(up1);
664 cst2 = isl_upoly_as_cst(up2);
666 if (isl_int_eq(cst1->d, cst2->d))
667 isl_int_add(cst1->n, cst1->n, cst2->n);
668 else {
669 isl_int_mul(cst1->n, cst1->n, cst2->d);
670 isl_int_addmul(cst1->n, cst2->n, cst1->d);
671 isl_int_mul(cst1->d, cst1->d, cst2->d);
674 isl_upoly_cst_reduce(cst1);
676 isl_upoly_free(up2);
677 return up1;
678 error:
679 isl_upoly_free(up1);
680 isl_upoly_free(up2);
681 return NULL;
684 static __isl_give struct isl_upoly *replace_by_zero(
685 __isl_take struct isl_upoly *up)
687 struct isl_ctx *ctx;
689 if (!up)
690 return NULL;
691 ctx = up->ctx;
692 isl_upoly_free(up);
693 return isl_upoly_zero(ctx);
696 static __isl_give struct isl_upoly *replace_by_constant_term(
697 __isl_take struct isl_upoly *up)
699 struct isl_upoly_rec *rec;
700 struct isl_upoly *cst;
702 if (!up)
703 return NULL;
705 rec = isl_upoly_as_rec(up);
706 if (!rec)
707 goto error;
708 cst = isl_upoly_copy(rec->p[0]);
709 isl_upoly_free(up);
710 return cst;
711 error:
712 isl_upoly_free(up);
713 return NULL;
716 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
717 __isl_take struct isl_upoly *up2)
719 int i;
720 struct isl_upoly_rec *rec1, *rec2;
722 if (!up1 || !up2)
723 goto error;
725 if (isl_upoly_is_nan(up1)) {
726 isl_upoly_free(up2);
727 return up1;
730 if (isl_upoly_is_nan(up2)) {
731 isl_upoly_free(up1);
732 return up2;
735 if (isl_upoly_is_zero(up1)) {
736 isl_upoly_free(up1);
737 return up2;
740 if (isl_upoly_is_zero(up2)) {
741 isl_upoly_free(up2);
742 return up1;
745 if (up1->var < up2->var)
746 return isl_upoly_sum(up2, up1);
748 if (up2->var < up1->var) {
749 struct isl_upoly_rec *rec;
750 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
751 isl_upoly_free(up1);
752 return up2;
754 up1 = isl_upoly_cow(up1);
755 rec = isl_upoly_as_rec(up1);
756 if (!rec)
757 goto error;
758 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
759 if (rec->n == 1)
760 up1 = replace_by_constant_term(up1);
761 return up1;
764 if (isl_upoly_is_cst(up1))
765 return isl_upoly_sum_cst(up1, up2);
767 rec1 = isl_upoly_as_rec(up1);
768 rec2 = isl_upoly_as_rec(up2);
769 if (!rec1 || !rec2)
770 goto error;
772 if (rec1->n < rec2->n)
773 return isl_upoly_sum(up2, up1);
775 up1 = isl_upoly_cow(up1);
776 rec1 = isl_upoly_as_rec(up1);
777 if (!rec1)
778 goto error;
780 for (i = rec2->n - 1; i >= 0; --i) {
781 rec1->p[i] = isl_upoly_sum(rec1->p[i],
782 isl_upoly_copy(rec2->p[i]));
783 if (!rec1->p[i])
784 goto error;
785 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
786 isl_upoly_free(rec1->p[i]);
787 rec1->n--;
791 if (rec1->n == 0)
792 up1 = replace_by_zero(up1);
793 else if (rec1->n == 1)
794 up1 = replace_by_constant_term(up1);
796 isl_upoly_free(up2);
798 return up1;
799 error:
800 isl_upoly_free(up1);
801 isl_upoly_free(up2);
802 return NULL;
805 __isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
806 __isl_take struct isl_upoly *up, isl_int v)
808 struct isl_upoly_cst *cst;
810 up = isl_upoly_cow(up);
811 if (!up)
812 return NULL;
814 cst = isl_upoly_as_cst(up);
816 isl_int_addmul(cst->n, cst->d, v);
818 return up;
821 __isl_give struct isl_upoly *isl_upoly_add_isl_int(
822 __isl_take struct isl_upoly *up, isl_int v)
824 struct isl_upoly_rec *rec;
826 if (!up)
827 return NULL;
829 if (isl_upoly_is_cst(up))
830 return isl_upoly_cst_add_isl_int(up, v);
832 up = isl_upoly_cow(up);
833 rec = isl_upoly_as_rec(up);
834 if (!rec)
835 goto error;
837 rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
838 if (!rec->p[0])
839 goto error;
841 return up;
842 error:
843 isl_upoly_free(up);
844 return NULL;
847 __isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
848 __isl_take struct isl_upoly *up, isl_int v)
850 struct isl_upoly_cst *cst;
852 if (isl_upoly_is_zero(up))
853 return up;
855 up = isl_upoly_cow(up);
856 if (!up)
857 return NULL;
859 cst = isl_upoly_as_cst(up);
861 isl_int_mul(cst->n, cst->n, v);
863 return up;
866 __isl_give struct isl_upoly *isl_upoly_mul_isl_int(
867 __isl_take struct isl_upoly *up, isl_int v)
869 int i;
870 struct isl_upoly_rec *rec;
872 if (!up)
873 return NULL;
875 if (isl_upoly_is_cst(up))
876 return isl_upoly_cst_mul_isl_int(up, v);
878 up = isl_upoly_cow(up);
879 rec = isl_upoly_as_rec(up);
880 if (!rec)
881 goto error;
883 for (i = 0; i < rec->n; ++i) {
884 rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
885 if (!rec->p[i])
886 goto error;
889 return up;
890 error:
891 isl_upoly_free(up);
892 return NULL;
895 /* Multiply the constant polynomial "up" by "v".
897 static __isl_give struct isl_upoly *isl_upoly_cst_scale_val(
898 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
900 struct isl_upoly_cst *cst;
902 if (isl_upoly_is_zero(up))
903 return up;
905 up = isl_upoly_cow(up);
906 if (!up)
907 return NULL;
909 cst = isl_upoly_as_cst(up);
911 isl_int_mul(cst->n, cst->n, v->n);
912 isl_int_mul(cst->d, cst->d, v->d);
913 isl_upoly_cst_reduce(cst);
915 return up;
918 /* Multiply the polynomial "up" by "v".
920 static __isl_give struct isl_upoly *isl_upoly_scale_val(
921 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
923 int i;
924 struct isl_upoly_rec *rec;
926 if (!up)
927 return NULL;
929 if (isl_upoly_is_cst(up))
930 return isl_upoly_cst_scale_val(up, v);
932 up = isl_upoly_cow(up);
933 rec = isl_upoly_as_rec(up);
934 if (!rec)
935 goto error;
937 for (i = 0; i < rec->n; ++i) {
938 rec->p[i] = isl_upoly_scale_val(rec->p[i], v);
939 if (!rec->p[i])
940 goto error;
943 return up;
944 error:
945 isl_upoly_free(up);
946 return NULL;
949 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
950 __isl_take struct isl_upoly *up2)
952 struct isl_upoly_cst *cst1;
953 struct isl_upoly_cst *cst2;
955 up1 = isl_upoly_cow(up1);
956 if (!up1 || !up2)
957 goto error;
959 cst1 = isl_upoly_as_cst(up1);
960 cst2 = isl_upoly_as_cst(up2);
962 isl_int_mul(cst1->n, cst1->n, cst2->n);
963 isl_int_mul(cst1->d, cst1->d, cst2->d);
965 isl_upoly_cst_reduce(cst1);
967 isl_upoly_free(up2);
968 return up1;
969 error:
970 isl_upoly_free(up1);
971 isl_upoly_free(up2);
972 return NULL;
975 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
976 __isl_take struct isl_upoly *up2)
978 struct isl_upoly_rec *rec1;
979 struct isl_upoly_rec *rec2;
980 struct isl_upoly_rec *res = NULL;
981 int i, j;
982 int size;
984 rec1 = isl_upoly_as_rec(up1);
985 rec2 = isl_upoly_as_rec(up2);
986 if (!rec1 || !rec2)
987 goto error;
988 size = rec1->n + rec2->n - 1;
989 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
990 if (!res)
991 goto error;
993 for (i = 0; i < rec1->n; ++i) {
994 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
995 isl_upoly_copy(rec1->p[i]));
996 if (!res->p[i])
997 goto error;
998 res->n++;
1000 for (; i < size; ++i) {
1001 res->p[i] = isl_upoly_zero(up1->ctx);
1002 if (!res->p[i])
1003 goto error;
1004 res->n++;
1006 for (i = 0; i < rec1->n; ++i) {
1007 for (j = 1; j < rec2->n; ++j) {
1008 struct isl_upoly *up;
1009 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
1010 isl_upoly_copy(rec1->p[i]));
1011 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
1012 if (!res->p[i + j])
1013 goto error;
1017 isl_upoly_free(up1);
1018 isl_upoly_free(up2);
1020 return &res->up;
1021 error:
1022 isl_upoly_free(up1);
1023 isl_upoly_free(up2);
1024 isl_upoly_free(&res->up);
1025 return NULL;
1028 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
1029 __isl_take struct isl_upoly *up2)
1031 if (!up1 || !up2)
1032 goto error;
1034 if (isl_upoly_is_nan(up1)) {
1035 isl_upoly_free(up2);
1036 return up1;
1039 if (isl_upoly_is_nan(up2)) {
1040 isl_upoly_free(up1);
1041 return up2;
1044 if (isl_upoly_is_zero(up1)) {
1045 isl_upoly_free(up2);
1046 return up1;
1049 if (isl_upoly_is_zero(up2)) {
1050 isl_upoly_free(up1);
1051 return up2;
1054 if (isl_upoly_is_one(up1)) {
1055 isl_upoly_free(up1);
1056 return up2;
1059 if (isl_upoly_is_one(up2)) {
1060 isl_upoly_free(up2);
1061 return up1;
1064 if (up1->var < up2->var)
1065 return isl_upoly_mul(up2, up1);
1067 if (up2->var < up1->var) {
1068 int i;
1069 struct isl_upoly_rec *rec;
1070 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
1071 isl_ctx *ctx = up1->ctx;
1072 isl_upoly_free(up1);
1073 isl_upoly_free(up2);
1074 return isl_upoly_nan(ctx);
1076 up1 = isl_upoly_cow(up1);
1077 rec = isl_upoly_as_rec(up1);
1078 if (!rec)
1079 goto error;
1081 for (i = 0; i < rec->n; ++i) {
1082 rec->p[i] = isl_upoly_mul(rec->p[i],
1083 isl_upoly_copy(up2));
1084 if (!rec->p[i])
1085 goto error;
1087 isl_upoly_free(up2);
1088 return up1;
1091 if (isl_upoly_is_cst(up1))
1092 return isl_upoly_mul_cst(up1, up2);
1094 return isl_upoly_mul_rec(up1, up2);
1095 error:
1096 isl_upoly_free(up1);
1097 isl_upoly_free(up2);
1098 return NULL;
1101 __isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
1102 unsigned power)
1104 struct isl_upoly *res;
1106 if (!up)
1107 return NULL;
1108 if (power == 1)
1109 return up;
1111 if (power % 2)
1112 res = isl_upoly_copy(up);
1113 else
1114 res = isl_upoly_one(up->ctx);
1116 while (power >>= 1) {
1117 up = isl_upoly_mul(up, isl_upoly_copy(up));
1118 if (power % 2)
1119 res = isl_upoly_mul(res, isl_upoly_copy(up));
1122 isl_upoly_free(up);
1123 return res;
1126 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
1127 unsigned n_div, __isl_take struct isl_upoly *up)
1129 struct isl_qpolynomial *qp = NULL;
1130 unsigned total;
1132 if (!dim || !up)
1133 goto error;
1135 if (!isl_space_is_set(dim))
1136 isl_die(isl_space_get_ctx(dim), isl_error_invalid,
1137 "domain of polynomial should be a set", goto error);
1139 total = isl_space_dim(dim, isl_dim_all);
1141 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
1142 if (!qp)
1143 goto error;
1145 qp->ref = 1;
1146 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
1147 if (!qp->div)
1148 goto error;
1150 qp->dim = dim;
1151 qp->upoly = up;
1153 return qp;
1154 error:
1155 isl_space_free(dim);
1156 isl_upoly_free(up);
1157 isl_qpolynomial_free(qp);
1158 return NULL;
1161 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1163 if (!qp)
1164 return NULL;
1166 qp->ref++;
1167 return qp;
1170 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1172 struct isl_qpolynomial *dup;
1174 if (!qp)
1175 return NULL;
1177 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1178 isl_upoly_copy(qp->upoly));
1179 if (!dup)
1180 return NULL;
1181 isl_mat_free(dup->div);
1182 dup->div = isl_mat_copy(qp->div);
1183 if (!dup->div)
1184 goto error;
1186 return dup;
1187 error:
1188 isl_qpolynomial_free(dup);
1189 return NULL;
1192 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1194 if (!qp)
1195 return NULL;
1197 if (qp->ref == 1)
1198 return qp;
1199 qp->ref--;
1200 return isl_qpolynomial_dup(qp);
1203 __isl_null isl_qpolynomial *isl_qpolynomial_free(
1204 __isl_take isl_qpolynomial *qp)
1206 if (!qp)
1207 return NULL;
1209 if (--qp->ref > 0)
1210 return NULL;
1212 isl_space_free(qp->dim);
1213 isl_mat_free(qp->div);
1214 isl_upoly_free(qp->upoly);
1216 free(qp);
1217 return NULL;
1220 __isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
1222 int i;
1223 struct isl_upoly_rec *rec;
1224 struct isl_upoly_cst *cst;
1226 rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1227 if (!rec)
1228 return NULL;
1229 for (i = 0; i < 1 + power; ++i) {
1230 rec->p[i] = isl_upoly_zero(ctx);
1231 if (!rec->p[i])
1232 goto error;
1233 rec->n++;
1235 cst = isl_upoly_as_cst(rec->p[power]);
1236 isl_int_set_si(cst->n, 1);
1238 return &rec->up;
1239 error:
1240 isl_upoly_free(&rec->up);
1241 return NULL;
1244 /* r array maps original positions to new positions.
1246 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
1247 int *r)
1249 int i;
1250 struct isl_upoly_rec *rec;
1251 struct isl_upoly *base;
1252 struct isl_upoly *res;
1254 if (isl_upoly_is_cst(up))
1255 return up;
1257 rec = isl_upoly_as_rec(up);
1258 if (!rec)
1259 goto error;
1261 isl_assert(up->ctx, rec->n >= 1, goto error);
1263 base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
1264 res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
1266 for (i = rec->n - 2; i >= 0; --i) {
1267 res = isl_upoly_mul(res, isl_upoly_copy(base));
1268 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
1271 isl_upoly_free(base);
1272 isl_upoly_free(up);
1274 return res;
1275 error:
1276 isl_upoly_free(up);
1277 return NULL;
1280 static isl_bool compatible_divs(__isl_keep isl_mat *div1,
1281 __isl_keep isl_mat *div2)
1283 int n_row, n_col;
1284 isl_bool equal;
1286 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1287 div1->n_col >= div2->n_col,
1288 return isl_bool_error);
1290 if (div1->n_row == div2->n_row)
1291 return isl_mat_is_equal(div1, div2);
1293 n_row = div1->n_row;
1294 n_col = div1->n_col;
1295 div1->n_row = div2->n_row;
1296 div1->n_col = div2->n_col;
1298 equal = isl_mat_is_equal(div1, div2);
1300 div1->n_row = n_row;
1301 div1->n_col = n_col;
1303 return equal;
1306 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1308 int li, lj;
1310 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1311 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1313 if (li != lj)
1314 return li - lj;
1316 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1319 struct isl_div_sort_info {
1320 isl_mat *div;
1321 int row;
1324 static int div_sort_cmp(const void *p1, const void *p2)
1326 const struct isl_div_sort_info *i1, *i2;
1327 i1 = (const struct isl_div_sort_info *) p1;
1328 i2 = (const struct isl_div_sort_info *) p2;
1330 return cmp_row(i1->div, i1->row, i2->row);
1333 /* Sort divs and remove duplicates.
1335 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1337 int i;
1338 int skip;
1339 int len;
1340 struct isl_div_sort_info *array = NULL;
1341 int *pos = NULL, *at = NULL;
1342 int *reordering = NULL;
1343 unsigned div_pos;
1345 if (!qp)
1346 return NULL;
1347 if (qp->div->n_row <= 1)
1348 return qp;
1350 div_pos = isl_space_dim(qp->dim, isl_dim_all);
1352 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1353 qp->div->n_row);
1354 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1355 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1356 len = qp->div->n_col - 2;
1357 reordering = isl_alloc_array(qp->div->ctx, int, len);
1358 if (!array || !pos || !at || !reordering)
1359 goto error;
1361 for (i = 0; i < qp->div->n_row; ++i) {
1362 array[i].div = qp->div;
1363 array[i].row = i;
1364 pos[i] = i;
1365 at[i] = i;
1368 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1369 div_sort_cmp);
1371 for (i = 0; i < div_pos; ++i)
1372 reordering[i] = i;
1374 for (i = 0; i < qp->div->n_row; ++i) {
1375 if (pos[array[i].row] == i)
1376 continue;
1377 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1378 pos[at[i]] = pos[array[i].row];
1379 at[pos[array[i].row]] = at[i];
1380 at[i] = array[i].row;
1381 pos[array[i].row] = i;
1384 skip = 0;
1385 for (i = 0; i < len - div_pos; ++i) {
1386 if (i > 0 &&
1387 isl_seq_eq(qp->div->row[i - skip - 1],
1388 qp->div->row[i - skip], qp->div->n_col)) {
1389 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1390 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1391 2 + div_pos + i - skip);
1392 qp->div = isl_mat_drop_cols(qp->div,
1393 2 + div_pos + i - skip, 1);
1394 skip++;
1396 reordering[div_pos + array[i].row] = div_pos + i - skip;
1399 qp->upoly = reorder(qp->upoly, reordering);
1401 if (!qp->upoly || !qp->div)
1402 goto error;
1404 free(at);
1405 free(pos);
1406 free(array);
1407 free(reordering);
1409 return qp;
1410 error:
1411 free(at);
1412 free(pos);
1413 free(array);
1414 free(reordering);
1415 isl_qpolynomial_free(qp);
1416 return NULL;
1419 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1420 int *exp, int first)
1422 int i;
1423 struct isl_upoly_rec *rec;
1425 if (isl_upoly_is_cst(up))
1426 return up;
1428 if (up->var < first)
1429 return up;
1431 if (exp[up->var - first] == up->var - first)
1432 return up;
1434 up = isl_upoly_cow(up);
1435 if (!up)
1436 goto error;
1438 up->var = exp[up->var - first] + first;
1440 rec = isl_upoly_as_rec(up);
1441 if (!rec)
1442 goto error;
1444 for (i = 0; i < rec->n; ++i) {
1445 rec->p[i] = expand(rec->p[i], exp, first);
1446 if (!rec->p[i])
1447 goto error;
1450 return up;
1451 error:
1452 isl_upoly_free(up);
1453 return NULL;
1456 static __isl_give isl_qpolynomial *with_merged_divs(
1457 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1458 __isl_take isl_qpolynomial *qp2),
1459 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1461 int *exp1 = NULL;
1462 int *exp2 = NULL;
1463 isl_mat *div = NULL;
1464 int n_div1, n_div2;
1466 qp1 = isl_qpolynomial_cow(qp1);
1467 qp2 = isl_qpolynomial_cow(qp2);
1469 if (!qp1 || !qp2)
1470 goto error;
1472 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1473 qp1->div->n_col >= qp2->div->n_col, goto error);
1475 n_div1 = qp1->div->n_row;
1476 n_div2 = qp2->div->n_row;
1477 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1);
1478 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2);
1479 if ((n_div1 && !exp1) || (n_div2 && !exp2))
1480 goto error;
1482 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1483 if (!div)
1484 goto error;
1486 isl_mat_free(qp1->div);
1487 qp1->div = isl_mat_copy(div);
1488 isl_mat_free(qp2->div);
1489 qp2->div = isl_mat_copy(div);
1491 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1492 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1494 if (!qp1->upoly || !qp2->upoly)
1495 goto error;
1497 isl_mat_free(div);
1498 free(exp1);
1499 free(exp2);
1501 return fn(qp1, qp2);
1502 error:
1503 isl_mat_free(div);
1504 free(exp1);
1505 free(exp2);
1506 isl_qpolynomial_free(qp1);
1507 isl_qpolynomial_free(qp2);
1508 return NULL;
1511 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1512 __isl_take isl_qpolynomial *qp2)
1514 isl_bool compatible;
1516 qp1 = isl_qpolynomial_cow(qp1);
1518 if (!qp1 || !qp2)
1519 goto error;
1521 if (qp1->div->n_row < qp2->div->n_row)
1522 return isl_qpolynomial_add(qp2, qp1);
1524 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1525 compatible = compatible_divs(qp1->div, qp2->div);
1526 if (compatible < 0)
1527 goto error;
1528 if (!compatible)
1529 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1531 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1532 if (!qp1->upoly)
1533 goto error;
1535 isl_qpolynomial_free(qp2);
1537 return qp1;
1538 error:
1539 isl_qpolynomial_free(qp1);
1540 isl_qpolynomial_free(qp2);
1541 return NULL;
1544 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1545 __isl_keep isl_set *dom,
1546 __isl_take isl_qpolynomial *qp1,
1547 __isl_take isl_qpolynomial *qp2)
1549 qp1 = isl_qpolynomial_add(qp1, qp2);
1550 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1551 return qp1;
1554 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1555 __isl_take isl_qpolynomial *qp2)
1557 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1560 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1561 __isl_take isl_qpolynomial *qp, isl_int v)
1563 if (isl_int_is_zero(v))
1564 return qp;
1566 qp = isl_qpolynomial_cow(qp);
1567 if (!qp)
1568 return NULL;
1570 qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1571 if (!qp->upoly)
1572 goto error;
1574 return qp;
1575 error:
1576 isl_qpolynomial_free(qp);
1577 return NULL;
1581 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1583 if (!qp)
1584 return NULL;
1586 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1589 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1590 __isl_take isl_qpolynomial *qp, isl_int v)
1592 if (isl_int_is_one(v))
1593 return qp;
1595 if (qp && isl_int_is_zero(v)) {
1596 isl_qpolynomial *zero;
1597 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1598 isl_qpolynomial_free(qp);
1599 return zero;
1602 qp = isl_qpolynomial_cow(qp);
1603 if (!qp)
1604 return NULL;
1606 qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1607 if (!qp->upoly)
1608 goto error;
1610 return qp;
1611 error:
1612 isl_qpolynomial_free(qp);
1613 return NULL;
1616 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1617 __isl_take isl_qpolynomial *qp, isl_int v)
1619 return isl_qpolynomial_mul_isl_int(qp, v);
1622 /* Multiply "qp" by "v".
1624 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1625 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1627 if (!qp || !v)
1628 goto error;
1630 if (!isl_val_is_rat(v))
1631 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1632 "expecting rational factor", goto error);
1634 if (isl_val_is_one(v)) {
1635 isl_val_free(v);
1636 return qp;
1639 if (isl_val_is_zero(v)) {
1640 isl_space *space;
1642 space = isl_qpolynomial_get_domain_space(qp);
1643 isl_qpolynomial_free(qp);
1644 isl_val_free(v);
1645 return isl_qpolynomial_zero_on_domain(space);
1648 qp = isl_qpolynomial_cow(qp);
1649 if (!qp)
1650 goto error;
1652 qp->upoly = isl_upoly_scale_val(qp->upoly, v);
1653 if (!qp->upoly)
1654 qp = isl_qpolynomial_free(qp);
1656 isl_val_free(v);
1657 return qp;
1658 error:
1659 isl_val_free(v);
1660 isl_qpolynomial_free(qp);
1661 return NULL;
1664 /* Divide "qp" by "v".
1666 __isl_give isl_qpolynomial *isl_qpolynomial_scale_down_val(
1667 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1669 if (!qp || !v)
1670 goto error;
1672 if (!isl_val_is_rat(v))
1673 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1674 "expecting rational factor", goto error);
1675 if (isl_val_is_zero(v))
1676 isl_die(isl_val_get_ctx(v), isl_error_invalid,
1677 "cannot scale down by zero", goto error);
1679 return isl_qpolynomial_scale_val(qp, isl_val_inv(v));
1680 error:
1681 isl_val_free(v);
1682 isl_qpolynomial_free(qp);
1683 return NULL;
1686 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1687 __isl_take isl_qpolynomial *qp2)
1689 isl_bool compatible;
1691 qp1 = isl_qpolynomial_cow(qp1);
1693 if (!qp1 || !qp2)
1694 goto error;
1696 if (qp1->div->n_row < qp2->div->n_row)
1697 return isl_qpolynomial_mul(qp2, qp1);
1699 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1700 compatible = compatible_divs(qp1->div, qp2->div);
1701 if (compatible < 0)
1702 goto error;
1703 if (!compatible)
1704 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1706 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1707 if (!qp1->upoly)
1708 goto error;
1710 isl_qpolynomial_free(qp2);
1712 return qp1;
1713 error:
1714 isl_qpolynomial_free(qp1);
1715 isl_qpolynomial_free(qp2);
1716 return NULL;
1719 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1720 unsigned power)
1722 qp = isl_qpolynomial_cow(qp);
1724 if (!qp)
1725 return NULL;
1727 qp->upoly = isl_upoly_pow(qp->upoly, power);
1728 if (!qp->upoly)
1729 goto error;
1731 return qp;
1732 error:
1733 isl_qpolynomial_free(qp);
1734 return NULL;
1737 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1738 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1740 int i;
1742 if (power == 1)
1743 return pwqp;
1745 pwqp = isl_pw_qpolynomial_cow(pwqp);
1746 if (!pwqp)
1747 return NULL;
1749 for (i = 0; i < pwqp->n; ++i) {
1750 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1751 if (!pwqp->p[i].qp)
1752 return isl_pw_qpolynomial_free(pwqp);
1755 return pwqp;
1758 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1759 __isl_take isl_space *dim)
1761 if (!dim)
1762 return NULL;
1763 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1766 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1767 __isl_take isl_space *dim)
1769 if (!dim)
1770 return NULL;
1771 return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1774 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1775 __isl_take isl_space *dim)
1777 if (!dim)
1778 return NULL;
1779 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1782 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1783 __isl_take isl_space *dim)
1785 if (!dim)
1786 return NULL;
1787 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1790 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1791 __isl_take isl_space *dim)
1793 if (!dim)
1794 return NULL;
1795 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1798 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1799 __isl_take isl_space *dim,
1800 isl_int v)
1802 struct isl_qpolynomial *qp;
1803 struct isl_upoly_cst *cst;
1805 if (!dim)
1806 return NULL;
1808 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1809 if (!qp)
1810 return NULL;
1812 cst = isl_upoly_as_cst(qp->upoly);
1813 isl_int_set(cst->n, v);
1815 return qp;
1818 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1819 isl_int *n, isl_int *d)
1821 struct isl_upoly_cst *cst;
1823 if (!qp)
1824 return -1;
1826 if (!isl_upoly_is_cst(qp->upoly))
1827 return 0;
1829 cst = isl_upoly_as_cst(qp->upoly);
1830 if (!cst)
1831 return -1;
1833 if (n)
1834 isl_int_set(*n, cst->n);
1835 if (d)
1836 isl_int_set(*d, cst->d);
1838 return 1;
1841 /* Return the constant term of "up".
1843 static __isl_give isl_val *isl_upoly_get_constant_val(
1844 __isl_keep struct isl_upoly *up)
1846 struct isl_upoly_cst *cst;
1848 if (!up)
1849 return NULL;
1851 while (!isl_upoly_is_cst(up)) {
1852 struct isl_upoly_rec *rec;
1854 rec = isl_upoly_as_rec(up);
1855 if (!rec)
1856 return NULL;
1857 up = rec->p[0];
1860 cst = isl_upoly_as_cst(up);
1861 if (!cst)
1862 return NULL;
1863 return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1866 /* Return the constant term of "qp".
1868 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1869 __isl_keep isl_qpolynomial *qp)
1871 if (!qp)
1872 return NULL;
1874 return isl_upoly_get_constant_val(qp->upoly);
1877 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1879 int is_cst;
1880 struct isl_upoly_rec *rec;
1882 if (!up)
1883 return -1;
1885 if (up->var < 0)
1886 return 1;
1888 rec = isl_upoly_as_rec(up);
1889 if (!rec)
1890 return -1;
1892 if (rec->n > 2)
1893 return 0;
1895 isl_assert(up->ctx, rec->n > 1, return -1);
1897 is_cst = isl_upoly_is_cst(rec->p[1]);
1898 if (is_cst < 0)
1899 return -1;
1900 if (!is_cst)
1901 return 0;
1903 return isl_upoly_is_affine(rec->p[0]);
1906 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1908 if (!qp)
1909 return -1;
1911 if (qp->div->n_row > 0)
1912 return 0;
1914 return isl_upoly_is_affine(qp->upoly);
1917 static void update_coeff(__isl_keep isl_vec *aff,
1918 __isl_keep struct isl_upoly_cst *cst, int pos)
1920 isl_int gcd;
1921 isl_int f;
1923 if (isl_int_is_zero(cst->n))
1924 return;
1926 isl_int_init(gcd);
1927 isl_int_init(f);
1928 isl_int_gcd(gcd, cst->d, aff->el[0]);
1929 isl_int_divexact(f, cst->d, gcd);
1930 isl_int_divexact(gcd, aff->el[0], gcd);
1931 isl_seq_scale(aff->el, aff->el, f, aff->size);
1932 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1933 isl_int_clear(gcd);
1934 isl_int_clear(f);
1937 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1938 __isl_keep isl_vec *aff)
1940 struct isl_upoly_cst *cst;
1941 struct isl_upoly_rec *rec;
1943 if (!up || !aff)
1944 return -1;
1946 if (up->var < 0) {
1947 struct isl_upoly_cst *cst;
1949 cst = isl_upoly_as_cst(up);
1950 if (!cst)
1951 return -1;
1952 update_coeff(aff, cst, 0);
1953 return 0;
1956 rec = isl_upoly_as_rec(up);
1957 if (!rec)
1958 return -1;
1959 isl_assert(up->ctx, rec->n == 2, return -1);
1961 cst = isl_upoly_as_cst(rec->p[1]);
1962 if (!cst)
1963 return -1;
1964 update_coeff(aff, cst, 1 + up->var);
1966 return isl_upoly_update_affine(rec->p[0], aff);
1969 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1970 __isl_keep isl_qpolynomial *qp)
1972 isl_vec *aff;
1973 unsigned d;
1975 if (!qp)
1976 return NULL;
1978 d = isl_space_dim(qp->dim, isl_dim_all);
1979 aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1980 if (!aff)
1981 return NULL;
1983 isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1984 isl_int_set_si(aff->el[0], 1);
1986 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1987 goto error;
1989 return aff;
1990 error:
1991 isl_vec_free(aff);
1992 return NULL;
1995 /* Compare two quasi-polynomials.
1997 * Return -1 if "qp1" is "smaller" than "qp2", 1 if "qp1" is "greater"
1998 * than "qp2" and 0 if they are equal.
2000 int isl_qpolynomial_plain_cmp(__isl_keep isl_qpolynomial *qp1,
2001 __isl_keep isl_qpolynomial *qp2)
2003 int cmp;
2005 if (qp1 == qp2)
2006 return 0;
2007 if (!qp1)
2008 return -1;
2009 if (!qp2)
2010 return 1;
2012 cmp = isl_space_cmp(qp1->dim, qp2->dim);
2013 if (cmp != 0)
2014 return cmp;
2016 cmp = isl_local_cmp(qp1->div, qp2->div);
2017 if (cmp != 0)
2018 return cmp;
2020 return isl_upoly_plain_cmp(qp1->upoly, qp2->upoly);
2023 /* Is "qp1" obviously equal to "qp2"?
2025 * NaN is not equal to anything, not even to another NaN.
2027 isl_bool isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
2028 __isl_keep isl_qpolynomial *qp2)
2030 isl_bool equal;
2032 if (!qp1 || !qp2)
2033 return isl_bool_error;
2035 if (isl_qpolynomial_is_nan(qp1) || isl_qpolynomial_is_nan(qp2))
2036 return isl_bool_false;
2038 equal = isl_space_is_equal(qp1->dim, qp2->dim);
2039 if (equal < 0 || !equal)
2040 return equal;
2042 equal = isl_mat_is_equal(qp1->div, qp2->div);
2043 if (equal < 0 || !equal)
2044 return equal;
2046 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
2049 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
2051 int i;
2052 struct isl_upoly_rec *rec;
2054 if (isl_upoly_is_cst(up)) {
2055 struct isl_upoly_cst *cst;
2056 cst = isl_upoly_as_cst(up);
2057 if (!cst)
2058 return;
2059 isl_int_lcm(*d, *d, cst->d);
2060 return;
2063 rec = isl_upoly_as_rec(up);
2064 if (!rec)
2065 return;
2067 for (i = 0; i < rec->n; ++i)
2068 upoly_update_den(rec->p[i], d);
2071 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
2073 isl_int_set_si(*d, 1);
2074 if (!qp)
2075 return;
2076 upoly_update_den(qp->upoly, d);
2079 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
2080 __isl_take isl_space *dim, int pos, int power)
2082 struct isl_ctx *ctx;
2084 if (!dim)
2085 return NULL;
2087 ctx = dim->ctx;
2089 return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
2092 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
2093 enum isl_dim_type type, unsigned pos)
2095 if (!dim)
2096 return NULL;
2098 isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
2099 isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
2101 if (type == isl_dim_set)
2102 pos += isl_space_dim(dim, isl_dim_param);
2104 return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
2105 error:
2106 isl_space_free(dim);
2107 return NULL;
2110 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
2111 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
2113 int i;
2114 struct isl_upoly_rec *rec;
2115 struct isl_upoly *base, *res;
2117 if (!up)
2118 return NULL;
2120 if (isl_upoly_is_cst(up))
2121 return up;
2123 if (up->var < first)
2124 return up;
2126 rec = isl_upoly_as_rec(up);
2127 if (!rec)
2128 goto error;
2130 isl_assert(up->ctx, rec->n >= 1, goto error);
2132 if (up->var >= first + n)
2133 base = isl_upoly_var_pow(up->ctx, up->var, 1);
2134 else
2135 base = isl_upoly_copy(subs[up->var - first]);
2137 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
2138 for (i = rec->n - 2; i >= 0; --i) {
2139 struct isl_upoly *t;
2140 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
2141 res = isl_upoly_mul(res, isl_upoly_copy(base));
2142 res = isl_upoly_sum(res, t);
2145 isl_upoly_free(base);
2146 isl_upoly_free(up);
2148 return res;
2149 error:
2150 isl_upoly_free(up);
2151 return NULL;
2154 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
2155 isl_int denom, unsigned len)
2157 int i;
2158 struct isl_upoly *up;
2160 isl_assert(ctx, len >= 1, return NULL);
2162 up = isl_upoly_rat_cst(ctx, f[0], denom);
2163 for (i = 0; i < len - 1; ++i) {
2164 struct isl_upoly *t;
2165 struct isl_upoly *c;
2167 if (isl_int_is_zero(f[1 + i]))
2168 continue;
2170 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
2171 t = isl_upoly_var_pow(ctx, i, 1);
2172 t = isl_upoly_mul(c, t);
2173 up = isl_upoly_sum(up, t);
2176 return up;
2179 /* Remove common factor of non-constant terms and denominator.
2181 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2183 isl_ctx *ctx = qp->div->ctx;
2184 unsigned total = qp->div->n_col - 2;
2186 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2187 isl_int_gcd(ctx->normalize_gcd,
2188 ctx->normalize_gcd, qp->div->row[div][0]);
2189 if (isl_int_is_one(ctx->normalize_gcd))
2190 return;
2192 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2193 ctx->normalize_gcd, total);
2194 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2195 ctx->normalize_gcd);
2196 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2197 ctx->normalize_gcd);
2200 /* Replace the integer division identified by "div" by the polynomial "s".
2201 * The integer division is assumed not to appear in the definition
2202 * of any other integer divisions.
2204 static __isl_give isl_qpolynomial *substitute_div(
2205 __isl_take isl_qpolynomial *qp,
2206 int div, __isl_take struct isl_upoly *s)
2208 int i;
2209 int total;
2210 int *reordering;
2212 if (!qp || !s)
2213 goto error;
2215 qp = isl_qpolynomial_cow(qp);
2216 if (!qp)
2217 goto error;
2219 total = isl_space_dim(qp->dim, isl_dim_all);
2220 qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
2221 if (!qp->upoly)
2222 goto error;
2224 reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
2225 if (!reordering)
2226 goto error;
2227 for (i = 0; i < total + div; ++i)
2228 reordering[i] = i;
2229 for (i = total + div + 1; i < total + qp->div->n_row; ++i)
2230 reordering[i] = i - 1;
2231 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2232 qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
2233 qp->upoly = reorder(qp->upoly, reordering);
2234 free(reordering);
2236 if (!qp->upoly || !qp->div)
2237 goto error;
2239 isl_upoly_free(s);
2240 return qp;
2241 error:
2242 isl_qpolynomial_free(qp);
2243 isl_upoly_free(s);
2244 return NULL;
2247 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2248 * divisions because d is equal to 1 by their definition, i.e., e.
2250 static __isl_give isl_qpolynomial *substitute_non_divs(
2251 __isl_take isl_qpolynomial *qp)
2253 int i, j;
2254 int total;
2255 struct isl_upoly *s;
2257 if (!qp)
2258 return NULL;
2260 total = isl_space_dim(qp->dim, isl_dim_all);
2261 for (i = 0; qp && i < qp->div->n_row; ++i) {
2262 if (!isl_int_is_one(qp->div->row[i][0]))
2263 continue;
2264 for (j = i + 1; j < qp->div->n_row; ++j) {
2265 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2266 continue;
2267 isl_seq_combine(qp->div->row[j] + 1,
2268 qp->div->ctx->one, qp->div->row[j] + 1,
2269 qp->div->row[j][2 + total + i],
2270 qp->div->row[i] + 1, 1 + total + i);
2271 isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2272 normalize_div(qp, j);
2274 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2275 qp->div->row[i][0], qp->div->n_col - 1);
2276 qp = substitute_div(qp, i, s);
2277 --i;
2280 return qp;
2283 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2284 * with d the denominator. When replacing the coefficient e of x by
2285 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2286 * inside the division, so we need to add floor(e/d) * x outside.
2287 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2288 * to adjust the coefficient of x in each later div that depends on the
2289 * current div "div" and also in the affine expressions in the rows of "mat"
2290 * (if they too depend on "div").
2292 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2293 __isl_keep isl_mat **mat)
2295 int i, j;
2296 isl_int v;
2297 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2299 isl_int_init(v);
2300 for (i = 0; i < 1 + total + div; ++i) {
2301 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2302 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2303 continue;
2304 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2305 isl_int_fdiv_r(qp->div->row[div][1 + i],
2306 qp->div->row[div][1 + i], qp->div->row[div][0]);
2307 *mat = isl_mat_col_addmul(*mat, i, v, 1 + total + div);
2308 for (j = div + 1; j < qp->div->n_row; ++j) {
2309 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2310 continue;
2311 isl_int_addmul(qp->div->row[j][1 + i],
2312 v, qp->div->row[j][2 + total + div]);
2315 isl_int_clear(v);
2318 /* Check if the last non-zero coefficient is bigger that half of the
2319 * denominator. If so, we will invert the div to further reduce the number
2320 * of distinct divs that may appear.
2321 * If the last non-zero coefficient is exactly half the denominator,
2322 * then we continue looking for earlier coefficients that are bigger
2323 * than half the denominator.
2325 static int needs_invert(__isl_keep isl_mat *div, int row)
2327 int i;
2328 int cmp;
2330 for (i = div->n_col - 1; i >= 1; --i) {
2331 if (isl_int_is_zero(div->row[row][i]))
2332 continue;
2333 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2334 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2335 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2336 if (cmp)
2337 return cmp > 0;
2338 if (i == 1)
2339 return 1;
2342 return 0;
2345 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2346 * We only invert the coefficients of e (and the coefficient of q in
2347 * later divs and in the rows of "mat"). After calling this function, the
2348 * coefficients of e should be reduced again.
2350 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2351 __isl_keep isl_mat **mat)
2353 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2355 isl_seq_neg(qp->div->row[div] + 1,
2356 qp->div->row[div] + 1, qp->div->n_col - 1);
2357 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2358 isl_int_add(qp->div->row[div][1],
2359 qp->div->row[div][1], qp->div->row[div][0]);
2360 *mat = isl_mat_col_neg(*mat, 1 + total + div);
2361 isl_mat_col_mul(qp->div, 2 + total + div,
2362 qp->div->ctx->negone, 2 + total + div);
2365 /* Reduce all divs of "qp" to have coefficients
2366 * in the interval [0, d-1], with d the denominator and such that the
2367 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2368 * The modifications to the integer divisions need to be reflected
2369 * in the factors of the polynomial that refer to the original
2370 * integer divisions. To this end, the modifications are collected
2371 * as a set of affine expressions and then plugged into the polynomial.
2373 * After the reduction, some divs may have become redundant or identical,
2374 * so we call substitute_non_divs and sort_divs. If these functions
2375 * eliminate divs or merge two or more divs into one, the coefficients
2376 * of the enclosing divs may have to be reduced again, so we call
2377 * ourselves recursively if the number of divs decreases.
2379 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2381 int i;
2382 isl_ctx *ctx;
2383 isl_mat *mat;
2384 struct isl_upoly **s;
2385 unsigned o_div, n_div, total;
2387 if (!qp)
2388 return NULL;
2390 total = isl_qpolynomial_domain_dim(qp, isl_dim_all);
2391 n_div = isl_qpolynomial_domain_dim(qp, isl_dim_div);
2392 o_div = isl_qpolynomial_domain_offset(qp, isl_dim_div);
2393 ctx = isl_qpolynomial_get_ctx(qp);
2394 mat = isl_mat_zero(ctx, n_div, 1 + total);
2396 for (i = 0; i < n_div; ++i)
2397 mat = isl_mat_set_element_si(mat, i, o_div + i, 1);
2399 for (i = 0; i < qp->div->n_row; ++i) {
2400 normalize_div(qp, i);
2401 reduce_div(qp, i, &mat);
2402 if (needs_invert(qp->div, i)) {
2403 invert_div(qp, i, &mat);
2404 reduce_div(qp, i, &mat);
2407 if (!mat)
2408 goto error;
2410 s = isl_alloc_array(ctx, struct isl_upoly *, n_div);
2411 if (n_div && !s)
2412 goto error;
2413 for (i = 0; i < n_div; ++i)
2414 s[i] = isl_upoly_from_affine(ctx, mat->row[i], ctx->one,
2415 1 + total);
2416 qp->upoly = isl_upoly_subs(qp->upoly, o_div - 1, n_div, s);
2417 for (i = 0; i < n_div; ++i)
2418 isl_upoly_free(s[i]);
2419 free(s);
2420 if (!qp->upoly)
2421 goto error;
2423 isl_mat_free(mat);
2425 qp = substitute_non_divs(qp);
2426 qp = sort_divs(qp);
2427 if (qp && isl_qpolynomial_domain_dim(qp, isl_dim_div) < n_div)
2428 return reduce_divs(qp);
2430 return qp;
2431 error:
2432 isl_qpolynomial_free(qp);
2433 isl_mat_free(mat);
2434 return NULL;
2437 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2438 __isl_take isl_space *dim, const isl_int n, const isl_int d)
2440 struct isl_qpolynomial *qp;
2441 struct isl_upoly_cst *cst;
2443 if (!dim)
2444 return NULL;
2446 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2447 if (!qp)
2448 return NULL;
2450 cst = isl_upoly_as_cst(qp->upoly);
2451 isl_int_set(cst->n, n);
2452 isl_int_set(cst->d, d);
2454 return qp;
2457 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2459 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2460 __isl_take isl_space *domain, __isl_take isl_val *val)
2462 isl_qpolynomial *qp;
2463 struct isl_upoly_cst *cst;
2465 if (!domain || !val)
2466 goto error;
2468 qp = isl_qpolynomial_alloc(isl_space_copy(domain), 0,
2469 isl_upoly_zero(domain->ctx));
2470 if (!qp)
2471 goto error;
2473 cst = isl_upoly_as_cst(qp->upoly);
2474 isl_int_set(cst->n, val->n);
2475 isl_int_set(cst->d, val->d);
2477 isl_space_free(domain);
2478 isl_val_free(val);
2479 return qp;
2480 error:
2481 isl_space_free(domain);
2482 isl_val_free(val);
2483 return NULL;
2486 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2488 struct isl_upoly_rec *rec;
2489 int i;
2491 if (!up)
2492 return -1;
2494 if (isl_upoly_is_cst(up))
2495 return 0;
2497 if (up->var < d)
2498 active[up->var] = 1;
2500 rec = isl_upoly_as_rec(up);
2501 for (i = 0; i < rec->n; ++i)
2502 if (up_set_active(rec->p[i], active, d) < 0)
2503 return -1;
2505 return 0;
2508 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2510 int i, j;
2511 int d = isl_space_dim(qp->dim, isl_dim_all);
2513 if (!qp || !active)
2514 return -1;
2516 for (i = 0; i < d; ++i)
2517 for (j = 0; j < qp->div->n_row; ++j) {
2518 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2519 continue;
2520 active[i] = 1;
2521 break;
2524 return up_set_active(qp->upoly, active, d);
2527 isl_bool isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2528 enum isl_dim_type type, unsigned first, unsigned n)
2530 int i;
2531 int *active = NULL;
2532 isl_bool involves = isl_bool_false;
2534 if (!qp)
2535 return isl_bool_error;
2536 if (n == 0)
2537 return isl_bool_false;
2539 isl_assert(qp->dim->ctx,
2540 first + n <= isl_qpolynomial_dim(qp, type),
2541 return isl_bool_error);
2542 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2543 type == isl_dim_in, return isl_bool_error);
2545 active = isl_calloc_array(qp->dim->ctx, int,
2546 isl_space_dim(qp->dim, isl_dim_all));
2547 if (set_active(qp, active) < 0)
2548 goto error;
2550 if (type == isl_dim_in)
2551 first += isl_space_dim(qp->dim, isl_dim_param);
2552 for (i = 0; i < n; ++i)
2553 if (active[first + i]) {
2554 involves = isl_bool_true;
2555 break;
2558 free(active);
2560 return involves;
2561 error:
2562 free(active);
2563 return isl_bool_error;
2566 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2567 * of the divs that do appear in the quasi-polynomial.
2569 static __isl_give isl_qpolynomial *remove_redundant_divs(
2570 __isl_take isl_qpolynomial *qp)
2572 int i, j;
2573 int d;
2574 int len;
2575 int skip;
2576 int *active = NULL;
2577 int *reordering = NULL;
2578 int redundant = 0;
2579 int n_div;
2580 isl_ctx *ctx;
2582 if (!qp)
2583 return NULL;
2584 if (qp->div->n_row == 0)
2585 return qp;
2587 d = isl_space_dim(qp->dim, isl_dim_all);
2588 len = qp->div->n_col - 2;
2589 ctx = isl_qpolynomial_get_ctx(qp);
2590 active = isl_calloc_array(ctx, int, len);
2591 if (!active)
2592 goto error;
2594 if (up_set_active(qp->upoly, active, len) < 0)
2595 goto error;
2597 for (i = qp->div->n_row - 1; i >= 0; --i) {
2598 if (!active[d + i]) {
2599 redundant = 1;
2600 continue;
2602 for (j = 0; j < i; ++j) {
2603 if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2604 continue;
2605 active[d + j] = 1;
2606 break;
2610 if (!redundant) {
2611 free(active);
2612 return qp;
2615 reordering = isl_alloc_array(qp->div->ctx, int, len);
2616 if (!reordering)
2617 goto error;
2619 for (i = 0; i < d; ++i)
2620 reordering[i] = i;
2622 skip = 0;
2623 n_div = qp->div->n_row;
2624 for (i = 0; i < n_div; ++i) {
2625 if (!active[d + i]) {
2626 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2627 qp->div = isl_mat_drop_cols(qp->div,
2628 2 + d + i - skip, 1);
2629 skip++;
2631 reordering[d + i] = d + i - skip;
2634 qp->upoly = reorder(qp->upoly, reordering);
2636 if (!qp->upoly || !qp->div)
2637 goto error;
2639 free(active);
2640 free(reordering);
2642 return qp;
2643 error:
2644 free(active);
2645 free(reordering);
2646 isl_qpolynomial_free(qp);
2647 return NULL;
2650 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2651 unsigned first, unsigned n)
2653 int i;
2654 struct isl_upoly_rec *rec;
2656 if (!up)
2657 return NULL;
2658 if (n == 0 || up->var < 0 || up->var < first)
2659 return up;
2660 if (up->var < first + n) {
2661 up = replace_by_constant_term(up);
2662 return isl_upoly_drop(up, first, n);
2664 up = isl_upoly_cow(up);
2665 if (!up)
2666 return NULL;
2667 up->var -= n;
2668 rec = isl_upoly_as_rec(up);
2669 if (!rec)
2670 goto error;
2672 for (i = 0; i < rec->n; ++i) {
2673 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2674 if (!rec->p[i])
2675 goto error;
2678 return up;
2679 error:
2680 isl_upoly_free(up);
2681 return NULL;
2684 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2685 __isl_take isl_qpolynomial *qp,
2686 enum isl_dim_type type, unsigned pos, const char *s)
2688 qp = isl_qpolynomial_cow(qp);
2689 if (!qp)
2690 return NULL;
2691 if (type == isl_dim_out)
2692 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2693 "cannot set name of output/set dimension",
2694 return isl_qpolynomial_free(qp));
2695 if (type == isl_dim_in)
2696 type = isl_dim_set;
2697 qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2698 if (!qp->dim)
2699 goto error;
2700 return qp;
2701 error:
2702 isl_qpolynomial_free(qp);
2703 return NULL;
2706 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2707 __isl_take isl_qpolynomial *qp,
2708 enum isl_dim_type type, unsigned first, unsigned n)
2710 if (!qp)
2711 return NULL;
2712 if (type == isl_dim_out)
2713 isl_die(qp->dim->ctx, isl_error_invalid,
2714 "cannot drop output/set dimension",
2715 goto error);
2716 if (type == isl_dim_in)
2717 type = isl_dim_set;
2718 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2719 return qp;
2721 qp = isl_qpolynomial_cow(qp);
2722 if (!qp)
2723 return NULL;
2725 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2726 goto error);
2727 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2728 type == isl_dim_set, goto error);
2730 qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2731 if (!qp->dim)
2732 goto error;
2734 if (type == isl_dim_set)
2735 first += isl_space_dim(qp->dim, isl_dim_param);
2737 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2738 if (!qp->div)
2739 goto error;
2741 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2742 if (!qp->upoly)
2743 goto error;
2745 return qp;
2746 error:
2747 isl_qpolynomial_free(qp);
2748 return NULL;
2751 /* Project the domain of the quasi-polynomial onto its parameter space.
2752 * The quasi-polynomial may not involve any of the domain dimensions.
2754 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2755 __isl_take isl_qpolynomial *qp)
2757 isl_space *space;
2758 unsigned n;
2759 int involves;
2761 n = isl_qpolynomial_dim(qp, isl_dim_in);
2762 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2763 if (involves < 0)
2764 return isl_qpolynomial_free(qp);
2765 if (involves)
2766 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2767 "polynomial involves some of the domain dimensions",
2768 return isl_qpolynomial_free(qp));
2769 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2770 space = isl_qpolynomial_get_domain_space(qp);
2771 space = isl_space_params(space);
2772 qp = isl_qpolynomial_reset_domain_space(qp, space);
2773 return qp;
2776 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2777 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2779 int i, j, k;
2780 isl_int denom;
2781 unsigned total;
2782 unsigned n_div;
2783 struct isl_upoly *up;
2785 if (!eq)
2786 goto error;
2787 if (eq->n_eq == 0) {
2788 isl_basic_set_free(eq);
2789 return qp;
2792 qp = isl_qpolynomial_cow(qp);
2793 if (!qp)
2794 goto error;
2795 qp->div = isl_mat_cow(qp->div);
2796 if (!qp->div)
2797 goto error;
2799 total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2800 n_div = eq->n_div;
2801 isl_int_init(denom);
2802 for (i = 0; i < eq->n_eq; ++i) {
2803 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2804 if (j < 0 || j == 0 || j >= total)
2805 continue;
2807 for (k = 0; k < qp->div->n_row; ++k) {
2808 if (isl_int_is_zero(qp->div->row[k][1 + j]))
2809 continue;
2810 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2811 &qp->div->row[k][0]);
2812 normalize_div(qp, k);
2815 if (isl_int_is_pos(eq->eq[i][j]))
2816 isl_seq_neg(eq->eq[i], eq->eq[i], total);
2817 isl_int_abs(denom, eq->eq[i][j]);
2818 isl_int_set_si(eq->eq[i][j], 0);
2820 up = isl_upoly_from_affine(qp->dim->ctx,
2821 eq->eq[i], denom, total);
2822 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2823 isl_upoly_free(up);
2825 isl_int_clear(denom);
2827 if (!qp->upoly)
2828 goto error;
2830 isl_basic_set_free(eq);
2832 qp = substitute_non_divs(qp);
2833 qp = sort_divs(qp);
2835 return qp;
2836 error:
2837 isl_basic_set_free(eq);
2838 isl_qpolynomial_free(qp);
2839 return NULL;
2842 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2844 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2845 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2847 if (!qp || !eq)
2848 goto error;
2849 if (qp->div->n_row > 0)
2850 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2851 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2852 error:
2853 isl_basic_set_free(eq);
2854 isl_qpolynomial_free(qp);
2855 return NULL;
2858 static __isl_give isl_basic_set *add_div_constraints(
2859 __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2861 int i;
2862 unsigned total;
2864 if (!bset || !div)
2865 goto error;
2867 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2868 if (!bset)
2869 goto error;
2870 total = isl_basic_set_total_dim(bset);
2871 for (i = 0; i < div->n_row; ++i)
2872 if (isl_basic_set_add_div_constraints_var(bset,
2873 total - div->n_row + i, div->row[i]) < 0)
2874 goto error;
2876 isl_mat_free(div);
2877 return bset;
2878 error:
2879 isl_mat_free(div);
2880 isl_basic_set_free(bset);
2881 return NULL;
2884 /* Look for equalities among the variables shared by context and qp
2885 * and the integer divisions of qp, if any.
2886 * The equalities are then used to eliminate variables and/or integer
2887 * divisions from qp.
2889 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2890 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2892 isl_basic_set *aff;
2894 if (!qp)
2895 goto error;
2896 if (qp->div->n_row > 0) {
2897 isl_basic_set *bset;
2898 context = isl_set_add_dims(context, isl_dim_set,
2899 qp->div->n_row);
2900 bset = isl_basic_set_universe(isl_set_get_space(context));
2901 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2902 context = isl_set_intersect(context,
2903 isl_set_from_basic_set(bset));
2906 aff = isl_set_affine_hull(context);
2907 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2908 error:
2909 isl_qpolynomial_free(qp);
2910 isl_set_free(context);
2911 return NULL;
2914 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2915 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2917 isl_space *space = isl_qpolynomial_get_domain_space(qp);
2918 isl_set *dom_context = isl_set_universe(space);
2919 dom_context = isl_set_intersect_params(dom_context, context);
2920 return isl_qpolynomial_gist(qp, dom_context);
2923 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2924 __isl_take isl_qpolynomial *qp)
2926 isl_set *dom;
2928 if (!qp)
2929 return NULL;
2930 if (isl_qpolynomial_is_zero(qp)) {
2931 isl_space *dim = isl_qpolynomial_get_space(qp);
2932 isl_qpolynomial_free(qp);
2933 return isl_pw_qpolynomial_zero(dim);
2936 dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2937 return isl_pw_qpolynomial_alloc(dom, qp);
2940 #define isl_qpolynomial_involves_nan isl_qpolynomial_is_nan
2942 #undef PW
2943 #define PW isl_pw_qpolynomial
2944 #undef EL
2945 #define EL isl_qpolynomial
2946 #undef EL_IS_ZERO
2947 #define EL_IS_ZERO is_zero
2948 #undef ZERO
2949 #define ZERO zero
2950 #undef IS_ZERO
2951 #define IS_ZERO is_zero
2952 #undef FIELD
2953 #define FIELD qp
2954 #undef DEFAULT_IS_ZERO
2955 #define DEFAULT_IS_ZERO 1
2957 #define NO_PULLBACK
2959 #include <isl_pw_templ.c>
2960 #include <isl_pw_eval.c>
2962 #undef UNION
2963 #define UNION isl_union_pw_qpolynomial
2964 #undef PART
2965 #define PART isl_pw_qpolynomial
2966 #undef PARTS
2967 #define PARTS pw_qpolynomial
2969 #include <isl_union_single.c>
2970 #include <isl_union_eval.c>
2971 #include <isl_union_neg.c>
2973 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2975 if (!pwqp)
2976 return -1;
2978 if (pwqp->n != -1)
2979 return 0;
2981 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2982 return 0;
2984 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2987 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2988 __isl_take isl_pw_qpolynomial *pwqp1,
2989 __isl_take isl_pw_qpolynomial *pwqp2)
2991 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2994 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2995 __isl_take isl_pw_qpolynomial *pwqp1,
2996 __isl_take isl_pw_qpolynomial *pwqp2)
2998 int i, j, n;
2999 struct isl_pw_qpolynomial *res;
3001 if (!pwqp1 || !pwqp2)
3002 goto error;
3004 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
3005 goto error);
3007 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
3008 isl_pw_qpolynomial_free(pwqp2);
3009 return pwqp1;
3012 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
3013 isl_pw_qpolynomial_free(pwqp1);
3014 return pwqp2;
3017 if (isl_pw_qpolynomial_is_one(pwqp1)) {
3018 isl_pw_qpolynomial_free(pwqp1);
3019 return pwqp2;
3022 if (isl_pw_qpolynomial_is_one(pwqp2)) {
3023 isl_pw_qpolynomial_free(pwqp2);
3024 return pwqp1;
3027 n = pwqp1->n * pwqp2->n;
3028 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
3030 for (i = 0; i < pwqp1->n; ++i) {
3031 for (j = 0; j < pwqp2->n; ++j) {
3032 struct isl_set *common;
3033 struct isl_qpolynomial *prod;
3034 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
3035 isl_set_copy(pwqp2->p[j].set));
3036 if (isl_set_plain_is_empty(common)) {
3037 isl_set_free(common);
3038 continue;
3041 prod = isl_qpolynomial_mul(
3042 isl_qpolynomial_copy(pwqp1->p[i].qp),
3043 isl_qpolynomial_copy(pwqp2->p[j].qp));
3045 res = isl_pw_qpolynomial_add_piece(res, common, prod);
3049 isl_pw_qpolynomial_free(pwqp1);
3050 isl_pw_qpolynomial_free(pwqp2);
3052 return res;
3053 error:
3054 isl_pw_qpolynomial_free(pwqp1);
3055 isl_pw_qpolynomial_free(pwqp2);
3056 return NULL;
3059 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
3060 __isl_take isl_vec *vec)
3062 int i;
3063 struct isl_upoly_rec *rec;
3064 isl_val *res;
3065 isl_val *base;
3067 if (isl_upoly_is_cst(up)) {
3068 isl_vec_free(vec);
3069 res = isl_upoly_get_constant_val(up);
3070 isl_upoly_free(up);
3071 return res;
3074 rec = isl_upoly_as_rec(up);
3075 if (!rec || !vec)
3076 goto error;
3078 isl_assert(up->ctx, rec->n >= 1, goto error);
3080 base = isl_val_rat_from_isl_int(up->ctx,
3081 vec->el[1 + up->var], vec->el[0]);
3083 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
3084 isl_vec_copy(vec));
3086 for (i = rec->n - 2; i >= 0; --i) {
3087 res = isl_val_mul(res, isl_val_copy(base));
3088 res = isl_val_add(res,
3089 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
3090 isl_vec_copy(vec)));
3093 isl_val_free(base);
3094 isl_upoly_free(up);
3095 isl_vec_free(vec);
3096 return res;
3097 error:
3098 isl_upoly_free(up);
3099 isl_vec_free(vec);
3100 return NULL;
3103 /* Evaluate "qp" in the void point "pnt".
3104 * In particular, return the value NaN.
3106 static __isl_give isl_val *eval_void(__isl_take isl_qpolynomial *qp,
3107 __isl_take isl_point *pnt)
3109 isl_ctx *ctx;
3111 ctx = isl_point_get_ctx(pnt);
3112 isl_qpolynomial_free(qp);
3113 isl_point_free(pnt);
3114 return isl_val_nan(ctx);
3117 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
3118 __isl_take isl_point *pnt)
3120 isl_bool is_void;
3121 isl_vec *ext;
3122 isl_val *v;
3124 if (!qp || !pnt)
3125 goto error;
3126 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
3127 is_void = isl_point_is_void(pnt);
3128 if (is_void < 0)
3129 goto error;
3130 if (is_void)
3131 return eval_void(qp, pnt);
3133 ext = isl_local_extend_point_vec(qp->div, isl_vec_copy(pnt->vec));
3135 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
3137 isl_qpolynomial_free(qp);
3138 isl_point_free(pnt);
3140 return v;
3141 error:
3142 isl_qpolynomial_free(qp);
3143 isl_point_free(pnt);
3144 return NULL;
3147 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
3148 __isl_keep struct isl_upoly_cst *cst2)
3150 int cmp;
3151 isl_int t;
3152 isl_int_init(t);
3153 isl_int_mul(t, cst1->n, cst2->d);
3154 isl_int_submul(t, cst2->n, cst1->d);
3155 cmp = isl_int_sgn(t);
3156 isl_int_clear(t);
3157 return cmp;
3160 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
3161 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
3162 unsigned first, unsigned n)
3164 unsigned total;
3165 unsigned g_pos;
3166 int *exp;
3168 if (!qp)
3169 return NULL;
3170 if (type == isl_dim_out)
3171 isl_die(qp->div->ctx, isl_error_invalid,
3172 "cannot insert output/set dimensions",
3173 goto error);
3174 if (type == isl_dim_in)
3175 type = isl_dim_set;
3176 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3177 return qp;
3179 qp = isl_qpolynomial_cow(qp);
3180 if (!qp)
3181 return NULL;
3183 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
3184 goto error);
3186 g_pos = pos(qp->dim, type) + first;
3188 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3189 if (!qp->div)
3190 goto error;
3192 total = qp->div->n_col - 2;
3193 if (total > g_pos) {
3194 int i;
3195 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3196 if (!exp)
3197 goto error;
3198 for (i = 0; i < total - g_pos; ++i)
3199 exp[i] = i + n;
3200 qp->upoly = expand(qp->upoly, exp, g_pos);
3201 free(exp);
3202 if (!qp->upoly)
3203 goto error;
3206 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3207 if (!qp->dim)
3208 goto error;
3210 return qp;
3211 error:
3212 isl_qpolynomial_free(qp);
3213 return NULL;
3216 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3217 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3219 unsigned pos;
3221 pos = isl_qpolynomial_dim(qp, type);
3223 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3226 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3227 __isl_take isl_pw_qpolynomial *pwqp,
3228 enum isl_dim_type type, unsigned n)
3230 unsigned pos;
3232 pos = isl_pw_qpolynomial_dim(pwqp, type);
3234 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3237 static int *reordering_move(isl_ctx *ctx,
3238 unsigned len, unsigned dst, unsigned src, unsigned n)
3240 int i;
3241 int *reordering;
3243 reordering = isl_alloc_array(ctx, int, len);
3244 if (!reordering)
3245 return NULL;
3247 if (dst <= src) {
3248 for (i = 0; i < dst; ++i)
3249 reordering[i] = i;
3250 for (i = 0; i < n; ++i)
3251 reordering[src + i] = dst + i;
3252 for (i = 0; i < src - dst; ++i)
3253 reordering[dst + i] = dst + n + i;
3254 for (i = 0; i < len - src - n; ++i)
3255 reordering[src + n + i] = src + n + i;
3256 } else {
3257 for (i = 0; i < src; ++i)
3258 reordering[i] = i;
3259 for (i = 0; i < n; ++i)
3260 reordering[src + i] = dst + i;
3261 for (i = 0; i < dst - src; ++i)
3262 reordering[src + n + i] = src + i;
3263 for (i = 0; i < len - dst - n; ++i)
3264 reordering[dst + n + i] = dst + n + i;
3267 return reordering;
3270 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3271 __isl_take isl_qpolynomial *qp,
3272 enum isl_dim_type dst_type, unsigned dst_pos,
3273 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3275 unsigned g_dst_pos;
3276 unsigned g_src_pos;
3277 int *reordering;
3279 if (!qp)
3280 return NULL;
3282 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3283 isl_die(qp->dim->ctx, isl_error_invalid,
3284 "cannot move output/set dimension",
3285 goto error);
3286 if (dst_type == isl_dim_in)
3287 dst_type = isl_dim_set;
3288 if (src_type == isl_dim_in)
3289 src_type = isl_dim_set;
3291 if (n == 0 &&
3292 !isl_space_is_named_or_nested(qp->dim, src_type) &&
3293 !isl_space_is_named_or_nested(qp->dim, dst_type))
3294 return qp;
3296 qp = isl_qpolynomial_cow(qp);
3297 if (!qp)
3298 return NULL;
3300 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3301 goto error);
3303 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3304 g_src_pos = pos(qp->dim, src_type) + src_pos;
3305 if (dst_type > src_type)
3306 g_dst_pos -= n;
3308 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3309 if (!qp->div)
3310 goto error;
3311 qp = sort_divs(qp);
3312 if (!qp)
3313 goto error;
3315 reordering = reordering_move(qp->dim->ctx,
3316 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3317 if (!reordering)
3318 goto error;
3320 qp->upoly = reorder(qp->upoly, reordering);
3321 free(reordering);
3322 if (!qp->upoly)
3323 goto error;
3325 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3326 if (!qp->dim)
3327 goto error;
3329 return qp;
3330 error:
3331 isl_qpolynomial_free(qp);
3332 return NULL;
3335 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3336 isl_int *f, isl_int denom)
3338 struct isl_upoly *up;
3340 dim = isl_space_domain(dim);
3341 if (!dim)
3342 return NULL;
3344 up = isl_upoly_from_affine(dim->ctx, f, denom,
3345 1 + isl_space_dim(dim, isl_dim_all));
3347 return isl_qpolynomial_alloc(dim, 0, up);
3350 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3352 isl_ctx *ctx;
3353 struct isl_upoly *up;
3354 isl_qpolynomial *qp;
3356 if (!aff)
3357 return NULL;
3359 ctx = isl_aff_get_ctx(aff);
3360 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3361 aff->v->size - 1);
3363 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3364 aff->ls->div->n_row, up);
3365 if (!qp)
3366 goto error;
3368 isl_mat_free(qp->div);
3369 qp->div = isl_mat_copy(aff->ls->div);
3370 qp->div = isl_mat_cow(qp->div);
3371 if (!qp->div)
3372 goto error;
3374 isl_aff_free(aff);
3375 qp = reduce_divs(qp);
3376 qp = remove_redundant_divs(qp);
3377 return qp;
3378 error:
3379 isl_aff_free(aff);
3380 return isl_qpolynomial_free(qp);
3383 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3384 __isl_take isl_pw_aff *pwaff)
3386 int i;
3387 isl_pw_qpolynomial *pwqp;
3389 if (!pwaff)
3390 return NULL;
3392 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3393 pwaff->n);
3395 for (i = 0; i < pwaff->n; ++i) {
3396 isl_set *dom;
3397 isl_qpolynomial *qp;
3399 dom = isl_set_copy(pwaff->p[i].set);
3400 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3401 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3404 isl_pw_aff_free(pwaff);
3405 return pwqp;
3408 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3409 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3411 isl_aff *aff;
3413 aff = isl_constraint_get_bound(c, type, pos);
3414 isl_constraint_free(c);
3415 return isl_qpolynomial_from_aff(aff);
3418 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3419 * in "qp" by subs[i].
3421 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3422 __isl_take isl_qpolynomial *qp,
3423 enum isl_dim_type type, unsigned first, unsigned n,
3424 __isl_keep isl_qpolynomial **subs)
3426 int i;
3427 struct isl_upoly **ups;
3429 if (n == 0)
3430 return qp;
3432 qp = isl_qpolynomial_cow(qp);
3433 if (!qp)
3434 return NULL;
3436 if (type == isl_dim_out)
3437 isl_die(qp->dim->ctx, isl_error_invalid,
3438 "cannot substitute output/set dimension",
3439 goto error);
3440 if (type == isl_dim_in)
3441 type = isl_dim_set;
3443 for (i = 0; i < n; ++i)
3444 if (!subs[i])
3445 goto error;
3447 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3448 goto error);
3450 for (i = 0; i < n; ++i)
3451 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3452 goto error);
3454 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3455 for (i = 0; i < n; ++i)
3456 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3458 first += pos(qp->dim, type);
3460 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3461 if (!ups)
3462 goto error;
3463 for (i = 0; i < n; ++i)
3464 ups[i] = subs[i]->upoly;
3466 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3468 free(ups);
3470 if (!qp->upoly)
3471 goto error;
3473 return qp;
3474 error:
3475 isl_qpolynomial_free(qp);
3476 return NULL;
3479 /* Extend "bset" with extra set dimensions for each integer division
3480 * in "qp" and then call "fn" with the extended bset and the polynomial
3481 * that results from replacing each of the integer divisions by the
3482 * corresponding extra set dimension.
3484 isl_stat isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3485 __isl_keep isl_basic_set *bset,
3486 isl_stat (*fn)(__isl_take isl_basic_set *bset,
3487 __isl_take isl_qpolynomial *poly, void *user), void *user)
3489 isl_space *dim;
3490 isl_mat *div;
3491 isl_qpolynomial *poly;
3493 if (!qp || !bset)
3494 return isl_stat_error;
3495 if (qp->div->n_row == 0)
3496 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3497 user);
3499 div = isl_mat_copy(qp->div);
3500 dim = isl_space_copy(qp->dim);
3501 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3502 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3503 bset = isl_basic_set_copy(bset);
3504 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3505 bset = add_div_constraints(bset, div);
3507 return fn(bset, poly, user);
3510 /* Return total degree in variables first (inclusive) up to last (exclusive).
3512 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3514 int deg = -1;
3515 int i;
3516 struct isl_upoly_rec *rec;
3518 if (!up)
3519 return -2;
3520 if (isl_upoly_is_zero(up))
3521 return -1;
3522 if (isl_upoly_is_cst(up) || up->var < first)
3523 return 0;
3525 rec = isl_upoly_as_rec(up);
3526 if (!rec)
3527 return -2;
3529 for (i = 0; i < rec->n; ++i) {
3530 int d;
3532 if (isl_upoly_is_zero(rec->p[i]))
3533 continue;
3534 d = isl_upoly_degree(rec->p[i], first, last);
3535 if (up->var < last)
3536 d += i;
3537 if (d > deg)
3538 deg = d;
3541 return deg;
3544 /* Return total degree in set variables.
3546 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3548 unsigned ovar;
3549 unsigned nvar;
3551 if (!poly)
3552 return -2;
3554 ovar = isl_space_offset(poly->dim, isl_dim_set);
3555 nvar = isl_space_dim(poly->dim, isl_dim_set);
3556 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3559 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3560 unsigned pos, int deg)
3562 int i;
3563 struct isl_upoly_rec *rec;
3565 if (!up)
3566 return NULL;
3568 if (isl_upoly_is_cst(up) || up->var < pos) {
3569 if (deg == 0)
3570 return isl_upoly_copy(up);
3571 else
3572 return isl_upoly_zero(up->ctx);
3575 rec = isl_upoly_as_rec(up);
3576 if (!rec)
3577 return NULL;
3579 if (up->var == pos) {
3580 if (deg < rec->n)
3581 return isl_upoly_copy(rec->p[deg]);
3582 else
3583 return isl_upoly_zero(up->ctx);
3586 up = isl_upoly_copy(up);
3587 up = isl_upoly_cow(up);
3588 rec = isl_upoly_as_rec(up);
3589 if (!rec)
3590 goto error;
3592 for (i = 0; i < rec->n; ++i) {
3593 struct isl_upoly *t;
3594 t = isl_upoly_coeff(rec->p[i], pos, deg);
3595 if (!t)
3596 goto error;
3597 isl_upoly_free(rec->p[i]);
3598 rec->p[i] = t;
3601 return up;
3602 error:
3603 isl_upoly_free(up);
3604 return NULL;
3607 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3609 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3610 __isl_keep isl_qpolynomial *qp,
3611 enum isl_dim_type type, unsigned t_pos, int deg)
3613 unsigned g_pos;
3614 struct isl_upoly *up;
3615 isl_qpolynomial *c;
3617 if (!qp)
3618 return NULL;
3620 if (type == isl_dim_out)
3621 isl_die(qp->div->ctx, isl_error_invalid,
3622 "output/set dimension does not have a coefficient",
3623 return NULL);
3624 if (type == isl_dim_in)
3625 type = isl_dim_set;
3627 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3628 return NULL);
3630 g_pos = pos(qp->dim, type) + t_pos;
3631 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3633 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3634 if (!c)
3635 return NULL;
3636 isl_mat_free(c->div);
3637 c->div = isl_mat_copy(qp->div);
3638 if (!c->div)
3639 goto error;
3640 return c;
3641 error:
3642 isl_qpolynomial_free(c);
3643 return NULL;
3646 /* Homogenize the polynomial in the variables first (inclusive) up to
3647 * last (exclusive) by inserting powers of variable first.
3648 * Variable first is assumed not to appear in the input.
3650 __isl_give struct isl_upoly *isl_upoly_homogenize(
3651 __isl_take struct isl_upoly *up, int deg, int target,
3652 int first, int last)
3654 int i;
3655 struct isl_upoly_rec *rec;
3657 if (!up)
3658 return NULL;
3659 if (isl_upoly_is_zero(up))
3660 return up;
3661 if (deg == target)
3662 return up;
3663 if (isl_upoly_is_cst(up) || up->var < first) {
3664 struct isl_upoly *hom;
3666 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3667 if (!hom)
3668 goto error;
3669 rec = isl_upoly_as_rec(hom);
3670 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3672 return hom;
3675 up = isl_upoly_cow(up);
3676 rec = isl_upoly_as_rec(up);
3677 if (!rec)
3678 goto error;
3680 for (i = 0; i < rec->n; ++i) {
3681 if (isl_upoly_is_zero(rec->p[i]))
3682 continue;
3683 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3684 up->var < last ? deg + i : i, target,
3685 first, last);
3686 if (!rec->p[i])
3687 goto error;
3690 return up;
3691 error:
3692 isl_upoly_free(up);
3693 return NULL;
3696 /* Homogenize the polynomial in the set variables by introducing
3697 * powers of an extra set variable at position 0.
3699 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3700 __isl_take isl_qpolynomial *poly)
3702 unsigned ovar;
3703 unsigned nvar;
3704 int deg = isl_qpolynomial_degree(poly);
3706 if (deg < -1)
3707 goto error;
3709 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3710 poly = isl_qpolynomial_cow(poly);
3711 if (!poly)
3712 goto error;
3714 ovar = isl_space_offset(poly->dim, isl_dim_set);
3715 nvar = isl_space_dim(poly->dim, isl_dim_set);
3716 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3717 ovar, ovar + nvar);
3718 if (!poly->upoly)
3719 goto error;
3721 return poly;
3722 error:
3723 isl_qpolynomial_free(poly);
3724 return NULL;
3727 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3728 __isl_take isl_mat *div)
3730 isl_term *term;
3731 int n;
3733 if (!dim || !div)
3734 goto error;
3736 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3738 term = isl_calloc(dim->ctx, struct isl_term,
3739 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3740 if (!term)
3741 goto error;
3743 term->ref = 1;
3744 term->dim = dim;
3745 term->div = div;
3746 isl_int_init(term->n);
3747 isl_int_init(term->d);
3749 return term;
3750 error:
3751 isl_space_free(dim);
3752 isl_mat_free(div);
3753 return NULL;
3756 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3758 if (!term)
3759 return NULL;
3761 term->ref++;
3762 return term;
3765 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3767 int i;
3768 isl_term *dup;
3769 unsigned total;
3771 if (!term)
3772 return NULL;
3774 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3776 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3777 if (!dup)
3778 return NULL;
3780 isl_int_set(dup->n, term->n);
3781 isl_int_set(dup->d, term->d);
3783 for (i = 0; i < total; ++i)
3784 dup->pow[i] = term->pow[i];
3786 return dup;
3789 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3791 if (!term)
3792 return NULL;
3794 if (term->ref == 1)
3795 return term;
3796 term->ref--;
3797 return isl_term_dup(term);
3800 void isl_term_free(__isl_take isl_term *term)
3802 if (!term)
3803 return;
3805 if (--term->ref > 0)
3806 return;
3808 isl_space_free(term->dim);
3809 isl_mat_free(term->div);
3810 isl_int_clear(term->n);
3811 isl_int_clear(term->d);
3812 free(term);
3815 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3817 if (!term)
3818 return 0;
3820 switch (type) {
3821 case isl_dim_param:
3822 case isl_dim_in:
3823 case isl_dim_out: return isl_space_dim(term->dim, type);
3824 case isl_dim_div: return term->div->n_row;
3825 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3826 term->div->n_row;
3827 default: return 0;
3831 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3833 return term ? term->dim->ctx : NULL;
3836 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3838 if (!term)
3839 return;
3840 isl_int_set(*n, term->n);
3843 /* Return the coefficient of the term "term".
3845 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3847 if (!term)
3848 return NULL;
3850 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3851 term->n, term->d);
3854 int isl_term_get_exp(__isl_keep isl_term *term,
3855 enum isl_dim_type type, unsigned pos)
3857 if (!term)
3858 return -1;
3860 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3862 if (type >= isl_dim_set)
3863 pos += isl_space_dim(term->dim, isl_dim_param);
3864 if (type >= isl_dim_div)
3865 pos += isl_space_dim(term->dim, isl_dim_set);
3867 return term->pow[pos];
3870 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3872 isl_local_space *ls;
3873 isl_aff *aff;
3875 if (!term)
3876 return NULL;
3878 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3879 return NULL);
3881 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3882 isl_mat_copy(term->div));
3883 aff = isl_aff_alloc(ls);
3884 if (!aff)
3885 return NULL;
3887 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3889 aff = isl_aff_normalize(aff);
3891 return aff;
3894 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3895 isl_stat (*fn)(__isl_take isl_term *term, void *user),
3896 __isl_take isl_term *term, void *user)
3898 int i;
3899 struct isl_upoly_rec *rec;
3901 if (!up || !term)
3902 goto error;
3904 if (isl_upoly_is_zero(up))
3905 return term;
3907 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3908 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3909 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3911 if (isl_upoly_is_cst(up)) {
3912 struct isl_upoly_cst *cst;
3913 cst = isl_upoly_as_cst(up);
3914 if (!cst)
3915 goto error;
3916 term = isl_term_cow(term);
3917 if (!term)
3918 goto error;
3919 isl_int_set(term->n, cst->n);
3920 isl_int_set(term->d, cst->d);
3921 if (fn(isl_term_copy(term), user) < 0)
3922 goto error;
3923 return term;
3926 rec = isl_upoly_as_rec(up);
3927 if (!rec)
3928 goto error;
3930 for (i = 0; i < rec->n; ++i) {
3931 term = isl_term_cow(term);
3932 if (!term)
3933 goto error;
3934 term->pow[up->var] = i;
3935 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3936 if (!term)
3937 goto error;
3939 term->pow[up->var] = 0;
3941 return term;
3942 error:
3943 isl_term_free(term);
3944 return NULL;
3947 isl_stat isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3948 isl_stat (*fn)(__isl_take isl_term *term, void *user), void *user)
3950 isl_term *term;
3952 if (!qp)
3953 return isl_stat_error;
3955 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3956 if (!term)
3957 return isl_stat_error;
3959 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3961 isl_term_free(term);
3963 return term ? isl_stat_ok : isl_stat_error;
3966 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3968 struct isl_upoly *up;
3969 isl_qpolynomial *qp;
3970 int i, n;
3972 if (!term)
3973 return NULL;
3975 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3977 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3978 for (i = 0; i < n; ++i) {
3979 if (!term->pow[i])
3980 continue;
3981 up = isl_upoly_mul(up,
3982 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3985 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3986 if (!qp)
3987 goto error;
3988 isl_mat_free(qp->div);
3989 qp->div = isl_mat_copy(term->div);
3990 if (!qp->div)
3991 goto error;
3993 isl_term_free(term);
3994 return qp;
3995 error:
3996 isl_qpolynomial_free(qp);
3997 isl_term_free(term);
3998 return NULL;
4001 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
4002 __isl_take isl_space *dim)
4004 int i;
4005 int extra;
4006 unsigned total;
4008 if (!qp || !dim)
4009 goto error;
4011 if (isl_space_is_equal(qp->dim, dim)) {
4012 isl_space_free(dim);
4013 return qp;
4016 qp = isl_qpolynomial_cow(qp);
4017 if (!qp)
4018 goto error;
4020 extra = isl_space_dim(dim, isl_dim_set) -
4021 isl_space_dim(qp->dim, isl_dim_set);
4022 total = isl_space_dim(qp->dim, isl_dim_all);
4023 if (qp->div->n_row) {
4024 int *exp;
4026 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
4027 if (!exp)
4028 goto error;
4029 for (i = 0; i < qp->div->n_row; ++i)
4030 exp[i] = extra + i;
4031 qp->upoly = expand(qp->upoly, exp, total);
4032 free(exp);
4033 if (!qp->upoly)
4034 goto error;
4036 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
4037 if (!qp->div)
4038 goto error;
4039 for (i = 0; i < qp->div->n_row; ++i)
4040 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
4042 isl_space_free(qp->dim);
4043 qp->dim = dim;
4045 return qp;
4046 error:
4047 isl_space_free(dim);
4048 isl_qpolynomial_free(qp);
4049 return NULL;
4052 /* For each parameter or variable that does not appear in qp,
4053 * first eliminate the variable from all constraints and then set it to zero.
4055 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
4056 __isl_keep isl_qpolynomial *qp)
4058 int *active = NULL;
4059 int i;
4060 int d;
4061 unsigned nparam;
4062 unsigned nvar;
4064 if (!set || !qp)
4065 goto error;
4067 d = isl_space_dim(set->dim, isl_dim_all);
4068 active = isl_calloc_array(set->ctx, int, d);
4069 if (set_active(qp, active) < 0)
4070 goto error;
4072 for (i = 0; i < d; ++i)
4073 if (!active[i])
4074 break;
4076 if (i == d) {
4077 free(active);
4078 return set;
4081 nparam = isl_space_dim(set->dim, isl_dim_param);
4082 nvar = isl_space_dim(set->dim, isl_dim_set);
4083 for (i = 0; i < nparam; ++i) {
4084 if (active[i])
4085 continue;
4086 set = isl_set_eliminate(set, isl_dim_param, i, 1);
4087 set = isl_set_fix_si(set, isl_dim_param, i, 0);
4089 for (i = 0; i < nvar; ++i) {
4090 if (active[nparam + i])
4091 continue;
4092 set = isl_set_eliminate(set, isl_dim_set, i, 1);
4093 set = isl_set_fix_si(set, isl_dim_set, i, 0);
4096 free(active);
4098 return set;
4099 error:
4100 free(active);
4101 isl_set_free(set);
4102 return NULL;
4105 struct isl_opt_data {
4106 isl_qpolynomial *qp;
4107 int first;
4108 isl_val *opt;
4109 int max;
4112 static isl_stat opt_fn(__isl_take isl_point *pnt, void *user)
4114 struct isl_opt_data *data = (struct isl_opt_data *)user;
4115 isl_val *val;
4117 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
4118 if (data->first) {
4119 data->first = 0;
4120 data->opt = val;
4121 } else if (data->max) {
4122 data->opt = isl_val_max(data->opt, val);
4123 } else {
4124 data->opt = isl_val_min(data->opt, val);
4127 return isl_stat_ok;
4130 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
4131 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
4133 struct isl_opt_data data = { NULL, 1, NULL, max };
4135 if (!set || !qp)
4136 goto error;
4138 if (isl_upoly_is_cst(qp->upoly)) {
4139 isl_set_free(set);
4140 data.opt = isl_qpolynomial_get_constant_val(qp);
4141 isl_qpolynomial_free(qp);
4142 return data.opt;
4145 set = fix_inactive(set, qp);
4147 data.qp = qp;
4148 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
4149 goto error;
4151 if (data.first)
4152 data.opt = isl_val_zero(isl_set_get_ctx(set));
4154 isl_set_free(set);
4155 isl_qpolynomial_free(qp);
4156 return data.opt;
4157 error:
4158 isl_set_free(set);
4159 isl_qpolynomial_free(qp);
4160 isl_val_free(data.opt);
4161 return NULL;
4164 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
4165 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
4167 int i;
4168 int n_sub;
4169 isl_ctx *ctx;
4170 struct isl_upoly **subs;
4171 isl_mat *mat, *diag;
4173 qp = isl_qpolynomial_cow(qp);
4174 if (!qp || !morph)
4175 goto error;
4177 ctx = qp->dim->ctx;
4178 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
4180 n_sub = morph->inv->n_row - 1;
4181 if (morph->inv->n_row != morph->inv->n_col)
4182 n_sub += qp->div->n_row;
4183 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
4184 if (n_sub && !subs)
4185 goto error;
4187 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4188 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4189 morph->inv->row[0][0], morph->inv->n_col);
4190 if (morph->inv->n_row != morph->inv->n_col)
4191 for (i = 0; i < qp->div->n_row; ++i)
4192 subs[morph->inv->n_row - 1 + i] =
4193 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4195 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4197 for (i = 0; i < n_sub; ++i)
4198 isl_upoly_free(subs[i]);
4199 free(subs);
4201 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4202 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4203 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4204 mat = isl_mat_diagonal(mat, diag);
4205 qp->div = isl_mat_product(qp->div, mat);
4206 isl_space_free(qp->dim);
4207 qp->dim = isl_space_copy(morph->ran->dim);
4209 if (!qp->upoly || !qp->div || !qp->dim)
4210 goto error;
4212 isl_morph_free(morph);
4214 return qp;
4215 error:
4216 isl_qpolynomial_free(qp);
4217 isl_morph_free(morph);
4218 return NULL;
4221 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4222 __isl_take isl_union_pw_qpolynomial *upwqp1,
4223 __isl_take isl_union_pw_qpolynomial *upwqp2)
4225 return isl_union_pw_qpolynomial_match_bin_op(upwqp1, upwqp2,
4226 &isl_pw_qpolynomial_mul);
4229 /* Reorder the columns of the given div definitions according to the
4230 * given reordering.
4232 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4233 __isl_take isl_reordering *r)
4235 int i, j;
4236 isl_mat *mat;
4237 int extra;
4239 if (!div || !r)
4240 goto error;
4242 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4243 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4244 if (!mat)
4245 goto error;
4247 for (i = 0; i < div->n_row; ++i) {
4248 isl_seq_cpy(mat->row[i], div->row[i], 2);
4249 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4250 for (j = 0; j < r->len; ++j)
4251 isl_int_set(mat->row[i][2 + r->pos[j]],
4252 div->row[i][2 + j]);
4255 isl_reordering_free(r);
4256 isl_mat_free(div);
4257 return mat;
4258 error:
4259 isl_reordering_free(r);
4260 isl_mat_free(div);
4261 return NULL;
4264 /* Reorder the dimension of "qp" according to the given reordering.
4266 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4267 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4269 qp = isl_qpolynomial_cow(qp);
4270 if (!qp)
4271 goto error;
4273 r = isl_reordering_extend(r, qp->div->n_row);
4274 if (!r)
4275 goto error;
4277 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4278 if (!qp->div)
4279 goto error;
4281 qp->upoly = reorder(qp->upoly, r->pos);
4282 if (!qp->upoly)
4283 goto error;
4285 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4287 isl_reordering_free(r);
4288 return qp;
4289 error:
4290 isl_qpolynomial_free(qp);
4291 isl_reordering_free(r);
4292 return NULL;
4295 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4296 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4298 isl_bool equal_params;
4300 if (!qp || !model)
4301 goto error;
4303 equal_params = isl_space_has_equal_params(qp->dim, model);
4304 if (equal_params < 0)
4305 goto error;
4306 if (!equal_params) {
4307 isl_reordering *exp;
4309 model = isl_space_drop_dims(model, isl_dim_in,
4310 0, isl_space_dim(model, isl_dim_in));
4311 model = isl_space_drop_dims(model, isl_dim_out,
4312 0, isl_space_dim(model, isl_dim_out));
4313 exp = isl_parameter_alignment_reordering(qp->dim, model);
4314 exp = isl_reordering_extend_space(exp,
4315 isl_qpolynomial_get_domain_space(qp));
4316 qp = isl_qpolynomial_realign_domain(qp, exp);
4319 isl_space_free(model);
4320 return qp;
4321 error:
4322 isl_space_free(model);
4323 isl_qpolynomial_free(qp);
4324 return NULL;
4327 struct isl_split_periods_data {
4328 int max_periods;
4329 isl_pw_qpolynomial *res;
4332 /* Create a slice where the integer division "div" has the fixed value "v".
4333 * In particular, if "div" refers to floor(f/m), then create a slice
4335 * m v <= f <= m v + (m - 1)
4337 * or
4339 * f - m v >= 0
4340 * -f + m v + (m - 1) >= 0
4342 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4343 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4345 int total;
4346 isl_basic_set *bset = NULL;
4347 int k;
4349 if (!dim || !qp)
4350 goto error;
4352 total = isl_space_dim(dim, isl_dim_all);
4353 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4355 k = isl_basic_set_alloc_inequality(bset);
4356 if (k < 0)
4357 goto error;
4358 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4359 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4361 k = isl_basic_set_alloc_inequality(bset);
4362 if (k < 0)
4363 goto error;
4364 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4365 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4366 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4367 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4369 isl_space_free(dim);
4370 return isl_set_from_basic_set(bset);
4371 error:
4372 isl_basic_set_free(bset);
4373 isl_space_free(dim);
4374 return NULL;
4377 static isl_stat split_periods(__isl_take isl_set *set,
4378 __isl_take isl_qpolynomial *qp, void *user);
4380 /* Create a slice of the domain "set" such that integer division "div"
4381 * has the fixed value "v" and add the results to data->res,
4382 * replacing the integer division by "v" in "qp".
4384 static isl_stat set_div(__isl_take isl_set *set,
4385 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4386 struct isl_split_periods_data *data)
4388 int i;
4389 int total;
4390 isl_set *slice;
4391 struct isl_upoly *cst;
4393 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4394 set = isl_set_intersect(set, slice);
4396 if (!qp)
4397 goto error;
4399 total = isl_space_dim(qp->dim, isl_dim_all);
4401 for (i = div + 1; i < qp->div->n_row; ++i) {
4402 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4403 continue;
4404 isl_int_addmul(qp->div->row[i][1],
4405 qp->div->row[i][2 + total + div], v);
4406 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4409 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4410 qp = substitute_div(qp, div, cst);
4412 return split_periods(set, qp, data);
4413 error:
4414 isl_set_free(set);
4415 isl_qpolynomial_free(qp);
4416 return -1;
4419 /* Split the domain "set" such that integer division "div"
4420 * has a fixed value (ranging from "min" to "max") on each slice
4421 * and add the results to data->res.
4423 static isl_stat split_div(__isl_take isl_set *set,
4424 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4425 struct isl_split_periods_data *data)
4427 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4428 isl_set *set_i = isl_set_copy(set);
4429 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4431 if (set_div(set_i, qp_i, div, min, data) < 0)
4432 goto error;
4434 isl_set_free(set);
4435 isl_qpolynomial_free(qp);
4436 return isl_stat_ok;
4437 error:
4438 isl_set_free(set);
4439 isl_qpolynomial_free(qp);
4440 return isl_stat_error;
4443 /* If "qp" refers to any integer division
4444 * that can only attain "max_periods" distinct values on "set"
4445 * then split the domain along those distinct values.
4446 * Add the results (or the original if no splitting occurs)
4447 * to data->res.
4449 static isl_stat split_periods(__isl_take isl_set *set,
4450 __isl_take isl_qpolynomial *qp, void *user)
4452 int i;
4453 isl_pw_qpolynomial *pwqp;
4454 struct isl_split_periods_data *data;
4455 isl_int min, max;
4456 int total;
4457 isl_stat r = isl_stat_ok;
4459 data = (struct isl_split_periods_data *)user;
4461 if (!set || !qp)
4462 goto error;
4464 if (qp->div->n_row == 0) {
4465 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4466 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4467 return isl_stat_ok;
4470 isl_int_init(min);
4471 isl_int_init(max);
4472 total = isl_space_dim(qp->dim, isl_dim_all);
4473 for (i = 0; i < qp->div->n_row; ++i) {
4474 enum isl_lp_result lp_res;
4476 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4477 qp->div->n_row) != -1)
4478 continue;
4480 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4481 set->ctx->one, &min, NULL, NULL);
4482 if (lp_res == isl_lp_error)
4483 goto error2;
4484 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4485 continue;
4486 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4488 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4489 set->ctx->one, &max, NULL, NULL);
4490 if (lp_res == isl_lp_error)
4491 goto error2;
4492 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4493 continue;
4494 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4496 isl_int_sub(max, max, min);
4497 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4498 isl_int_add(max, max, min);
4499 break;
4503 if (i < qp->div->n_row) {
4504 r = split_div(set, qp, i, min, max, data);
4505 } else {
4506 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4507 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4510 isl_int_clear(max);
4511 isl_int_clear(min);
4513 return r;
4514 error2:
4515 isl_int_clear(max);
4516 isl_int_clear(min);
4517 error:
4518 isl_set_free(set);
4519 isl_qpolynomial_free(qp);
4520 return isl_stat_error;
4523 /* If any quasi-polynomial in pwqp refers to any integer division
4524 * that can only attain "max_periods" distinct values on its domain
4525 * then split the domain along those distinct values.
4527 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4528 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4530 struct isl_split_periods_data data;
4532 data.max_periods = max_periods;
4533 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4535 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4536 goto error;
4538 isl_pw_qpolynomial_free(pwqp);
4540 return data.res;
4541 error:
4542 isl_pw_qpolynomial_free(data.res);
4543 isl_pw_qpolynomial_free(pwqp);
4544 return NULL;
4547 /* Construct a piecewise quasipolynomial that is constant on the given
4548 * domain. In particular, it is
4549 * 0 if cst == 0
4550 * 1 if cst == 1
4551 * infinity if cst == -1
4553 * If cst == -1, then explicitly check whether the domain is empty and,
4554 * if so, return 0 instead.
4556 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4557 __isl_take isl_basic_set *bset, int cst)
4559 isl_space *dim;
4560 isl_qpolynomial *qp;
4562 if (cst < 0 && isl_basic_set_is_empty(bset) == isl_bool_true)
4563 cst = 0;
4564 if (!bset)
4565 return NULL;
4567 bset = isl_basic_set_params(bset);
4568 dim = isl_basic_set_get_space(bset);
4569 if (cst < 0)
4570 qp = isl_qpolynomial_infty_on_domain(dim);
4571 else if (cst == 0)
4572 qp = isl_qpolynomial_zero_on_domain(dim);
4573 else
4574 qp = isl_qpolynomial_one_on_domain(dim);
4575 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4578 /* Factor bset, call fn on each of the factors and return the product.
4580 * If no factors can be found, simply call fn on the input.
4581 * Otherwise, construct the factors based on the factorizer,
4582 * call fn on each factor and compute the product.
4584 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4585 __isl_take isl_basic_set *bset,
4586 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4588 int i, n;
4589 isl_space *space;
4590 isl_set *set;
4591 isl_factorizer *f;
4592 isl_qpolynomial *qp;
4593 isl_pw_qpolynomial *pwqp;
4594 unsigned nparam;
4595 unsigned nvar;
4597 f = isl_basic_set_factorizer(bset);
4598 if (!f)
4599 goto error;
4600 if (f->n_group == 0) {
4601 isl_factorizer_free(f);
4602 return fn(bset);
4605 nparam = isl_basic_set_dim(bset, isl_dim_param);
4606 nvar = isl_basic_set_dim(bset, isl_dim_set);
4608 space = isl_basic_set_get_space(bset);
4609 space = isl_space_params(space);
4610 set = isl_set_universe(isl_space_copy(space));
4611 qp = isl_qpolynomial_one_on_domain(space);
4612 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4614 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4616 for (i = 0, n = 0; i < f->n_group; ++i) {
4617 isl_basic_set *bset_i;
4618 isl_pw_qpolynomial *pwqp_i;
4620 bset_i = isl_basic_set_copy(bset);
4621 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4622 nparam + n + f->len[i], nvar - n - f->len[i]);
4623 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4624 nparam, n);
4625 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4626 n + f->len[i], nvar - n - f->len[i]);
4627 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4629 pwqp_i = fn(bset_i);
4630 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4632 n += f->len[i];
4635 isl_basic_set_free(bset);
4636 isl_factorizer_free(f);
4638 return pwqp;
4639 error:
4640 isl_basic_set_free(bset);
4641 return NULL;
4644 /* Factor bset, call fn on each of the factors and return the product.
4645 * The function is assumed to evaluate to zero on empty domains,
4646 * to one on zero-dimensional domains and to infinity on unbounded domains
4647 * and will not be called explicitly on zero-dimensional or unbounded domains.
4649 * We first check for some special cases and remove all equalities.
4650 * Then we hand over control to compressed_multiplicative_call.
4652 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4653 __isl_take isl_basic_set *bset,
4654 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4656 isl_bool bounded;
4657 isl_morph *morph;
4658 isl_pw_qpolynomial *pwqp;
4660 if (!bset)
4661 return NULL;
4663 if (isl_basic_set_plain_is_empty(bset))
4664 return constant_on_domain(bset, 0);
4666 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4667 return constant_on_domain(bset, 1);
4669 bounded = isl_basic_set_is_bounded(bset);
4670 if (bounded < 0)
4671 goto error;
4672 if (!bounded)
4673 return constant_on_domain(bset, -1);
4675 if (bset->n_eq == 0)
4676 return compressed_multiplicative_call(bset, fn);
4678 morph = isl_basic_set_full_compression(bset);
4679 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4681 pwqp = compressed_multiplicative_call(bset, fn);
4683 morph = isl_morph_dom_params(morph);
4684 morph = isl_morph_ran_params(morph);
4685 morph = isl_morph_inverse(morph);
4687 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4689 return pwqp;
4690 error:
4691 isl_basic_set_free(bset);
4692 return NULL;
4695 /* Drop all floors in "qp", turning each integer division [a/m] into
4696 * a rational division a/m. If "down" is set, then the integer division
4697 * is replaced by (a-(m-1))/m instead.
4699 static __isl_give isl_qpolynomial *qp_drop_floors(
4700 __isl_take isl_qpolynomial *qp, int down)
4702 int i;
4703 struct isl_upoly *s;
4705 if (!qp)
4706 return NULL;
4707 if (qp->div->n_row == 0)
4708 return qp;
4710 qp = isl_qpolynomial_cow(qp);
4711 if (!qp)
4712 return NULL;
4714 for (i = qp->div->n_row - 1; i >= 0; --i) {
4715 if (down) {
4716 isl_int_sub(qp->div->row[i][1],
4717 qp->div->row[i][1], qp->div->row[i][0]);
4718 isl_int_add_ui(qp->div->row[i][1],
4719 qp->div->row[i][1], 1);
4721 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4722 qp->div->row[i][0], qp->div->n_col - 1);
4723 qp = substitute_div(qp, i, s);
4724 if (!qp)
4725 return NULL;
4728 return qp;
4731 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4732 * a rational division a/m.
4734 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4735 __isl_take isl_pw_qpolynomial *pwqp)
4737 int i;
4739 if (!pwqp)
4740 return NULL;
4742 if (isl_pw_qpolynomial_is_zero(pwqp))
4743 return pwqp;
4745 pwqp = isl_pw_qpolynomial_cow(pwqp);
4746 if (!pwqp)
4747 return NULL;
4749 for (i = 0; i < pwqp->n; ++i) {
4750 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4751 if (!pwqp->p[i].qp)
4752 goto error;
4755 return pwqp;
4756 error:
4757 isl_pw_qpolynomial_free(pwqp);
4758 return NULL;
4761 /* Adjust all the integer divisions in "qp" such that they are at least
4762 * one over the given orthant (identified by "signs"). This ensures
4763 * that they will still be non-negative even after subtracting (m-1)/m.
4765 * In particular, f is replaced by f' + v, changing f = [a/m]
4766 * to f' = [(a - m v)/m].
4767 * If the constant term k in a is smaller than m,
4768 * the constant term of v is set to floor(k/m) - 1.
4769 * For any other term, if the coefficient c and the variable x have
4770 * the same sign, then no changes are needed.
4771 * Otherwise, if the variable is positive (and c is negative),
4772 * then the coefficient of x in v is set to floor(c/m).
4773 * If the variable is negative (and c is positive),
4774 * then the coefficient of x in v is set to ceil(c/m).
4776 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4777 int *signs)
4779 int i, j;
4780 int total;
4781 isl_vec *v = NULL;
4782 struct isl_upoly *s;
4784 qp = isl_qpolynomial_cow(qp);
4785 if (!qp)
4786 return NULL;
4787 qp->div = isl_mat_cow(qp->div);
4788 if (!qp->div)
4789 goto error;
4791 total = isl_space_dim(qp->dim, isl_dim_all);
4792 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4794 for (i = 0; i < qp->div->n_row; ++i) {
4795 isl_int *row = qp->div->row[i];
4796 v = isl_vec_clr(v);
4797 if (!v)
4798 goto error;
4799 if (isl_int_lt(row[1], row[0])) {
4800 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4801 isl_int_sub_ui(v->el[0], v->el[0], 1);
4802 isl_int_submul(row[1], row[0], v->el[0]);
4804 for (j = 0; j < total; ++j) {
4805 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4806 continue;
4807 if (signs[j] < 0)
4808 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4809 else
4810 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4811 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4813 for (j = 0; j < i; ++j) {
4814 if (isl_int_sgn(row[2 + total + j]) >= 0)
4815 continue;
4816 isl_int_fdiv_q(v->el[1 + total + j],
4817 row[2 + total + j], row[0]);
4818 isl_int_submul(row[2 + total + j],
4819 row[0], v->el[1 + total + j]);
4821 for (j = i + 1; j < qp->div->n_row; ++j) {
4822 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4823 continue;
4824 isl_seq_combine(qp->div->row[j] + 1,
4825 qp->div->ctx->one, qp->div->row[j] + 1,
4826 qp->div->row[j][2 + total + i], v->el, v->size);
4828 isl_int_set_si(v->el[1 + total + i], 1);
4829 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4830 qp->div->ctx->one, v->size);
4831 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4832 isl_upoly_free(s);
4833 if (!qp->upoly)
4834 goto error;
4837 isl_vec_free(v);
4838 return qp;
4839 error:
4840 isl_vec_free(v);
4841 isl_qpolynomial_free(qp);
4842 return NULL;
4845 struct isl_to_poly_data {
4846 int sign;
4847 isl_pw_qpolynomial *res;
4848 isl_qpolynomial *qp;
4851 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4852 * We first make all integer divisions positive and then split the
4853 * quasipolynomials into terms with sign data->sign (the direction
4854 * of the requested approximation) and terms with the opposite sign.
4855 * In the first set of terms, each integer division [a/m] is
4856 * overapproximated by a/m, while in the second it is underapproximated
4857 * by (a-(m-1))/m.
4859 static isl_stat to_polynomial_on_orthant(__isl_take isl_set *orthant,
4860 int *signs, void *user)
4862 struct isl_to_poly_data *data = user;
4863 isl_pw_qpolynomial *t;
4864 isl_qpolynomial *qp, *up, *down;
4866 qp = isl_qpolynomial_copy(data->qp);
4867 qp = make_divs_pos(qp, signs);
4869 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4870 up = qp_drop_floors(up, 0);
4871 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4872 down = qp_drop_floors(down, 1);
4874 isl_qpolynomial_free(qp);
4875 qp = isl_qpolynomial_add(up, down);
4877 t = isl_pw_qpolynomial_alloc(orthant, qp);
4878 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4880 return isl_stat_ok;
4883 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4884 * the polynomial will be an overapproximation. If "sign" is negative,
4885 * it will be an underapproximation. If "sign" is zero, the approximation
4886 * will lie somewhere in between.
4888 * In particular, is sign == 0, we simply drop the floors, turning
4889 * the integer divisions into rational divisions.
4890 * Otherwise, we split the domains into orthants, make all integer divisions
4891 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4892 * depending on the requested sign and the sign of the term in which
4893 * the integer division appears.
4895 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4896 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4898 int i;
4899 struct isl_to_poly_data data;
4901 if (sign == 0)
4902 return pwqp_drop_floors(pwqp);
4904 if (!pwqp)
4905 return NULL;
4907 data.sign = sign;
4908 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4910 for (i = 0; i < pwqp->n; ++i) {
4911 if (pwqp->p[i].qp->div->n_row == 0) {
4912 isl_pw_qpolynomial *t;
4913 t = isl_pw_qpolynomial_alloc(
4914 isl_set_copy(pwqp->p[i].set),
4915 isl_qpolynomial_copy(pwqp->p[i].qp));
4916 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4917 continue;
4919 data.qp = pwqp->p[i].qp;
4920 if (isl_set_foreach_orthant(pwqp->p[i].set,
4921 &to_polynomial_on_orthant, &data) < 0)
4922 goto error;
4925 isl_pw_qpolynomial_free(pwqp);
4927 return data.res;
4928 error:
4929 isl_pw_qpolynomial_free(pwqp);
4930 isl_pw_qpolynomial_free(data.res);
4931 return NULL;
4934 static __isl_give isl_pw_qpolynomial *poly_entry(
4935 __isl_take isl_pw_qpolynomial *pwqp, void *user)
4937 int *sign = user;
4939 return isl_pw_qpolynomial_to_polynomial(pwqp, *sign);
4942 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4943 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4945 return isl_union_pw_qpolynomial_transform_inplace(upwqp,
4946 &poly_entry, &sign);
4949 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4950 __isl_take isl_qpolynomial *qp)
4952 int i, k;
4953 isl_space *dim;
4954 isl_vec *aff = NULL;
4955 isl_basic_map *bmap = NULL;
4956 unsigned pos;
4957 unsigned n_div;
4959 if (!qp)
4960 return NULL;
4961 if (!isl_upoly_is_affine(qp->upoly))
4962 isl_die(qp->dim->ctx, isl_error_invalid,
4963 "input quasi-polynomial not affine", goto error);
4964 aff = isl_qpolynomial_extract_affine(qp);
4965 if (!aff)
4966 goto error;
4967 dim = isl_qpolynomial_get_space(qp);
4968 pos = 1 + isl_space_offset(dim, isl_dim_out);
4969 n_div = qp->div->n_row;
4970 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4972 for (i = 0; i < n_div; ++i) {
4973 k = isl_basic_map_alloc_div(bmap);
4974 if (k < 0)
4975 goto error;
4976 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4977 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4978 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4979 goto error;
4981 k = isl_basic_map_alloc_equality(bmap);
4982 if (k < 0)
4983 goto error;
4984 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4985 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4986 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4988 isl_vec_free(aff);
4989 isl_qpolynomial_free(qp);
4990 bmap = isl_basic_map_finalize(bmap);
4991 return bmap;
4992 error:
4993 isl_vec_free(aff);
4994 isl_qpolynomial_free(qp);
4995 isl_basic_map_free(bmap);
4996 return NULL;