isl_sample.c: sample_bounded: use isl_basic_set_get_ctx
[isl.git] / isl_polynomial.c
blobdd2dbf66834b3da28127972c226c1752ac569ad0
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_space_private.h>
27 #include <isl_aff_private.h>
28 #include <isl_val_private.h>
29 #include <isl_config.h>
30 #include <isl/deprecated/polynomial_int.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 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
71 __isl_keep struct isl_upoly *up2)
73 int i;
74 struct isl_upoly_rec *rec1, *rec2;
76 if (!up1 || !up2)
77 return -1;
78 if (up1 == up2)
79 return 1;
80 if (up1->var != up2->var)
81 return 0;
82 if (isl_upoly_is_cst(up1)) {
83 struct isl_upoly_cst *cst1, *cst2;
84 cst1 = isl_upoly_as_cst(up1);
85 cst2 = isl_upoly_as_cst(up2);
86 if (!cst1 || !cst2)
87 return -1;
88 return isl_int_eq(cst1->n, cst2->n) &&
89 isl_int_eq(cst1->d, cst2->d);
92 rec1 = isl_upoly_as_rec(up1);
93 rec2 = isl_upoly_as_rec(up2);
94 if (!rec1 || !rec2)
95 return -1;
97 if (rec1->n != rec2->n)
98 return 0;
100 for (i = 0; i < rec1->n; ++i) {
101 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
102 if (eq < 0 || !eq)
103 return eq;
106 return 1;
109 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
111 struct isl_upoly_cst *cst;
113 if (!up)
114 return -1;
115 if (!isl_upoly_is_cst(up))
116 return 0;
118 cst = isl_upoly_as_cst(up);
119 if (!cst)
120 return -1;
122 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
125 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
127 struct isl_upoly_cst *cst;
129 if (!up)
130 return 0;
131 if (!isl_upoly_is_cst(up))
132 return 0;
134 cst = isl_upoly_as_cst(up);
135 if (!cst)
136 return 0;
138 return isl_int_sgn(cst->n);
141 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
143 struct isl_upoly_cst *cst;
145 if (!up)
146 return -1;
147 if (!isl_upoly_is_cst(up))
148 return 0;
150 cst = isl_upoly_as_cst(up);
151 if (!cst)
152 return -1;
154 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
157 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
159 struct isl_upoly_cst *cst;
161 if (!up)
162 return -1;
163 if (!isl_upoly_is_cst(up))
164 return 0;
166 cst = isl_upoly_as_cst(up);
167 if (!cst)
168 return -1;
170 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
173 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
175 struct isl_upoly_cst *cst;
177 if (!up)
178 return -1;
179 if (!isl_upoly_is_cst(up))
180 return 0;
182 cst = isl_upoly_as_cst(up);
183 if (!cst)
184 return -1;
186 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
189 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
191 struct isl_upoly_cst *cst;
193 if (!up)
194 return -1;
195 if (!isl_upoly_is_cst(up))
196 return 0;
198 cst = isl_upoly_as_cst(up);
199 if (!cst)
200 return -1;
202 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
205 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
207 struct isl_upoly_cst *cst;
209 if (!up)
210 return -1;
211 if (!isl_upoly_is_cst(up))
212 return 0;
214 cst = isl_upoly_as_cst(up);
215 if (!cst)
216 return -1;
218 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
221 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
223 struct isl_upoly_cst *cst;
225 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
226 if (!cst)
227 return NULL;
229 cst->up.ref = 1;
230 cst->up.ctx = ctx;
231 isl_ctx_ref(ctx);
232 cst->up.var = -1;
234 isl_int_init(cst->n);
235 isl_int_init(cst->d);
237 return cst;
240 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
242 struct isl_upoly_cst *cst;
244 cst = isl_upoly_cst_alloc(ctx);
245 if (!cst)
246 return NULL;
248 isl_int_set_si(cst->n, 0);
249 isl_int_set_si(cst->d, 1);
251 return &cst->up;
254 __isl_give struct isl_upoly *isl_upoly_one(struct isl_ctx *ctx)
256 struct isl_upoly_cst *cst;
258 cst = isl_upoly_cst_alloc(ctx);
259 if (!cst)
260 return NULL;
262 isl_int_set_si(cst->n, 1);
263 isl_int_set_si(cst->d, 1);
265 return &cst->up;
268 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
270 struct isl_upoly_cst *cst;
272 cst = isl_upoly_cst_alloc(ctx);
273 if (!cst)
274 return NULL;
276 isl_int_set_si(cst->n, 1);
277 isl_int_set_si(cst->d, 0);
279 return &cst->up;
282 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
284 struct isl_upoly_cst *cst;
286 cst = isl_upoly_cst_alloc(ctx);
287 if (!cst)
288 return NULL;
290 isl_int_set_si(cst->n, -1);
291 isl_int_set_si(cst->d, 0);
293 return &cst->up;
296 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
298 struct isl_upoly_cst *cst;
300 cst = isl_upoly_cst_alloc(ctx);
301 if (!cst)
302 return NULL;
304 isl_int_set_si(cst->n, 0);
305 isl_int_set_si(cst->d, 0);
307 return &cst->up;
310 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
311 isl_int n, isl_int d)
313 struct isl_upoly_cst *cst;
315 cst = isl_upoly_cst_alloc(ctx);
316 if (!cst)
317 return NULL;
319 isl_int_set(cst->n, n);
320 isl_int_set(cst->d, d);
322 return &cst->up;
325 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
326 int var, int size)
328 struct isl_upoly_rec *rec;
330 isl_assert(ctx, var >= 0, return NULL);
331 isl_assert(ctx, size >= 0, return NULL);
332 rec = isl_calloc(ctx, struct isl_upoly_rec,
333 sizeof(struct isl_upoly_rec) +
334 size * sizeof(struct isl_upoly *));
335 if (!rec)
336 return NULL;
338 rec->up.ref = 1;
339 rec->up.ctx = ctx;
340 isl_ctx_ref(ctx);
341 rec->up.var = var;
343 rec->n = 0;
344 rec->size = size;
346 return rec;
349 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
350 __isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
352 qp = isl_qpolynomial_cow(qp);
353 if (!qp || !dim)
354 goto error;
356 isl_space_free(qp->dim);
357 qp->dim = dim;
359 return qp;
360 error:
361 isl_qpolynomial_free(qp);
362 isl_space_free(dim);
363 return NULL;
366 /* Reset the space of "qp". This function is called from isl_pw_templ.c
367 * and doesn't know if the space of an element object is represented
368 * directly or through its domain. It therefore passes along both.
370 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
371 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
372 __isl_take isl_space *domain)
374 isl_space_free(space);
375 return isl_qpolynomial_reset_domain_space(qp, domain);
378 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
380 return qp ? qp->dim->ctx : NULL;
383 __isl_give isl_space *isl_qpolynomial_get_domain_space(
384 __isl_keep isl_qpolynomial *qp)
386 return qp ? isl_space_copy(qp->dim) : NULL;
389 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
391 isl_space *space;
392 if (!qp)
393 return NULL;
394 space = isl_space_copy(qp->dim);
395 space = isl_space_from_domain(space);
396 space = isl_space_add_dims(space, isl_dim_out, 1);
397 return space;
400 /* Externally, an isl_qpolynomial has a map space, but internally, the
401 * ls field corresponds to the domain of that space.
403 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
404 enum isl_dim_type type)
406 if (!qp)
407 return 0;
408 if (type == isl_dim_out)
409 return 1;
410 if (type == isl_dim_in)
411 type = isl_dim_set;
412 return isl_space_dim(qp->dim, type);
415 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
417 return qp ? isl_upoly_is_zero(qp->upoly) : -1;
420 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
422 return qp ? isl_upoly_is_one(qp->upoly) : -1;
425 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
427 return qp ? isl_upoly_is_nan(qp->upoly) : -1;
430 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
432 return qp ? isl_upoly_is_infty(qp->upoly) : -1;
435 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
437 return qp ? isl_upoly_is_neginfty(qp->upoly) : -1;
440 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
442 return qp ? isl_upoly_sgn(qp->upoly) : 0;
445 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
447 isl_int_clear(cst->n);
448 isl_int_clear(cst->d);
451 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
453 int i;
455 for (i = 0; i < rec->n; ++i)
456 isl_upoly_free(rec->p[i]);
459 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
461 if (!up)
462 return NULL;
464 up->ref++;
465 return up;
468 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
470 struct isl_upoly_cst *cst;
471 struct isl_upoly_cst *dup;
473 cst = isl_upoly_as_cst(up);
474 if (!cst)
475 return NULL;
477 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
478 if (!dup)
479 return NULL;
480 isl_int_set(dup->n, cst->n);
481 isl_int_set(dup->d, cst->d);
483 return &dup->up;
486 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
488 int i;
489 struct isl_upoly_rec *rec;
490 struct isl_upoly_rec *dup;
492 rec = isl_upoly_as_rec(up);
493 if (!rec)
494 return NULL;
496 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
497 if (!dup)
498 return NULL;
500 for (i = 0; i < rec->n; ++i) {
501 dup->p[i] = isl_upoly_copy(rec->p[i]);
502 if (!dup->p[i])
503 goto error;
504 dup->n++;
507 return &dup->up;
508 error:
509 isl_upoly_free(&dup->up);
510 return NULL;
513 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
515 if (!up)
516 return NULL;
518 if (isl_upoly_is_cst(up))
519 return isl_upoly_dup_cst(up);
520 else
521 return isl_upoly_dup_rec(up);
524 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
526 if (!up)
527 return NULL;
529 if (up->ref == 1)
530 return up;
531 up->ref--;
532 return isl_upoly_dup(up);
535 void isl_upoly_free(__isl_take struct isl_upoly *up)
537 if (!up)
538 return;
540 if (--up->ref > 0)
541 return;
543 if (up->var < 0)
544 upoly_free_cst((struct isl_upoly_cst *)up);
545 else
546 upoly_free_rec((struct isl_upoly_rec *)up);
548 isl_ctx_deref(up->ctx);
549 free(up);
552 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
554 isl_int gcd;
556 isl_int_init(gcd);
557 isl_int_gcd(gcd, cst->n, cst->d);
558 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
559 isl_int_divexact(cst->n, cst->n, gcd);
560 isl_int_divexact(cst->d, cst->d, gcd);
562 isl_int_clear(gcd);
565 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
566 __isl_take struct isl_upoly *up2)
568 struct isl_upoly_cst *cst1;
569 struct isl_upoly_cst *cst2;
571 up1 = isl_upoly_cow(up1);
572 if (!up1 || !up2)
573 goto error;
575 cst1 = isl_upoly_as_cst(up1);
576 cst2 = isl_upoly_as_cst(up2);
578 if (isl_int_eq(cst1->d, cst2->d))
579 isl_int_add(cst1->n, cst1->n, cst2->n);
580 else {
581 isl_int_mul(cst1->n, cst1->n, cst2->d);
582 isl_int_addmul(cst1->n, cst2->n, cst1->d);
583 isl_int_mul(cst1->d, cst1->d, cst2->d);
586 isl_upoly_cst_reduce(cst1);
588 isl_upoly_free(up2);
589 return up1;
590 error:
591 isl_upoly_free(up1);
592 isl_upoly_free(up2);
593 return NULL;
596 static __isl_give struct isl_upoly *replace_by_zero(
597 __isl_take struct isl_upoly *up)
599 struct isl_ctx *ctx;
601 if (!up)
602 return NULL;
603 ctx = up->ctx;
604 isl_upoly_free(up);
605 return isl_upoly_zero(ctx);
608 static __isl_give struct isl_upoly *replace_by_constant_term(
609 __isl_take struct isl_upoly *up)
611 struct isl_upoly_rec *rec;
612 struct isl_upoly *cst;
614 if (!up)
615 return NULL;
617 rec = isl_upoly_as_rec(up);
618 if (!rec)
619 goto error;
620 cst = isl_upoly_copy(rec->p[0]);
621 isl_upoly_free(up);
622 return cst;
623 error:
624 isl_upoly_free(up);
625 return NULL;
628 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
629 __isl_take struct isl_upoly *up2)
631 int i;
632 struct isl_upoly_rec *rec1, *rec2;
634 if (!up1 || !up2)
635 goto error;
637 if (isl_upoly_is_nan(up1)) {
638 isl_upoly_free(up2);
639 return up1;
642 if (isl_upoly_is_nan(up2)) {
643 isl_upoly_free(up1);
644 return up2;
647 if (isl_upoly_is_zero(up1)) {
648 isl_upoly_free(up1);
649 return up2;
652 if (isl_upoly_is_zero(up2)) {
653 isl_upoly_free(up2);
654 return up1;
657 if (up1->var < up2->var)
658 return isl_upoly_sum(up2, up1);
660 if (up2->var < up1->var) {
661 struct isl_upoly_rec *rec;
662 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
663 isl_upoly_free(up1);
664 return up2;
666 up1 = isl_upoly_cow(up1);
667 rec = isl_upoly_as_rec(up1);
668 if (!rec)
669 goto error;
670 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
671 if (rec->n == 1)
672 up1 = replace_by_constant_term(up1);
673 return up1;
676 if (isl_upoly_is_cst(up1))
677 return isl_upoly_sum_cst(up1, up2);
679 rec1 = isl_upoly_as_rec(up1);
680 rec2 = isl_upoly_as_rec(up2);
681 if (!rec1 || !rec2)
682 goto error;
684 if (rec1->n < rec2->n)
685 return isl_upoly_sum(up2, up1);
687 up1 = isl_upoly_cow(up1);
688 rec1 = isl_upoly_as_rec(up1);
689 if (!rec1)
690 goto error;
692 for (i = rec2->n - 1; i >= 0; --i) {
693 rec1->p[i] = isl_upoly_sum(rec1->p[i],
694 isl_upoly_copy(rec2->p[i]));
695 if (!rec1->p[i])
696 goto error;
697 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
698 isl_upoly_free(rec1->p[i]);
699 rec1->n--;
703 if (rec1->n == 0)
704 up1 = replace_by_zero(up1);
705 else if (rec1->n == 1)
706 up1 = replace_by_constant_term(up1);
708 isl_upoly_free(up2);
710 return up1;
711 error:
712 isl_upoly_free(up1);
713 isl_upoly_free(up2);
714 return NULL;
717 __isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
718 __isl_take struct isl_upoly *up, isl_int v)
720 struct isl_upoly_cst *cst;
722 up = isl_upoly_cow(up);
723 if (!up)
724 return NULL;
726 cst = isl_upoly_as_cst(up);
728 isl_int_addmul(cst->n, cst->d, v);
730 return up;
733 __isl_give struct isl_upoly *isl_upoly_add_isl_int(
734 __isl_take struct isl_upoly *up, isl_int v)
736 struct isl_upoly_rec *rec;
738 if (!up)
739 return NULL;
741 if (isl_upoly_is_cst(up))
742 return isl_upoly_cst_add_isl_int(up, v);
744 up = isl_upoly_cow(up);
745 rec = isl_upoly_as_rec(up);
746 if (!rec)
747 goto error;
749 rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
750 if (!rec->p[0])
751 goto error;
753 return up;
754 error:
755 isl_upoly_free(up);
756 return NULL;
759 __isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
760 __isl_take struct isl_upoly *up, isl_int v)
762 struct isl_upoly_cst *cst;
764 if (isl_upoly_is_zero(up))
765 return up;
767 up = isl_upoly_cow(up);
768 if (!up)
769 return NULL;
771 cst = isl_upoly_as_cst(up);
773 isl_int_mul(cst->n, cst->n, v);
775 return up;
778 __isl_give struct isl_upoly *isl_upoly_mul_isl_int(
779 __isl_take struct isl_upoly *up, isl_int v)
781 int i;
782 struct isl_upoly_rec *rec;
784 if (!up)
785 return NULL;
787 if (isl_upoly_is_cst(up))
788 return isl_upoly_cst_mul_isl_int(up, v);
790 up = isl_upoly_cow(up);
791 rec = isl_upoly_as_rec(up);
792 if (!rec)
793 goto error;
795 for (i = 0; i < rec->n; ++i) {
796 rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
797 if (!rec->p[i])
798 goto error;
801 return up;
802 error:
803 isl_upoly_free(up);
804 return NULL;
807 /* Multiply the constant polynomial "up" by "v".
809 static __isl_give struct isl_upoly *isl_upoly_cst_scale_val(
810 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
812 struct isl_upoly_cst *cst;
814 if (isl_upoly_is_zero(up))
815 return up;
817 up = isl_upoly_cow(up);
818 if (!up)
819 return NULL;
821 cst = isl_upoly_as_cst(up);
823 isl_int_mul(cst->n, cst->n, v->n);
824 isl_int_mul(cst->d, cst->d, v->d);
825 isl_upoly_cst_reduce(cst);
827 return up;
830 /* Multiply the polynomial "up" by "v".
832 static __isl_give struct isl_upoly *isl_upoly_scale_val(
833 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
835 int i;
836 struct isl_upoly_rec *rec;
838 if (!up)
839 return NULL;
841 if (isl_upoly_is_cst(up))
842 return isl_upoly_cst_scale_val(up, v);
844 up = isl_upoly_cow(up);
845 rec = isl_upoly_as_rec(up);
846 if (!rec)
847 goto error;
849 for (i = 0; i < rec->n; ++i) {
850 rec->p[i] = isl_upoly_scale_val(rec->p[i], v);
851 if (!rec->p[i])
852 goto error;
855 return up;
856 error:
857 isl_upoly_free(up);
858 return NULL;
861 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
862 __isl_take struct isl_upoly *up2)
864 struct isl_upoly_cst *cst1;
865 struct isl_upoly_cst *cst2;
867 up1 = isl_upoly_cow(up1);
868 if (!up1 || !up2)
869 goto error;
871 cst1 = isl_upoly_as_cst(up1);
872 cst2 = isl_upoly_as_cst(up2);
874 isl_int_mul(cst1->n, cst1->n, cst2->n);
875 isl_int_mul(cst1->d, cst1->d, cst2->d);
877 isl_upoly_cst_reduce(cst1);
879 isl_upoly_free(up2);
880 return up1;
881 error:
882 isl_upoly_free(up1);
883 isl_upoly_free(up2);
884 return NULL;
887 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
888 __isl_take struct isl_upoly *up2)
890 struct isl_upoly_rec *rec1;
891 struct isl_upoly_rec *rec2;
892 struct isl_upoly_rec *res = NULL;
893 int i, j;
894 int size;
896 rec1 = isl_upoly_as_rec(up1);
897 rec2 = isl_upoly_as_rec(up2);
898 if (!rec1 || !rec2)
899 goto error;
900 size = rec1->n + rec2->n - 1;
901 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
902 if (!res)
903 goto error;
905 for (i = 0; i < rec1->n; ++i) {
906 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
907 isl_upoly_copy(rec1->p[i]));
908 if (!res->p[i])
909 goto error;
910 res->n++;
912 for (; i < size; ++i) {
913 res->p[i] = isl_upoly_zero(up1->ctx);
914 if (!res->p[i])
915 goto error;
916 res->n++;
918 for (i = 0; i < rec1->n; ++i) {
919 for (j = 1; j < rec2->n; ++j) {
920 struct isl_upoly *up;
921 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
922 isl_upoly_copy(rec1->p[i]));
923 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
924 if (!res->p[i + j])
925 goto error;
929 isl_upoly_free(up1);
930 isl_upoly_free(up2);
932 return &res->up;
933 error:
934 isl_upoly_free(up1);
935 isl_upoly_free(up2);
936 isl_upoly_free(&res->up);
937 return NULL;
940 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
941 __isl_take struct isl_upoly *up2)
943 if (!up1 || !up2)
944 goto error;
946 if (isl_upoly_is_nan(up1)) {
947 isl_upoly_free(up2);
948 return up1;
951 if (isl_upoly_is_nan(up2)) {
952 isl_upoly_free(up1);
953 return up2;
956 if (isl_upoly_is_zero(up1)) {
957 isl_upoly_free(up2);
958 return up1;
961 if (isl_upoly_is_zero(up2)) {
962 isl_upoly_free(up1);
963 return up2;
966 if (isl_upoly_is_one(up1)) {
967 isl_upoly_free(up1);
968 return up2;
971 if (isl_upoly_is_one(up2)) {
972 isl_upoly_free(up2);
973 return up1;
976 if (up1->var < up2->var)
977 return isl_upoly_mul(up2, up1);
979 if (up2->var < up1->var) {
980 int i;
981 struct isl_upoly_rec *rec;
982 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
983 isl_ctx *ctx = up1->ctx;
984 isl_upoly_free(up1);
985 isl_upoly_free(up2);
986 return isl_upoly_nan(ctx);
988 up1 = isl_upoly_cow(up1);
989 rec = isl_upoly_as_rec(up1);
990 if (!rec)
991 goto error;
993 for (i = 0; i < rec->n; ++i) {
994 rec->p[i] = isl_upoly_mul(rec->p[i],
995 isl_upoly_copy(up2));
996 if (!rec->p[i])
997 goto error;
999 isl_upoly_free(up2);
1000 return up1;
1003 if (isl_upoly_is_cst(up1))
1004 return isl_upoly_mul_cst(up1, up2);
1006 return isl_upoly_mul_rec(up1, up2);
1007 error:
1008 isl_upoly_free(up1);
1009 isl_upoly_free(up2);
1010 return NULL;
1013 __isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
1014 unsigned power)
1016 struct isl_upoly *res;
1018 if (!up)
1019 return NULL;
1020 if (power == 1)
1021 return up;
1023 if (power % 2)
1024 res = isl_upoly_copy(up);
1025 else
1026 res = isl_upoly_one(up->ctx);
1028 while (power >>= 1) {
1029 up = isl_upoly_mul(up, isl_upoly_copy(up));
1030 if (power % 2)
1031 res = isl_upoly_mul(res, isl_upoly_copy(up));
1034 isl_upoly_free(up);
1035 return res;
1038 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
1039 unsigned n_div, __isl_take struct isl_upoly *up)
1041 struct isl_qpolynomial *qp = NULL;
1042 unsigned total;
1044 if (!dim || !up)
1045 goto error;
1047 if (!isl_space_is_set(dim))
1048 isl_die(isl_space_get_ctx(dim), isl_error_invalid,
1049 "domain of polynomial should be a set", goto error);
1051 total = isl_space_dim(dim, isl_dim_all);
1053 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
1054 if (!qp)
1055 goto error;
1057 qp->ref = 1;
1058 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
1059 if (!qp->div)
1060 goto error;
1062 qp->dim = dim;
1063 qp->upoly = up;
1065 return qp;
1066 error:
1067 isl_space_free(dim);
1068 isl_upoly_free(up);
1069 isl_qpolynomial_free(qp);
1070 return NULL;
1073 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1075 if (!qp)
1076 return NULL;
1078 qp->ref++;
1079 return qp;
1082 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1084 struct isl_qpolynomial *dup;
1086 if (!qp)
1087 return NULL;
1089 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1090 isl_upoly_copy(qp->upoly));
1091 if (!dup)
1092 return NULL;
1093 isl_mat_free(dup->div);
1094 dup->div = isl_mat_copy(qp->div);
1095 if (!dup->div)
1096 goto error;
1098 return dup;
1099 error:
1100 isl_qpolynomial_free(dup);
1101 return NULL;
1104 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1106 if (!qp)
1107 return NULL;
1109 if (qp->ref == 1)
1110 return qp;
1111 qp->ref--;
1112 return isl_qpolynomial_dup(qp);
1115 __isl_null isl_qpolynomial *isl_qpolynomial_free(
1116 __isl_take isl_qpolynomial *qp)
1118 if (!qp)
1119 return NULL;
1121 if (--qp->ref > 0)
1122 return NULL;
1124 isl_space_free(qp->dim);
1125 isl_mat_free(qp->div);
1126 isl_upoly_free(qp->upoly);
1128 free(qp);
1129 return NULL;
1132 __isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
1134 int i;
1135 struct isl_upoly_rec *rec;
1136 struct isl_upoly_cst *cst;
1138 rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1139 if (!rec)
1140 return NULL;
1141 for (i = 0; i < 1 + power; ++i) {
1142 rec->p[i] = isl_upoly_zero(ctx);
1143 if (!rec->p[i])
1144 goto error;
1145 rec->n++;
1147 cst = isl_upoly_as_cst(rec->p[power]);
1148 isl_int_set_si(cst->n, 1);
1150 return &rec->up;
1151 error:
1152 isl_upoly_free(&rec->up);
1153 return NULL;
1156 /* r array maps original positions to new positions.
1158 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
1159 int *r)
1161 int i;
1162 struct isl_upoly_rec *rec;
1163 struct isl_upoly *base;
1164 struct isl_upoly *res;
1166 if (isl_upoly_is_cst(up))
1167 return up;
1169 rec = isl_upoly_as_rec(up);
1170 if (!rec)
1171 goto error;
1173 isl_assert(up->ctx, rec->n >= 1, goto error);
1175 base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
1176 res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
1178 for (i = rec->n - 2; i >= 0; --i) {
1179 res = isl_upoly_mul(res, isl_upoly_copy(base));
1180 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
1183 isl_upoly_free(base);
1184 isl_upoly_free(up);
1186 return res;
1187 error:
1188 isl_upoly_free(up);
1189 return NULL;
1192 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
1194 int n_row, n_col;
1195 int equal;
1197 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1198 div1->n_col >= div2->n_col, return -1);
1200 if (div1->n_row == div2->n_row)
1201 return isl_mat_is_equal(div1, div2);
1203 n_row = div1->n_row;
1204 n_col = div1->n_col;
1205 div1->n_row = div2->n_row;
1206 div1->n_col = div2->n_col;
1208 equal = isl_mat_is_equal(div1, div2);
1210 div1->n_row = n_row;
1211 div1->n_col = n_col;
1213 return equal;
1216 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1218 int li, lj;
1220 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1221 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1223 if (li != lj)
1224 return li - lj;
1226 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1229 struct isl_div_sort_info {
1230 isl_mat *div;
1231 int row;
1234 static int div_sort_cmp(const void *p1, const void *p2)
1236 const struct isl_div_sort_info *i1, *i2;
1237 i1 = (const struct isl_div_sort_info *) p1;
1238 i2 = (const struct isl_div_sort_info *) p2;
1240 return cmp_row(i1->div, i1->row, i2->row);
1243 /* Sort divs and remove duplicates.
1245 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1247 int i;
1248 int skip;
1249 int len;
1250 struct isl_div_sort_info *array = NULL;
1251 int *pos = NULL, *at = NULL;
1252 int *reordering = NULL;
1253 unsigned div_pos;
1255 if (!qp)
1256 return NULL;
1257 if (qp->div->n_row <= 1)
1258 return qp;
1260 div_pos = isl_space_dim(qp->dim, isl_dim_all);
1262 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1263 qp->div->n_row);
1264 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1265 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1266 len = qp->div->n_col - 2;
1267 reordering = isl_alloc_array(qp->div->ctx, int, len);
1268 if (!array || !pos || !at || !reordering)
1269 goto error;
1271 for (i = 0; i < qp->div->n_row; ++i) {
1272 array[i].div = qp->div;
1273 array[i].row = i;
1274 pos[i] = i;
1275 at[i] = i;
1278 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1279 div_sort_cmp);
1281 for (i = 0; i < div_pos; ++i)
1282 reordering[i] = i;
1284 for (i = 0; i < qp->div->n_row; ++i) {
1285 if (pos[array[i].row] == i)
1286 continue;
1287 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1288 pos[at[i]] = pos[array[i].row];
1289 at[pos[array[i].row]] = at[i];
1290 at[i] = array[i].row;
1291 pos[array[i].row] = i;
1294 skip = 0;
1295 for (i = 0; i < len - div_pos; ++i) {
1296 if (i > 0 &&
1297 isl_seq_eq(qp->div->row[i - skip - 1],
1298 qp->div->row[i - skip], qp->div->n_col)) {
1299 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1300 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1301 2 + div_pos + i - skip);
1302 qp->div = isl_mat_drop_cols(qp->div,
1303 2 + div_pos + i - skip, 1);
1304 skip++;
1306 reordering[div_pos + array[i].row] = div_pos + i - skip;
1309 qp->upoly = reorder(qp->upoly, reordering);
1311 if (!qp->upoly || !qp->div)
1312 goto error;
1314 free(at);
1315 free(pos);
1316 free(array);
1317 free(reordering);
1319 return qp;
1320 error:
1321 free(at);
1322 free(pos);
1323 free(array);
1324 free(reordering);
1325 isl_qpolynomial_free(qp);
1326 return NULL;
1329 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1330 int *exp, int first)
1332 int i;
1333 struct isl_upoly_rec *rec;
1335 if (isl_upoly_is_cst(up))
1336 return up;
1338 if (up->var < first)
1339 return up;
1341 if (exp[up->var - first] == up->var - first)
1342 return up;
1344 up = isl_upoly_cow(up);
1345 if (!up)
1346 goto error;
1348 up->var = exp[up->var - first] + first;
1350 rec = isl_upoly_as_rec(up);
1351 if (!rec)
1352 goto error;
1354 for (i = 0; i < rec->n; ++i) {
1355 rec->p[i] = expand(rec->p[i], exp, first);
1356 if (!rec->p[i])
1357 goto error;
1360 return up;
1361 error:
1362 isl_upoly_free(up);
1363 return NULL;
1366 static __isl_give isl_qpolynomial *with_merged_divs(
1367 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1368 __isl_take isl_qpolynomial *qp2),
1369 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1371 int *exp1 = NULL;
1372 int *exp2 = NULL;
1373 isl_mat *div = NULL;
1374 int n_div1, n_div2;
1376 qp1 = isl_qpolynomial_cow(qp1);
1377 qp2 = isl_qpolynomial_cow(qp2);
1379 if (!qp1 || !qp2)
1380 goto error;
1382 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1383 qp1->div->n_col >= qp2->div->n_col, goto error);
1385 n_div1 = qp1->div->n_row;
1386 n_div2 = qp2->div->n_row;
1387 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1);
1388 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2);
1389 if ((n_div1 && !exp1) || (n_div2 && !exp2))
1390 goto error;
1392 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1393 if (!div)
1394 goto error;
1396 isl_mat_free(qp1->div);
1397 qp1->div = isl_mat_copy(div);
1398 isl_mat_free(qp2->div);
1399 qp2->div = isl_mat_copy(div);
1401 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1402 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1404 if (!qp1->upoly || !qp2->upoly)
1405 goto error;
1407 isl_mat_free(div);
1408 free(exp1);
1409 free(exp2);
1411 return fn(qp1, qp2);
1412 error:
1413 isl_mat_free(div);
1414 free(exp1);
1415 free(exp2);
1416 isl_qpolynomial_free(qp1);
1417 isl_qpolynomial_free(qp2);
1418 return NULL;
1421 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1422 __isl_take isl_qpolynomial *qp2)
1424 qp1 = isl_qpolynomial_cow(qp1);
1426 if (!qp1 || !qp2)
1427 goto error;
1429 if (qp1->div->n_row < qp2->div->n_row)
1430 return isl_qpolynomial_add(qp2, qp1);
1432 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1433 if (!compatible_divs(qp1->div, qp2->div))
1434 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1436 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1437 if (!qp1->upoly)
1438 goto error;
1440 isl_qpolynomial_free(qp2);
1442 return qp1;
1443 error:
1444 isl_qpolynomial_free(qp1);
1445 isl_qpolynomial_free(qp2);
1446 return NULL;
1449 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1450 __isl_keep isl_set *dom,
1451 __isl_take isl_qpolynomial *qp1,
1452 __isl_take isl_qpolynomial *qp2)
1454 qp1 = isl_qpolynomial_add(qp1, qp2);
1455 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1456 return qp1;
1459 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1460 __isl_take isl_qpolynomial *qp2)
1462 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1465 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1466 __isl_take isl_qpolynomial *qp, isl_int v)
1468 if (isl_int_is_zero(v))
1469 return qp;
1471 qp = isl_qpolynomial_cow(qp);
1472 if (!qp)
1473 return NULL;
1475 qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1476 if (!qp->upoly)
1477 goto error;
1479 return qp;
1480 error:
1481 isl_qpolynomial_free(qp);
1482 return NULL;
1486 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1488 if (!qp)
1489 return NULL;
1491 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1494 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1495 __isl_take isl_qpolynomial *qp, isl_int v)
1497 if (isl_int_is_one(v))
1498 return qp;
1500 if (qp && isl_int_is_zero(v)) {
1501 isl_qpolynomial *zero;
1502 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1503 isl_qpolynomial_free(qp);
1504 return zero;
1507 qp = isl_qpolynomial_cow(qp);
1508 if (!qp)
1509 return NULL;
1511 qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1512 if (!qp->upoly)
1513 goto error;
1515 return qp;
1516 error:
1517 isl_qpolynomial_free(qp);
1518 return NULL;
1521 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1522 __isl_take isl_qpolynomial *qp, isl_int v)
1524 return isl_qpolynomial_mul_isl_int(qp, v);
1527 /* Multiply "qp" by "v".
1529 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1530 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1532 if (!qp || !v)
1533 goto error;
1535 if (!isl_val_is_rat(v))
1536 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1537 "expecting rational factor", goto error);
1539 if (isl_val_is_one(v)) {
1540 isl_val_free(v);
1541 return qp;
1544 if (isl_val_is_zero(v)) {
1545 isl_space *space;
1547 space = isl_qpolynomial_get_domain_space(qp);
1548 isl_qpolynomial_free(qp);
1549 isl_val_free(v);
1550 return isl_qpolynomial_zero_on_domain(space);
1553 qp = isl_qpolynomial_cow(qp);
1554 if (!qp)
1555 goto error;
1557 qp->upoly = isl_upoly_scale_val(qp->upoly, v);
1558 if (!qp->upoly)
1559 qp = isl_qpolynomial_free(qp);
1561 isl_val_free(v);
1562 return qp;
1563 error:
1564 isl_val_free(v);
1565 isl_qpolynomial_free(qp);
1566 return NULL;
1569 /* Divide "qp" by "v".
1571 __isl_give isl_qpolynomial *isl_qpolynomial_scale_down_val(
1572 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1574 if (!qp || !v)
1575 goto error;
1577 if (!isl_val_is_rat(v))
1578 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1579 "expecting rational factor", goto error);
1580 if (isl_val_is_zero(v))
1581 isl_die(isl_val_get_ctx(v), isl_error_invalid,
1582 "cannot scale down by zero", goto error);
1584 return isl_qpolynomial_scale_val(qp, isl_val_inv(v));
1585 error:
1586 isl_val_free(v);
1587 isl_qpolynomial_free(qp);
1588 return NULL;
1591 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1592 __isl_take isl_qpolynomial *qp2)
1594 qp1 = isl_qpolynomial_cow(qp1);
1596 if (!qp1 || !qp2)
1597 goto error;
1599 if (qp1->div->n_row < qp2->div->n_row)
1600 return isl_qpolynomial_mul(qp2, qp1);
1602 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1603 if (!compatible_divs(qp1->div, qp2->div))
1604 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1606 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1607 if (!qp1->upoly)
1608 goto error;
1610 isl_qpolynomial_free(qp2);
1612 return qp1;
1613 error:
1614 isl_qpolynomial_free(qp1);
1615 isl_qpolynomial_free(qp2);
1616 return NULL;
1619 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1620 unsigned power)
1622 qp = isl_qpolynomial_cow(qp);
1624 if (!qp)
1625 return NULL;
1627 qp->upoly = isl_upoly_pow(qp->upoly, power);
1628 if (!qp->upoly)
1629 goto error;
1631 return qp;
1632 error:
1633 isl_qpolynomial_free(qp);
1634 return NULL;
1637 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1638 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1640 int i;
1642 if (power == 1)
1643 return pwqp;
1645 pwqp = isl_pw_qpolynomial_cow(pwqp);
1646 if (!pwqp)
1647 return NULL;
1649 for (i = 0; i < pwqp->n; ++i) {
1650 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1651 if (!pwqp->p[i].qp)
1652 return isl_pw_qpolynomial_free(pwqp);
1655 return pwqp;
1658 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1659 __isl_take isl_space *dim)
1661 if (!dim)
1662 return NULL;
1663 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1666 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1667 __isl_take isl_space *dim)
1669 if (!dim)
1670 return NULL;
1671 return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1674 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1675 __isl_take isl_space *dim)
1677 if (!dim)
1678 return NULL;
1679 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1682 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1683 __isl_take isl_space *dim)
1685 if (!dim)
1686 return NULL;
1687 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1690 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1691 __isl_take isl_space *dim)
1693 if (!dim)
1694 return NULL;
1695 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1698 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1699 __isl_take isl_space *dim,
1700 isl_int v)
1702 struct isl_qpolynomial *qp;
1703 struct isl_upoly_cst *cst;
1705 if (!dim)
1706 return NULL;
1708 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1709 if (!qp)
1710 return NULL;
1712 cst = isl_upoly_as_cst(qp->upoly);
1713 isl_int_set(cst->n, v);
1715 return qp;
1718 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1719 isl_int *n, isl_int *d)
1721 struct isl_upoly_cst *cst;
1723 if (!qp)
1724 return -1;
1726 if (!isl_upoly_is_cst(qp->upoly))
1727 return 0;
1729 cst = isl_upoly_as_cst(qp->upoly);
1730 if (!cst)
1731 return -1;
1733 if (n)
1734 isl_int_set(*n, cst->n);
1735 if (d)
1736 isl_int_set(*d, cst->d);
1738 return 1;
1741 /* Return the constant term of "up".
1743 static __isl_give isl_val *isl_upoly_get_constant_val(
1744 __isl_keep struct isl_upoly *up)
1746 struct isl_upoly_cst *cst;
1748 if (!up)
1749 return NULL;
1751 while (!isl_upoly_is_cst(up)) {
1752 struct isl_upoly_rec *rec;
1754 rec = isl_upoly_as_rec(up);
1755 if (!rec)
1756 return NULL;
1757 up = rec->p[0];
1760 cst = isl_upoly_as_cst(up);
1761 if (!cst)
1762 return NULL;
1763 return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1766 /* Return the constant term of "qp".
1768 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1769 __isl_keep isl_qpolynomial *qp)
1771 if (!qp)
1772 return NULL;
1774 return isl_upoly_get_constant_val(qp->upoly);
1777 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1779 int is_cst;
1780 struct isl_upoly_rec *rec;
1782 if (!up)
1783 return -1;
1785 if (up->var < 0)
1786 return 1;
1788 rec = isl_upoly_as_rec(up);
1789 if (!rec)
1790 return -1;
1792 if (rec->n > 2)
1793 return 0;
1795 isl_assert(up->ctx, rec->n > 1, return -1);
1797 is_cst = isl_upoly_is_cst(rec->p[1]);
1798 if (is_cst < 0)
1799 return -1;
1800 if (!is_cst)
1801 return 0;
1803 return isl_upoly_is_affine(rec->p[0]);
1806 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1808 if (!qp)
1809 return -1;
1811 if (qp->div->n_row > 0)
1812 return 0;
1814 return isl_upoly_is_affine(qp->upoly);
1817 static void update_coeff(__isl_keep isl_vec *aff,
1818 __isl_keep struct isl_upoly_cst *cst, int pos)
1820 isl_int gcd;
1821 isl_int f;
1823 if (isl_int_is_zero(cst->n))
1824 return;
1826 isl_int_init(gcd);
1827 isl_int_init(f);
1828 isl_int_gcd(gcd, cst->d, aff->el[0]);
1829 isl_int_divexact(f, cst->d, gcd);
1830 isl_int_divexact(gcd, aff->el[0], gcd);
1831 isl_seq_scale(aff->el, aff->el, f, aff->size);
1832 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1833 isl_int_clear(gcd);
1834 isl_int_clear(f);
1837 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1838 __isl_keep isl_vec *aff)
1840 struct isl_upoly_cst *cst;
1841 struct isl_upoly_rec *rec;
1843 if (!up || !aff)
1844 return -1;
1846 if (up->var < 0) {
1847 struct isl_upoly_cst *cst;
1849 cst = isl_upoly_as_cst(up);
1850 if (!cst)
1851 return -1;
1852 update_coeff(aff, cst, 0);
1853 return 0;
1856 rec = isl_upoly_as_rec(up);
1857 if (!rec)
1858 return -1;
1859 isl_assert(up->ctx, rec->n == 2, return -1);
1861 cst = isl_upoly_as_cst(rec->p[1]);
1862 if (!cst)
1863 return -1;
1864 update_coeff(aff, cst, 1 + up->var);
1866 return isl_upoly_update_affine(rec->p[0], aff);
1869 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1870 __isl_keep isl_qpolynomial *qp)
1872 isl_vec *aff;
1873 unsigned d;
1875 if (!qp)
1876 return NULL;
1878 d = isl_space_dim(qp->dim, isl_dim_all);
1879 aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1880 if (!aff)
1881 return NULL;
1883 isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1884 isl_int_set_si(aff->el[0], 1);
1886 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1887 goto error;
1889 return aff;
1890 error:
1891 isl_vec_free(aff);
1892 return NULL;
1895 int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
1896 __isl_keep isl_qpolynomial *qp2)
1898 int equal;
1900 if (!qp1 || !qp2)
1901 return -1;
1903 equal = isl_space_is_equal(qp1->dim, qp2->dim);
1904 if (equal < 0 || !equal)
1905 return equal;
1907 equal = isl_mat_is_equal(qp1->div, qp2->div);
1908 if (equal < 0 || !equal)
1909 return equal;
1911 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1914 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1916 int i;
1917 struct isl_upoly_rec *rec;
1919 if (isl_upoly_is_cst(up)) {
1920 struct isl_upoly_cst *cst;
1921 cst = isl_upoly_as_cst(up);
1922 if (!cst)
1923 return;
1924 isl_int_lcm(*d, *d, cst->d);
1925 return;
1928 rec = isl_upoly_as_rec(up);
1929 if (!rec)
1930 return;
1932 for (i = 0; i < rec->n; ++i)
1933 upoly_update_den(rec->p[i], d);
1936 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1938 isl_int_set_si(*d, 1);
1939 if (!qp)
1940 return;
1941 upoly_update_den(qp->upoly, d);
1944 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
1945 __isl_take isl_space *dim, int pos, int power)
1947 struct isl_ctx *ctx;
1949 if (!dim)
1950 return NULL;
1952 ctx = dim->ctx;
1954 return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
1957 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
1958 enum isl_dim_type type, unsigned pos)
1960 if (!dim)
1961 return NULL;
1963 isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
1964 isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
1966 if (type == isl_dim_set)
1967 pos += isl_space_dim(dim, isl_dim_param);
1969 return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
1970 error:
1971 isl_space_free(dim);
1972 return NULL;
1975 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1976 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1978 int i;
1979 struct isl_upoly_rec *rec;
1980 struct isl_upoly *base, *res;
1982 if (!up)
1983 return NULL;
1985 if (isl_upoly_is_cst(up))
1986 return up;
1988 if (up->var < first)
1989 return up;
1991 rec = isl_upoly_as_rec(up);
1992 if (!rec)
1993 goto error;
1995 isl_assert(up->ctx, rec->n >= 1, goto error);
1997 if (up->var >= first + n)
1998 base = isl_upoly_var_pow(up->ctx, up->var, 1);
1999 else
2000 base = isl_upoly_copy(subs[up->var - first]);
2002 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
2003 for (i = rec->n - 2; i >= 0; --i) {
2004 struct isl_upoly *t;
2005 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
2006 res = isl_upoly_mul(res, isl_upoly_copy(base));
2007 res = isl_upoly_sum(res, t);
2010 isl_upoly_free(base);
2011 isl_upoly_free(up);
2013 return res;
2014 error:
2015 isl_upoly_free(up);
2016 return NULL;
2019 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
2020 isl_int denom, unsigned len)
2022 int i;
2023 struct isl_upoly *up;
2025 isl_assert(ctx, len >= 1, return NULL);
2027 up = isl_upoly_rat_cst(ctx, f[0], denom);
2028 for (i = 0; i < len - 1; ++i) {
2029 struct isl_upoly *t;
2030 struct isl_upoly *c;
2032 if (isl_int_is_zero(f[1 + i]))
2033 continue;
2035 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
2036 t = isl_upoly_var_pow(ctx, i, 1);
2037 t = isl_upoly_mul(c, t);
2038 up = isl_upoly_sum(up, t);
2041 return up;
2044 /* Remove common factor of non-constant terms and denominator.
2046 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2048 isl_ctx *ctx = qp->div->ctx;
2049 unsigned total = qp->div->n_col - 2;
2051 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2052 isl_int_gcd(ctx->normalize_gcd,
2053 ctx->normalize_gcd, qp->div->row[div][0]);
2054 if (isl_int_is_one(ctx->normalize_gcd))
2055 return;
2057 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2058 ctx->normalize_gcd, total);
2059 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2060 ctx->normalize_gcd);
2061 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2062 ctx->normalize_gcd);
2065 /* Replace the integer division identified by "div" by the polynomial "s".
2066 * The integer division is assumed not to appear in the definition
2067 * of any other integer divisions.
2069 static __isl_give isl_qpolynomial *substitute_div(
2070 __isl_take isl_qpolynomial *qp,
2071 int div, __isl_take struct isl_upoly *s)
2073 int i;
2074 int total;
2075 int *reordering;
2077 if (!qp || !s)
2078 goto error;
2080 qp = isl_qpolynomial_cow(qp);
2081 if (!qp)
2082 goto error;
2084 total = isl_space_dim(qp->dim, isl_dim_all);
2085 qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
2086 if (!qp->upoly)
2087 goto error;
2089 reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
2090 if (!reordering)
2091 goto error;
2092 for (i = 0; i < total + div; ++i)
2093 reordering[i] = i;
2094 for (i = total + div + 1; i < total + qp->div->n_row; ++i)
2095 reordering[i] = i - 1;
2096 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2097 qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
2098 qp->upoly = reorder(qp->upoly, reordering);
2099 free(reordering);
2101 if (!qp->upoly || !qp->div)
2102 goto error;
2104 isl_upoly_free(s);
2105 return qp;
2106 error:
2107 isl_qpolynomial_free(qp);
2108 isl_upoly_free(s);
2109 return NULL;
2112 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2113 * divisions because d is equal to 1 by their definition, i.e., e.
2115 static __isl_give isl_qpolynomial *substitute_non_divs(
2116 __isl_take isl_qpolynomial *qp)
2118 int i, j;
2119 int total;
2120 struct isl_upoly *s;
2122 if (!qp)
2123 return NULL;
2125 total = isl_space_dim(qp->dim, isl_dim_all);
2126 for (i = 0; qp && i < qp->div->n_row; ++i) {
2127 if (!isl_int_is_one(qp->div->row[i][0]))
2128 continue;
2129 for (j = i + 1; j < qp->div->n_row; ++j) {
2130 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2131 continue;
2132 isl_seq_combine(qp->div->row[j] + 1,
2133 qp->div->ctx->one, qp->div->row[j] + 1,
2134 qp->div->row[j][2 + total + i],
2135 qp->div->row[i] + 1, 1 + total + i);
2136 isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2137 normalize_div(qp, j);
2139 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2140 qp->div->row[i][0], qp->div->n_col - 1);
2141 qp = substitute_div(qp, i, s);
2142 --i;
2145 return qp;
2148 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2149 * with d the denominator. When replacing the coefficient e of x by
2150 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2151 * inside the division, so we need to add floor(e/d) * x outside.
2152 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2153 * to adjust the coefficient of x in each later div that depends on the
2154 * current div "div" and also in the affine expression "aff"
2155 * (if it too depends on "div").
2157 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2158 __isl_keep isl_vec *aff)
2160 int i, j;
2161 isl_int v;
2162 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2164 isl_int_init(v);
2165 for (i = 0; i < 1 + total + div; ++i) {
2166 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2167 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2168 continue;
2169 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2170 isl_int_fdiv_r(qp->div->row[div][1 + i],
2171 qp->div->row[div][1 + i], qp->div->row[div][0]);
2172 if (!isl_int_is_zero(aff->el[1 + total + div]))
2173 isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
2174 for (j = div + 1; j < qp->div->n_row; ++j) {
2175 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2176 continue;
2177 isl_int_addmul(qp->div->row[j][1 + i],
2178 v, qp->div->row[j][2 + total + div]);
2181 isl_int_clear(v);
2184 /* Check if the last non-zero coefficient is bigger that half of the
2185 * denominator. If so, we will invert the div to further reduce the number
2186 * of distinct divs that may appear.
2187 * If the last non-zero coefficient is exactly half the denominator,
2188 * then we continue looking for earlier coefficients that are bigger
2189 * than half the denominator.
2191 static int needs_invert(__isl_keep isl_mat *div, int row)
2193 int i;
2194 int cmp;
2196 for (i = div->n_col - 1; i >= 1; --i) {
2197 if (isl_int_is_zero(div->row[row][i]))
2198 continue;
2199 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2200 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2201 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2202 if (cmp)
2203 return cmp > 0;
2204 if (i == 1)
2205 return 1;
2208 return 0;
2211 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2212 * We only invert the coefficients of e (and the coefficient of q in
2213 * later divs and in "aff"). After calling this function, the
2214 * coefficients of e should be reduced again.
2216 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2217 __isl_keep isl_vec *aff)
2219 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2221 isl_seq_neg(qp->div->row[div] + 1,
2222 qp->div->row[div] + 1, qp->div->n_col - 1);
2223 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2224 isl_int_add(qp->div->row[div][1],
2225 qp->div->row[div][1], qp->div->row[div][0]);
2226 if (!isl_int_is_zero(aff->el[1 + total + div]))
2227 isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
2228 isl_mat_col_mul(qp->div, 2 + total + div,
2229 qp->div->ctx->negone, 2 + total + div);
2232 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2233 * in the interval [0, d-1], with d the denominator and such that the
2234 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2236 * After the reduction, some divs may have become redundant or identical,
2237 * so we call substitute_non_divs and sort_divs. If these functions
2238 * eliminate divs or merge two or more divs into one, the coefficients
2239 * of the enclosing divs may have to be reduced again, so we call
2240 * ourselves recursively if the number of divs decreases.
2242 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2244 int i;
2245 isl_vec *aff = NULL;
2246 struct isl_upoly *s;
2247 unsigned n_div;
2249 if (!qp)
2250 return NULL;
2252 aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
2253 aff = isl_vec_clr(aff);
2254 if (!aff)
2255 goto error;
2257 isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
2259 for (i = 0; i < qp->div->n_row; ++i) {
2260 normalize_div(qp, i);
2261 reduce_div(qp, i, aff);
2262 if (needs_invert(qp->div, i)) {
2263 invert_div(qp, i, aff);
2264 reduce_div(qp, i, aff);
2268 s = isl_upoly_from_affine(qp->div->ctx, aff->el,
2269 qp->div->ctx->one, aff->size);
2270 qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
2271 isl_upoly_free(s);
2272 if (!qp->upoly)
2273 goto error;
2275 isl_vec_free(aff);
2277 n_div = qp->div->n_row;
2278 qp = substitute_non_divs(qp);
2279 qp = sort_divs(qp);
2280 if (qp && qp->div->n_row < n_div)
2281 return reduce_divs(qp);
2283 return qp;
2284 error:
2285 isl_qpolynomial_free(qp);
2286 isl_vec_free(aff);
2287 return NULL;
2290 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2291 __isl_take isl_space *dim, const isl_int n, const isl_int d)
2293 struct isl_qpolynomial *qp;
2294 struct isl_upoly_cst *cst;
2296 if (!dim)
2297 return NULL;
2299 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2300 if (!qp)
2301 return NULL;
2303 cst = isl_upoly_as_cst(qp->upoly);
2304 isl_int_set(cst->n, n);
2305 isl_int_set(cst->d, d);
2307 return qp;
2310 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2312 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2313 __isl_take isl_space *domain, __isl_take isl_val *val)
2315 isl_qpolynomial *qp;
2316 struct isl_upoly_cst *cst;
2318 if (!domain || !val)
2319 goto error;
2321 qp = isl_qpolynomial_alloc(isl_space_copy(domain), 0,
2322 isl_upoly_zero(domain->ctx));
2323 if (!qp)
2324 goto error;
2326 cst = isl_upoly_as_cst(qp->upoly);
2327 isl_int_set(cst->n, val->n);
2328 isl_int_set(cst->d, val->d);
2330 isl_space_free(domain);
2331 isl_val_free(val);
2332 return qp;
2333 error:
2334 isl_space_free(domain);
2335 isl_val_free(val);
2336 return NULL;
2339 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2341 struct isl_upoly_rec *rec;
2342 int i;
2344 if (!up)
2345 return -1;
2347 if (isl_upoly_is_cst(up))
2348 return 0;
2350 if (up->var < d)
2351 active[up->var] = 1;
2353 rec = isl_upoly_as_rec(up);
2354 for (i = 0; i < rec->n; ++i)
2355 if (up_set_active(rec->p[i], active, d) < 0)
2356 return -1;
2358 return 0;
2361 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2363 int i, j;
2364 int d = isl_space_dim(qp->dim, isl_dim_all);
2366 if (!qp || !active)
2367 return -1;
2369 for (i = 0; i < d; ++i)
2370 for (j = 0; j < qp->div->n_row; ++j) {
2371 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2372 continue;
2373 active[i] = 1;
2374 break;
2377 return up_set_active(qp->upoly, active, d);
2380 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2381 enum isl_dim_type type, unsigned first, unsigned n)
2383 int i;
2384 int *active = NULL;
2385 int involves = 0;
2387 if (!qp)
2388 return -1;
2389 if (n == 0)
2390 return 0;
2392 isl_assert(qp->dim->ctx,
2393 first + n <= isl_qpolynomial_dim(qp, type), return -1);
2394 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2395 type == isl_dim_in, return -1);
2397 active = isl_calloc_array(qp->dim->ctx, int,
2398 isl_space_dim(qp->dim, isl_dim_all));
2399 if (set_active(qp, active) < 0)
2400 goto error;
2402 if (type == isl_dim_in)
2403 first += isl_space_dim(qp->dim, isl_dim_param);
2404 for (i = 0; i < n; ++i)
2405 if (active[first + i]) {
2406 involves = 1;
2407 break;
2410 free(active);
2412 return involves;
2413 error:
2414 free(active);
2415 return -1;
2418 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2419 * of the divs that do appear in the quasi-polynomial.
2421 static __isl_give isl_qpolynomial *remove_redundant_divs(
2422 __isl_take isl_qpolynomial *qp)
2424 int i, j;
2425 int d;
2426 int len;
2427 int skip;
2428 int *active = NULL;
2429 int *reordering = NULL;
2430 int redundant = 0;
2431 int n_div;
2432 isl_ctx *ctx;
2434 if (!qp)
2435 return NULL;
2436 if (qp->div->n_row == 0)
2437 return qp;
2439 d = isl_space_dim(qp->dim, isl_dim_all);
2440 len = qp->div->n_col - 2;
2441 ctx = isl_qpolynomial_get_ctx(qp);
2442 active = isl_calloc_array(ctx, int, len);
2443 if (!active)
2444 goto error;
2446 if (up_set_active(qp->upoly, active, len) < 0)
2447 goto error;
2449 for (i = qp->div->n_row - 1; i >= 0; --i) {
2450 if (!active[d + i]) {
2451 redundant = 1;
2452 continue;
2454 for (j = 0; j < i; ++j) {
2455 if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2456 continue;
2457 active[d + j] = 1;
2458 break;
2462 if (!redundant) {
2463 free(active);
2464 return qp;
2467 reordering = isl_alloc_array(qp->div->ctx, int, len);
2468 if (!reordering)
2469 goto error;
2471 for (i = 0; i < d; ++i)
2472 reordering[i] = i;
2474 skip = 0;
2475 n_div = qp->div->n_row;
2476 for (i = 0; i < n_div; ++i) {
2477 if (!active[d + i]) {
2478 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2479 qp->div = isl_mat_drop_cols(qp->div,
2480 2 + d + i - skip, 1);
2481 skip++;
2483 reordering[d + i] = d + i - skip;
2486 qp->upoly = reorder(qp->upoly, reordering);
2488 if (!qp->upoly || !qp->div)
2489 goto error;
2491 free(active);
2492 free(reordering);
2494 return qp;
2495 error:
2496 free(active);
2497 free(reordering);
2498 isl_qpolynomial_free(qp);
2499 return NULL;
2502 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2503 unsigned first, unsigned n)
2505 int i;
2506 struct isl_upoly_rec *rec;
2508 if (!up)
2509 return NULL;
2510 if (n == 0 || up->var < 0 || up->var < first)
2511 return up;
2512 if (up->var < first + n) {
2513 up = replace_by_constant_term(up);
2514 return isl_upoly_drop(up, first, n);
2516 up = isl_upoly_cow(up);
2517 if (!up)
2518 return NULL;
2519 up->var -= n;
2520 rec = isl_upoly_as_rec(up);
2521 if (!rec)
2522 goto error;
2524 for (i = 0; i < rec->n; ++i) {
2525 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2526 if (!rec->p[i])
2527 goto error;
2530 return up;
2531 error:
2532 isl_upoly_free(up);
2533 return NULL;
2536 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2537 __isl_take isl_qpolynomial *qp,
2538 enum isl_dim_type type, unsigned pos, const char *s)
2540 qp = isl_qpolynomial_cow(qp);
2541 if (!qp)
2542 return NULL;
2543 qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2544 if (!qp->dim)
2545 goto error;
2546 return qp;
2547 error:
2548 isl_qpolynomial_free(qp);
2549 return NULL;
2552 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2553 __isl_take isl_qpolynomial *qp,
2554 enum isl_dim_type type, unsigned first, unsigned n)
2556 if (!qp)
2557 return NULL;
2558 if (type == isl_dim_out)
2559 isl_die(qp->dim->ctx, isl_error_invalid,
2560 "cannot drop output/set dimension",
2561 goto error);
2562 if (type == isl_dim_in)
2563 type = isl_dim_set;
2564 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2565 return qp;
2567 qp = isl_qpolynomial_cow(qp);
2568 if (!qp)
2569 return NULL;
2571 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2572 goto error);
2573 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2574 type == isl_dim_set, goto error);
2576 qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2577 if (!qp->dim)
2578 goto error;
2580 if (type == isl_dim_set)
2581 first += isl_space_dim(qp->dim, isl_dim_param);
2583 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2584 if (!qp->div)
2585 goto error;
2587 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2588 if (!qp->upoly)
2589 goto error;
2591 return qp;
2592 error:
2593 isl_qpolynomial_free(qp);
2594 return NULL;
2597 /* Project the domain of the quasi-polynomial onto its parameter space.
2598 * The quasi-polynomial may not involve any of the domain dimensions.
2600 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2601 __isl_take isl_qpolynomial *qp)
2603 isl_space *space;
2604 unsigned n;
2605 int involves;
2607 n = isl_qpolynomial_dim(qp, isl_dim_in);
2608 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2609 if (involves < 0)
2610 return isl_qpolynomial_free(qp);
2611 if (involves)
2612 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2613 "polynomial involves some of the domain dimensions",
2614 return isl_qpolynomial_free(qp));
2615 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2616 space = isl_qpolynomial_get_domain_space(qp);
2617 space = isl_space_params(space);
2618 qp = isl_qpolynomial_reset_domain_space(qp, space);
2619 return qp;
2622 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2623 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2625 int i, j, k;
2626 isl_int denom;
2627 unsigned total;
2628 unsigned n_div;
2629 struct isl_upoly *up;
2631 if (!eq)
2632 goto error;
2633 if (eq->n_eq == 0) {
2634 isl_basic_set_free(eq);
2635 return qp;
2638 qp = isl_qpolynomial_cow(qp);
2639 if (!qp)
2640 goto error;
2641 qp->div = isl_mat_cow(qp->div);
2642 if (!qp->div)
2643 goto error;
2645 total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2646 n_div = eq->n_div;
2647 isl_int_init(denom);
2648 for (i = 0; i < eq->n_eq; ++i) {
2649 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2650 if (j < 0 || j == 0 || j >= total)
2651 continue;
2653 for (k = 0; k < qp->div->n_row; ++k) {
2654 if (isl_int_is_zero(qp->div->row[k][1 + j]))
2655 continue;
2656 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2657 &qp->div->row[k][0]);
2658 normalize_div(qp, k);
2661 if (isl_int_is_pos(eq->eq[i][j]))
2662 isl_seq_neg(eq->eq[i], eq->eq[i], total);
2663 isl_int_abs(denom, eq->eq[i][j]);
2664 isl_int_set_si(eq->eq[i][j], 0);
2666 up = isl_upoly_from_affine(qp->dim->ctx,
2667 eq->eq[i], denom, total);
2668 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2669 isl_upoly_free(up);
2671 isl_int_clear(denom);
2673 if (!qp->upoly)
2674 goto error;
2676 isl_basic_set_free(eq);
2678 qp = substitute_non_divs(qp);
2679 qp = sort_divs(qp);
2681 return qp;
2682 error:
2683 isl_basic_set_free(eq);
2684 isl_qpolynomial_free(qp);
2685 return NULL;
2688 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2690 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2691 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2693 if (!qp || !eq)
2694 goto error;
2695 if (qp->div->n_row > 0)
2696 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2697 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2698 error:
2699 isl_basic_set_free(eq);
2700 isl_qpolynomial_free(qp);
2701 return NULL;
2704 static __isl_give isl_basic_set *add_div_constraints(
2705 __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2707 int i;
2708 unsigned total;
2710 if (!bset || !div)
2711 goto error;
2713 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2714 if (!bset)
2715 goto error;
2716 total = isl_basic_set_total_dim(bset);
2717 for (i = 0; i < div->n_row; ++i)
2718 if (isl_basic_set_add_div_constraints_var(bset,
2719 total - div->n_row + i, div->row[i]) < 0)
2720 goto error;
2722 isl_mat_free(div);
2723 return bset;
2724 error:
2725 isl_mat_free(div);
2726 isl_basic_set_free(bset);
2727 return NULL;
2730 /* Look for equalities among the variables shared by context and qp
2731 * and the integer divisions of qp, if any.
2732 * The equalities are then used to eliminate variables and/or integer
2733 * divisions from qp.
2735 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2736 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2738 isl_basic_set *aff;
2740 if (!qp)
2741 goto error;
2742 if (qp->div->n_row > 0) {
2743 isl_basic_set *bset;
2744 context = isl_set_add_dims(context, isl_dim_set,
2745 qp->div->n_row);
2746 bset = isl_basic_set_universe(isl_set_get_space(context));
2747 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2748 context = isl_set_intersect(context,
2749 isl_set_from_basic_set(bset));
2752 aff = isl_set_affine_hull(context);
2753 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2754 error:
2755 isl_qpolynomial_free(qp);
2756 isl_set_free(context);
2757 return NULL;
2760 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2761 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2763 isl_space *space = isl_qpolynomial_get_domain_space(qp);
2764 isl_set *dom_context = isl_set_universe(space);
2765 dom_context = isl_set_intersect_params(dom_context, context);
2766 return isl_qpolynomial_gist(qp, dom_context);
2769 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2770 __isl_take isl_qpolynomial *qp)
2772 isl_set *dom;
2774 if (!qp)
2775 return NULL;
2776 if (isl_qpolynomial_is_zero(qp)) {
2777 isl_space *dim = isl_qpolynomial_get_space(qp);
2778 isl_qpolynomial_free(qp);
2779 return isl_pw_qpolynomial_zero(dim);
2782 dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2783 return isl_pw_qpolynomial_alloc(dom, qp);
2786 #undef PW
2787 #define PW isl_pw_qpolynomial
2788 #undef EL
2789 #define EL isl_qpolynomial
2790 #undef EL_IS_ZERO
2791 #define EL_IS_ZERO is_zero
2792 #undef ZERO
2793 #define ZERO zero
2794 #undef IS_ZERO
2795 #define IS_ZERO is_zero
2796 #undef FIELD
2797 #define FIELD qp
2798 #undef DEFAULT_IS_ZERO
2799 #define DEFAULT_IS_ZERO 1
2801 #define NO_PULLBACK
2803 #include <isl_pw_templ.c>
2805 #undef UNION
2806 #define UNION isl_union_pw_qpolynomial
2807 #undef PART
2808 #define PART isl_pw_qpolynomial
2809 #undef PARTS
2810 #define PARTS pw_qpolynomial
2812 #include <isl_union_templ.c>
2814 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2816 if (!pwqp)
2817 return -1;
2819 if (pwqp->n != -1)
2820 return 0;
2822 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2823 return 0;
2825 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2828 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2829 __isl_take isl_pw_qpolynomial *pwqp1,
2830 __isl_take isl_pw_qpolynomial *pwqp2)
2832 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2835 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2836 __isl_take isl_pw_qpolynomial *pwqp1,
2837 __isl_take isl_pw_qpolynomial *pwqp2)
2839 int i, j, n;
2840 struct isl_pw_qpolynomial *res;
2842 if (!pwqp1 || !pwqp2)
2843 goto error;
2845 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2846 goto error);
2848 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2849 isl_pw_qpolynomial_free(pwqp2);
2850 return pwqp1;
2853 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2854 isl_pw_qpolynomial_free(pwqp1);
2855 return pwqp2;
2858 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2859 isl_pw_qpolynomial_free(pwqp1);
2860 return pwqp2;
2863 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2864 isl_pw_qpolynomial_free(pwqp2);
2865 return pwqp1;
2868 n = pwqp1->n * pwqp2->n;
2869 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2871 for (i = 0; i < pwqp1->n; ++i) {
2872 for (j = 0; j < pwqp2->n; ++j) {
2873 struct isl_set *common;
2874 struct isl_qpolynomial *prod;
2875 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2876 isl_set_copy(pwqp2->p[j].set));
2877 if (isl_set_plain_is_empty(common)) {
2878 isl_set_free(common);
2879 continue;
2882 prod = isl_qpolynomial_mul(
2883 isl_qpolynomial_copy(pwqp1->p[i].qp),
2884 isl_qpolynomial_copy(pwqp2->p[j].qp));
2886 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2890 isl_pw_qpolynomial_free(pwqp1);
2891 isl_pw_qpolynomial_free(pwqp2);
2893 return res;
2894 error:
2895 isl_pw_qpolynomial_free(pwqp1);
2896 isl_pw_qpolynomial_free(pwqp2);
2897 return NULL;
2900 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
2901 __isl_take isl_vec *vec)
2903 int i;
2904 struct isl_upoly_rec *rec;
2905 isl_val *res;
2906 isl_val *base;
2908 if (isl_upoly_is_cst(up)) {
2909 isl_vec_free(vec);
2910 res = isl_upoly_get_constant_val(up);
2911 isl_upoly_free(up);
2912 return res;
2915 rec = isl_upoly_as_rec(up);
2916 if (!rec)
2917 goto error;
2919 isl_assert(up->ctx, rec->n >= 1, goto error);
2921 base = isl_val_rat_from_isl_int(up->ctx,
2922 vec->el[1 + up->var], vec->el[0]);
2924 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2925 isl_vec_copy(vec));
2927 for (i = rec->n - 2; i >= 0; --i) {
2928 res = isl_val_mul(res, isl_val_copy(base));
2929 res = isl_val_add(res,
2930 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2931 isl_vec_copy(vec)));
2934 isl_val_free(base);
2935 isl_upoly_free(up);
2936 isl_vec_free(vec);
2937 return res;
2938 error:
2939 isl_upoly_free(up);
2940 isl_vec_free(vec);
2941 return NULL;
2944 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
2945 __isl_take isl_point *pnt)
2947 isl_vec *ext;
2948 isl_val *v;
2950 if (!qp || !pnt)
2951 goto error;
2952 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2954 if (qp->div->n_row == 0)
2955 ext = isl_vec_copy(pnt->vec);
2956 else {
2957 int i;
2958 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2959 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2960 if (!ext)
2961 goto error;
2963 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2964 for (i = 0; i < qp->div->n_row; ++i) {
2965 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2966 1 + dim + i, &ext->el[1+dim+i]);
2967 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2968 qp->div->row[i][0]);
2972 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2974 isl_qpolynomial_free(qp);
2975 isl_point_free(pnt);
2977 return v;
2978 error:
2979 isl_qpolynomial_free(qp);
2980 isl_point_free(pnt);
2981 return NULL;
2984 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2985 __isl_keep struct isl_upoly_cst *cst2)
2987 int cmp;
2988 isl_int t;
2989 isl_int_init(t);
2990 isl_int_mul(t, cst1->n, cst2->d);
2991 isl_int_submul(t, cst2->n, cst1->d);
2992 cmp = isl_int_sgn(t);
2993 isl_int_clear(t);
2994 return cmp;
2997 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2998 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2999 unsigned first, unsigned n)
3001 unsigned total;
3002 unsigned g_pos;
3003 int *exp;
3005 if (!qp)
3006 return NULL;
3007 if (type == isl_dim_out)
3008 isl_die(qp->div->ctx, isl_error_invalid,
3009 "cannot insert output/set dimensions",
3010 goto error);
3011 if (type == isl_dim_in)
3012 type = isl_dim_set;
3013 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3014 return qp;
3016 qp = isl_qpolynomial_cow(qp);
3017 if (!qp)
3018 return NULL;
3020 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
3021 goto error);
3023 g_pos = pos(qp->dim, type) + first;
3025 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3026 if (!qp->div)
3027 goto error;
3029 total = qp->div->n_col - 2;
3030 if (total > g_pos) {
3031 int i;
3032 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3033 if (!exp)
3034 goto error;
3035 for (i = 0; i < total - g_pos; ++i)
3036 exp[i] = i + n;
3037 qp->upoly = expand(qp->upoly, exp, g_pos);
3038 free(exp);
3039 if (!qp->upoly)
3040 goto error;
3043 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3044 if (!qp->dim)
3045 goto error;
3047 return qp;
3048 error:
3049 isl_qpolynomial_free(qp);
3050 return NULL;
3053 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3054 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3056 unsigned pos;
3058 pos = isl_qpolynomial_dim(qp, type);
3060 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3063 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3064 __isl_take isl_pw_qpolynomial *pwqp,
3065 enum isl_dim_type type, unsigned n)
3067 unsigned pos;
3069 pos = isl_pw_qpolynomial_dim(pwqp, type);
3071 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3074 static int *reordering_move(isl_ctx *ctx,
3075 unsigned len, unsigned dst, unsigned src, unsigned n)
3077 int i;
3078 int *reordering;
3080 reordering = isl_alloc_array(ctx, int, len);
3081 if (!reordering)
3082 return NULL;
3084 if (dst <= src) {
3085 for (i = 0; i < dst; ++i)
3086 reordering[i] = i;
3087 for (i = 0; i < n; ++i)
3088 reordering[src + i] = dst + i;
3089 for (i = 0; i < src - dst; ++i)
3090 reordering[dst + i] = dst + n + i;
3091 for (i = 0; i < len - src - n; ++i)
3092 reordering[src + n + i] = src + n + i;
3093 } else {
3094 for (i = 0; i < src; ++i)
3095 reordering[i] = i;
3096 for (i = 0; i < n; ++i)
3097 reordering[src + i] = dst + i;
3098 for (i = 0; i < dst - src; ++i)
3099 reordering[src + n + i] = src + i;
3100 for (i = 0; i < len - dst - n; ++i)
3101 reordering[dst + n + i] = dst + n + i;
3104 return reordering;
3107 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3108 __isl_take isl_qpolynomial *qp,
3109 enum isl_dim_type dst_type, unsigned dst_pos,
3110 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3112 unsigned g_dst_pos;
3113 unsigned g_src_pos;
3114 int *reordering;
3116 if (n == 0)
3117 return qp;
3119 qp = isl_qpolynomial_cow(qp);
3120 if (!qp)
3121 return NULL;
3123 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3124 isl_die(qp->dim->ctx, isl_error_invalid,
3125 "cannot move output/set dimension",
3126 goto error);
3127 if (dst_type == isl_dim_in)
3128 dst_type = isl_dim_set;
3129 if (src_type == isl_dim_in)
3130 src_type = isl_dim_set;
3132 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3133 goto error);
3135 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3136 g_src_pos = pos(qp->dim, src_type) + src_pos;
3137 if (dst_type > src_type)
3138 g_dst_pos -= n;
3140 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3141 if (!qp->div)
3142 goto error;
3143 qp = sort_divs(qp);
3144 if (!qp)
3145 goto error;
3147 reordering = reordering_move(qp->dim->ctx,
3148 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3149 if (!reordering)
3150 goto error;
3152 qp->upoly = reorder(qp->upoly, reordering);
3153 free(reordering);
3154 if (!qp->upoly)
3155 goto error;
3157 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3158 if (!qp->dim)
3159 goto error;
3161 return qp;
3162 error:
3163 isl_qpolynomial_free(qp);
3164 return NULL;
3167 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3168 isl_int *f, isl_int denom)
3170 struct isl_upoly *up;
3172 dim = isl_space_domain(dim);
3173 if (!dim)
3174 return NULL;
3176 up = isl_upoly_from_affine(dim->ctx, f, denom,
3177 1 + isl_space_dim(dim, isl_dim_all));
3179 return isl_qpolynomial_alloc(dim, 0, up);
3182 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3184 isl_ctx *ctx;
3185 struct isl_upoly *up;
3186 isl_qpolynomial *qp;
3188 if (!aff)
3189 return NULL;
3191 ctx = isl_aff_get_ctx(aff);
3192 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3193 aff->v->size - 1);
3195 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3196 aff->ls->div->n_row, up);
3197 if (!qp)
3198 goto error;
3200 isl_mat_free(qp->div);
3201 qp->div = isl_mat_copy(aff->ls->div);
3202 qp->div = isl_mat_cow(qp->div);
3203 if (!qp->div)
3204 goto error;
3206 isl_aff_free(aff);
3207 qp = reduce_divs(qp);
3208 qp = remove_redundant_divs(qp);
3209 return qp;
3210 error:
3211 isl_aff_free(aff);
3212 return isl_qpolynomial_free(qp);
3215 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3216 __isl_take isl_pw_aff *pwaff)
3218 int i;
3219 isl_pw_qpolynomial *pwqp;
3221 if (!pwaff)
3222 return NULL;
3224 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3225 pwaff->n);
3227 for (i = 0; i < pwaff->n; ++i) {
3228 isl_set *dom;
3229 isl_qpolynomial *qp;
3231 dom = isl_set_copy(pwaff->p[i].set);
3232 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3233 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3236 isl_pw_aff_free(pwaff);
3237 return pwqp;
3240 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3241 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3243 isl_aff *aff;
3245 aff = isl_constraint_get_bound(c, type, pos);
3246 isl_constraint_free(c);
3247 return isl_qpolynomial_from_aff(aff);
3250 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3251 * in "qp" by subs[i].
3253 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3254 __isl_take isl_qpolynomial *qp,
3255 enum isl_dim_type type, unsigned first, unsigned n,
3256 __isl_keep isl_qpolynomial **subs)
3258 int i;
3259 struct isl_upoly **ups;
3261 if (n == 0)
3262 return qp;
3264 qp = isl_qpolynomial_cow(qp);
3265 if (!qp)
3266 return NULL;
3268 if (type == isl_dim_out)
3269 isl_die(qp->dim->ctx, isl_error_invalid,
3270 "cannot substitute output/set dimension",
3271 goto error);
3272 if (type == isl_dim_in)
3273 type = isl_dim_set;
3275 for (i = 0; i < n; ++i)
3276 if (!subs[i])
3277 goto error;
3279 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3280 goto error);
3282 for (i = 0; i < n; ++i)
3283 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3284 goto error);
3286 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3287 for (i = 0; i < n; ++i)
3288 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3290 first += pos(qp->dim, type);
3292 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3293 if (!ups)
3294 goto error;
3295 for (i = 0; i < n; ++i)
3296 ups[i] = subs[i]->upoly;
3298 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3300 free(ups);
3302 if (!qp->upoly)
3303 goto error;
3305 return qp;
3306 error:
3307 isl_qpolynomial_free(qp);
3308 return NULL;
3311 /* Extend "bset" with extra set dimensions for each integer division
3312 * in "qp" and then call "fn" with the extended bset and the polynomial
3313 * that results from replacing each of the integer divisions by the
3314 * corresponding extra set dimension.
3316 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3317 __isl_keep isl_basic_set *bset,
3318 int (*fn)(__isl_take isl_basic_set *bset,
3319 __isl_take isl_qpolynomial *poly, void *user), void *user)
3321 isl_space *dim;
3322 isl_mat *div;
3323 isl_qpolynomial *poly;
3325 if (!qp || !bset)
3326 goto error;
3327 if (qp->div->n_row == 0)
3328 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3329 user);
3331 div = isl_mat_copy(qp->div);
3332 dim = isl_space_copy(qp->dim);
3333 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3334 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3335 bset = isl_basic_set_copy(bset);
3336 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3337 bset = add_div_constraints(bset, div);
3339 return fn(bset, poly, user);
3340 error:
3341 return -1;
3344 /* Return total degree in variables first (inclusive) up to last (exclusive).
3346 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3348 int deg = -1;
3349 int i;
3350 struct isl_upoly_rec *rec;
3352 if (!up)
3353 return -2;
3354 if (isl_upoly_is_zero(up))
3355 return -1;
3356 if (isl_upoly_is_cst(up) || up->var < first)
3357 return 0;
3359 rec = isl_upoly_as_rec(up);
3360 if (!rec)
3361 return -2;
3363 for (i = 0; i < rec->n; ++i) {
3364 int d;
3366 if (isl_upoly_is_zero(rec->p[i]))
3367 continue;
3368 d = isl_upoly_degree(rec->p[i], first, last);
3369 if (up->var < last)
3370 d += i;
3371 if (d > deg)
3372 deg = d;
3375 return deg;
3378 /* Return total degree in set variables.
3380 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3382 unsigned ovar;
3383 unsigned nvar;
3385 if (!poly)
3386 return -2;
3388 ovar = isl_space_offset(poly->dim, isl_dim_set);
3389 nvar = isl_space_dim(poly->dim, isl_dim_set);
3390 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3393 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3394 unsigned pos, int deg)
3396 int i;
3397 struct isl_upoly_rec *rec;
3399 if (!up)
3400 return NULL;
3402 if (isl_upoly_is_cst(up) || up->var < pos) {
3403 if (deg == 0)
3404 return isl_upoly_copy(up);
3405 else
3406 return isl_upoly_zero(up->ctx);
3409 rec = isl_upoly_as_rec(up);
3410 if (!rec)
3411 return NULL;
3413 if (up->var == pos) {
3414 if (deg < rec->n)
3415 return isl_upoly_copy(rec->p[deg]);
3416 else
3417 return isl_upoly_zero(up->ctx);
3420 up = isl_upoly_copy(up);
3421 up = isl_upoly_cow(up);
3422 rec = isl_upoly_as_rec(up);
3423 if (!rec)
3424 goto error;
3426 for (i = 0; i < rec->n; ++i) {
3427 struct isl_upoly *t;
3428 t = isl_upoly_coeff(rec->p[i], pos, deg);
3429 if (!t)
3430 goto error;
3431 isl_upoly_free(rec->p[i]);
3432 rec->p[i] = t;
3435 return up;
3436 error:
3437 isl_upoly_free(up);
3438 return NULL;
3441 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3443 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3444 __isl_keep isl_qpolynomial *qp,
3445 enum isl_dim_type type, unsigned t_pos, int deg)
3447 unsigned g_pos;
3448 struct isl_upoly *up;
3449 isl_qpolynomial *c;
3451 if (!qp)
3452 return NULL;
3454 if (type == isl_dim_out)
3455 isl_die(qp->div->ctx, isl_error_invalid,
3456 "output/set dimension does not have a coefficient",
3457 return NULL);
3458 if (type == isl_dim_in)
3459 type = isl_dim_set;
3461 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3462 return NULL);
3464 g_pos = pos(qp->dim, type) + t_pos;
3465 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3467 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3468 if (!c)
3469 return NULL;
3470 isl_mat_free(c->div);
3471 c->div = isl_mat_copy(qp->div);
3472 if (!c->div)
3473 goto error;
3474 return c;
3475 error:
3476 isl_qpolynomial_free(c);
3477 return NULL;
3480 /* Homogenize the polynomial in the variables first (inclusive) up to
3481 * last (exclusive) by inserting powers of variable first.
3482 * Variable first is assumed not to appear in the input.
3484 __isl_give struct isl_upoly *isl_upoly_homogenize(
3485 __isl_take struct isl_upoly *up, int deg, int target,
3486 int first, int last)
3488 int i;
3489 struct isl_upoly_rec *rec;
3491 if (!up)
3492 return NULL;
3493 if (isl_upoly_is_zero(up))
3494 return up;
3495 if (deg == target)
3496 return up;
3497 if (isl_upoly_is_cst(up) || up->var < first) {
3498 struct isl_upoly *hom;
3500 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3501 if (!hom)
3502 goto error;
3503 rec = isl_upoly_as_rec(hom);
3504 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3506 return hom;
3509 up = isl_upoly_cow(up);
3510 rec = isl_upoly_as_rec(up);
3511 if (!rec)
3512 goto error;
3514 for (i = 0; i < rec->n; ++i) {
3515 if (isl_upoly_is_zero(rec->p[i]))
3516 continue;
3517 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3518 up->var < last ? deg + i : i, target,
3519 first, last);
3520 if (!rec->p[i])
3521 goto error;
3524 return up;
3525 error:
3526 isl_upoly_free(up);
3527 return NULL;
3530 /* Homogenize the polynomial in the set variables by introducing
3531 * powers of an extra set variable at position 0.
3533 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3534 __isl_take isl_qpolynomial *poly)
3536 unsigned ovar;
3537 unsigned nvar;
3538 int deg = isl_qpolynomial_degree(poly);
3540 if (deg < -1)
3541 goto error;
3543 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3544 poly = isl_qpolynomial_cow(poly);
3545 if (!poly)
3546 goto error;
3548 ovar = isl_space_offset(poly->dim, isl_dim_set);
3549 nvar = isl_space_dim(poly->dim, isl_dim_set);
3550 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3551 ovar, ovar + nvar);
3552 if (!poly->upoly)
3553 goto error;
3555 return poly;
3556 error:
3557 isl_qpolynomial_free(poly);
3558 return NULL;
3561 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3562 __isl_take isl_mat *div)
3564 isl_term *term;
3565 int n;
3567 if (!dim || !div)
3568 goto error;
3570 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3572 term = isl_calloc(dim->ctx, struct isl_term,
3573 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3574 if (!term)
3575 goto error;
3577 term->ref = 1;
3578 term->dim = dim;
3579 term->div = div;
3580 isl_int_init(term->n);
3581 isl_int_init(term->d);
3583 return term;
3584 error:
3585 isl_space_free(dim);
3586 isl_mat_free(div);
3587 return NULL;
3590 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3592 if (!term)
3593 return NULL;
3595 term->ref++;
3596 return term;
3599 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3601 int i;
3602 isl_term *dup;
3603 unsigned total;
3605 if (!term)
3606 return NULL;
3608 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3610 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3611 if (!dup)
3612 return NULL;
3614 isl_int_set(dup->n, term->n);
3615 isl_int_set(dup->d, term->d);
3617 for (i = 0; i < total; ++i)
3618 dup->pow[i] = term->pow[i];
3620 return dup;
3623 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3625 if (!term)
3626 return NULL;
3628 if (term->ref == 1)
3629 return term;
3630 term->ref--;
3631 return isl_term_dup(term);
3634 void isl_term_free(__isl_take isl_term *term)
3636 if (!term)
3637 return;
3639 if (--term->ref > 0)
3640 return;
3642 isl_space_free(term->dim);
3643 isl_mat_free(term->div);
3644 isl_int_clear(term->n);
3645 isl_int_clear(term->d);
3646 free(term);
3649 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3651 if (!term)
3652 return 0;
3654 switch (type) {
3655 case isl_dim_param:
3656 case isl_dim_in:
3657 case isl_dim_out: return isl_space_dim(term->dim, type);
3658 case isl_dim_div: return term->div->n_row;
3659 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3660 term->div->n_row;
3661 default: return 0;
3665 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3667 return term ? term->dim->ctx : NULL;
3670 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3672 if (!term)
3673 return;
3674 isl_int_set(*n, term->n);
3677 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3679 if (!term)
3680 return;
3681 isl_int_set(*d, term->d);
3684 /* Return the coefficient of the term "term".
3686 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3688 if (!term)
3689 return NULL;
3691 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3692 term->n, term->d);
3695 int isl_term_get_exp(__isl_keep isl_term *term,
3696 enum isl_dim_type type, unsigned pos)
3698 if (!term)
3699 return -1;
3701 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3703 if (type >= isl_dim_set)
3704 pos += isl_space_dim(term->dim, isl_dim_param);
3705 if (type >= isl_dim_div)
3706 pos += isl_space_dim(term->dim, isl_dim_set);
3708 return term->pow[pos];
3711 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3713 isl_local_space *ls;
3714 isl_aff *aff;
3716 if (!term)
3717 return NULL;
3719 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3720 return NULL);
3722 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3723 isl_mat_copy(term->div));
3724 aff = isl_aff_alloc(ls);
3725 if (!aff)
3726 return NULL;
3728 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3730 aff = isl_aff_normalize(aff);
3732 return aff;
3735 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3736 int (*fn)(__isl_take isl_term *term, void *user),
3737 __isl_take isl_term *term, void *user)
3739 int i;
3740 struct isl_upoly_rec *rec;
3742 if (!up || !term)
3743 goto error;
3745 if (isl_upoly_is_zero(up))
3746 return term;
3748 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3749 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3750 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3752 if (isl_upoly_is_cst(up)) {
3753 struct isl_upoly_cst *cst;
3754 cst = isl_upoly_as_cst(up);
3755 if (!cst)
3756 goto error;
3757 term = isl_term_cow(term);
3758 if (!term)
3759 goto error;
3760 isl_int_set(term->n, cst->n);
3761 isl_int_set(term->d, cst->d);
3762 if (fn(isl_term_copy(term), user) < 0)
3763 goto error;
3764 return term;
3767 rec = isl_upoly_as_rec(up);
3768 if (!rec)
3769 goto error;
3771 for (i = 0; i < rec->n; ++i) {
3772 term = isl_term_cow(term);
3773 if (!term)
3774 goto error;
3775 term->pow[up->var] = i;
3776 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3777 if (!term)
3778 goto error;
3780 term->pow[up->var] = 0;
3782 return term;
3783 error:
3784 isl_term_free(term);
3785 return NULL;
3788 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3789 int (*fn)(__isl_take isl_term *term, void *user), void *user)
3791 isl_term *term;
3793 if (!qp)
3794 return -1;
3796 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3797 if (!term)
3798 return -1;
3800 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3802 isl_term_free(term);
3804 return term ? 0 : -1;
3807 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3809 struct isl_upoly *up;
3810 isl_qpolynomial *qp;
3811 int i, n;
3813 if (!term)
3814 return NULL;
3816 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3818 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3819 for (i = 0; i < n; ++i) {
3820 if (!term->pow[i])
3821 continue;
3822 up = isl_upoly_mul(up,
3823 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3826 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3827 if (!qp)
3828 goto error;
3829 isl_mat_free(qp->div);
3830 qp->div = isl_mat_copy(term->div);
3831 if (!qp->div)
3832 goto error;
3834 isl_term_free(term);
3835 return qp;
3836 error:
3837 isl_qpolynomial_free(qp);
3838 isl_term_free(term);
3839 return NULL;
3842 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3843 __isl_take isl_space *dim)
3845 int i;
3846 int extra;
3847 unsigned total;
3849 if (!qp || !dim)
3850 goto error;
3852 if (isl_space_is_equal(qp->dim, dim)) {
3853 isl_space_free(dim);
3854 return qp;
3857 qp = isl_qpolynomial_cow(qp);
3858 if (!qp)
3859 goto error;
3861 extra = isl_space_dim(dim, isl_dim_set) -
3862 isl_space_dim(qp->dim, isl_dim_set);
3863 total = isl_space_dim(qp->dim, isl_dim_all);
3864 if (qp->div->n_row) {
3865 int *exp;
3867 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3868 if (!exp)
3869 goto error;
3870 for (i = 0; i < qp->div->n_row; ++i)
3871 exp[i] = extra + i;
3872 qp->upoly = expand(qp->upoly, exp, total);
3873 free(exp);
3874 if (!qp->upoly)
3875 goto error;
3877 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3878 if (!qp->div)
3879 goto error;
3880 for (i = 0; i < qp->div->n_row; ++i)
3881 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3883 isl_space_free(qp->dim);
3884 qp->dim = dim;
3886 return qp;
3887 error:
3888 isl_space_free(dim);
3889 isl_qpolynomial_free(qp);
3890 return NULL;
3893 /* For each parameter or variable that does not appear in qp,
3894 * first eliminate the variable from all constraints and then set it to zero.
3896 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3897 __isl_keep isl_qpolynomial *qp)
3899 int *active = NULL;
3900 int i;
3901 int d;
3902 unsigned nparam;
3903 unsigned nvar;
3905 if (!set || !qp)
3906 goto error;
3908 d = isl_space_dim(set->dim, isl_dim_all);
3909 active = isl_calloc_array(set->ctx, int, d);
3910 if (set_active(qp, active) < 0)
3911 goto error;
3913 for (i = 0; i < d; ++i)
3914 if (!active[i])
3915 break;
3917 if (i == d) {
3918 free(active);
3919 return set;
3922 nparam = isl_space_dim(set->dim, isl_dim_param);
3923 nvar = isl_space_dim(set->dim, isl_dim_set);
3924 for (i = 0; i < nparam; ++i) {
3925 if (active[i])
3926 continue;
3927 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3928 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3930 for (i = 0; i < nvar; ++i) {
3931 if (active[nparam + i])
3932 continue;
3933 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3934 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3937 free(active);
3939 return set;
3940 error:
3941 free(active);
3942 isl_set_free(set);
3943 return NULL;
3946 struct isl_opt_data {
3947 isl_qpolynomial *qp;
3948 int first;
3949 isl_val *opt;
3950 int max;
3953 static int opt_fn(__isl_take isl_point *pnt, void *user)
3955 struct isl_opt_data *data = (struct isl_opt_data *)user;
3956 isl_val *val;
3958 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3959 if (data->first) {
3960 data->first = 0;
3961 data->opt = val;
3962 } else if (data->max) {
3963 data->opt = isl_val_max(data->opt, val);
3964 } else {
3965 data->opt = isl_val_min(data->opt, val);
3968 return 0;
3971 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
3972 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3974 struct isl_opt_data data = { NULL, 1, NULL, max };
3976 if (!set || !qp)
3977 goto error;
3979 if (isl_upoly_is_cst(qp->upoly)) {
3980 isl_set_free(set);
3981 data.opt = isl_qpolynomial_get_constant_val(qp);
3982 isl_qpolynomial_free(qp);
3983 return data.opt;
3986 set = fix_inactive(set, qp);
3988 data.qp = qp;
3989 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3990 goto error;
3992 if (data.first)
3993 data.opt = isl_val_zero(isl_set_get_ctx(set));
3995 isl_set_free(set);
3996 isl_qpolynomial_free(qp);
3997 return data.opt;
3998 error:
3999 isl_set_free(set);
4000 isl_qpolynomial_free(qp);
4001 isl_val_free(data.opt);
4002 return NULL;
4005 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
4006 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
4008 int i;
4009 int n_sub;
4010 isl_ctx *ctx;
4011 struct isl_upoly **subs;
4012 isl_mat *mat, *diag;
4014 qp = isl_qpolynomial_cow(qp);
4015 if (!qp || !morph)
4016 goto error;
4018 ctx = qp->dim->ctx;
4019 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
4021 n_sub = morph->inv->n_row - 1;
4022 if (morph->inv->n_row != morph->inv->n_col)
4023 n_sub += qp->div->n_row;
4024 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
4025 if (n_sub && !subs)
4026 goto error;
4028 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4029 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4030 morph->inv->row[0][0], morph->inv->n_col);
4031 if (morph->inv->n_row != morph->inv->n_col)
4032 for (i = 0; i < qp->div->n_row; ++i)
4033 subs[morph->inv->n_row - 1 + i] =
4034 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4036 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4038 for (i = 0; i < n_sub; ++i)
4039 isl_upoly_free(subs[i]);
4040 free(subs);
4042 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4043 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4044 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4045 mat = isl_mat_diagonal(mat, diag);
4046 qp->div = isl_mat_product(qp->div, mat);
4047 isl_space_free(qp->dim);
4048 qp->dim = isl_space_copy(morph->ran->dim);
4050 if (!qp->upoly || !qp->div || !qp->dim)
4051 goto error;
4053 isl_morph_free(morph);
4055 return qp;
4056 error:
4057 isl_qpolynomial_free(qp);
4058 isl_morph_free(morph);
4059 return NULL;
4062 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4063 __isl_take isl_union_pw_qpolynomial *upwqp1,
4064 __isl_take isl_union_pw_qpolynomial *upwqp2)
4066 return isl_union_pw_qpolynomial_match_bin_op(upwqp1, upwqp2,
4067 &isl_pw_qpolynomial_mul);
4070 /* Reorder the columns of the given div definitions according to the
4071 * given reordering.
4073 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4074 __isl_take isl_reordering *r)
4076 int i, j;
4077 isl_mat *mat;
4078 int extra;
4080 if (!div || !r)
4081 goto error;
4083 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4084 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4085 if (!mat)
4086 goto error;
4088 for (i = 0; i < div->n_row; ++i) {
4089 isl_seq_cpy(mat->row[i], div->row[i], 2);
4090 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4091 for (j = 0; j < r->len; ++j)
4092 isl_int_set(mat->row[i][2 + r->pos[j]],
4093 div->row[i][2 + j]);
4096 isl_reordering_free(r);
4097 isl_mat_free(div);
4098 return mat;
4099 error:
4100 isl_reordering_free(r);
4101 isl_mat_free(div);
4102 return NULL;
4105 /* Reorder the dimension of "qp" according to the given reordering.
4107 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4108 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4110 qp = isl_qpolynomial_cow(qp);
4111 if (!qp)
4112 goto error;
4114 r = isl_reordering_extend(r, qp->div->n_row);
4115 if (!r)
4116 goto error;
4118 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4119 if (!qp->div)
4120 goto error;
4122 qp->upoly = reorder(qp->upoly, r->pos);
4123 if (!qp->upoly)
4124 goto error;
4126 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4128 isl_reordering_free(r);
4129 return qp;
4130 error:
4131 isl_qpolynomial_free(qp);
4132 isl_reordering_free(r);
4133 return NULL;
4136 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4137 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4139 if (!qp || !model)
4140 goto error;
4142 if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4143 isl_reordering *exp;
4145 model = isl_space_drop_dims(model, isl_dim_in,
4146 0, isl_space_dim(model, isl_dim_in));
4147 model = isl_space_drop_dims(model, isl_dim_out,
4148 0, isl_space_dim(model, isl_dim_out));
4149 exp = isl_parameter_alignment_reordering(qp->dim, model);
4150 exp = isl_reordering_extend_space(exp,
4151 isl_qpolynomial_get_domain_space(qp));
4152 qp = isl_qpolynomial_realign_domain(qp, exp);
4155 isl_space_free(model);
4156 return qp;
4157 error:
4158 isl_space_free(model);
4159 isl_qpolynomial_free(qp);
4160 return NULL;
4163 struct isl_split_periods_data {
4164 int max_periods;
4165 isl_pw_qpolynomial *res;
4168 /* Create a slice where the integer division "div" has the fixed value "v".
4169 * In particular, if "div" refers to floor(f/m), then create a slice
4171 * m v <= f <= m v + (m - 1)
4173 * or
4175 * f - m v >= 0
4176 * -f + m v + (m - 1) >= 0
4178 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4179 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4181 int total;
4182 isl_basic_set *bset = NULL;
4183 int k;
4185 if (!dim || !qp)
4186 goto error;
4188 total = isl_space_dim(dim, isl_dim_all);
4189 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4191 k = isl_basic_set_alloc_inequality(bset);
4192 if (k < 0)
4193 goto error;
4194 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4195 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4197 k = isl_basic_set_alloc_inequality(bset);
4198 if (k < 0)
4199 goto error;
4200 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4201 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4202 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4203 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4205 isl_space_free(dim);
4206 return isl_set_from_basic_set(bset);
4207 error:
4208 isl_basic_set_free(bset);
4209 isl_space_free(dim);
4210 return NULL;
4213 static int split_periods(__isl_take isl_set *set,
4214 __isl_take isl_qpolynomial *qp, void *user);
4216 /* Create a slice of the domain "set" such that integer division "div"
4217 * has the fixed value "v" and add the results to data->res,
4218 * replacing the integer division by "v" in "qp".
4220 static int set_div(__isl_take isl_set *set,
4221 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4222 struct isl_split_periods_data *data)
4224 int i;
4225 int total;
4226 isl_set *slice;
4227 struct isl_upoly *cst;
4229 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4230 set = isl_set_intersect(set, slice);
4232 if (!qp)
4233 goto error;
4235 total = isl_space_dim(qp->dim, isl_dim_all);
4237 for (i = div + 1; i < qp->div->n_row; ++i) {
4238 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4239 continue;
4240 isl_int_addmul(qp->div->row[i][1],
4241 qp->div->row[i][2 + total + div], v);
4242 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4245 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4246 qp = substitute_div(qp, div, cst);
4248 return split_periods(set, qp, data);
4249 error:
4250 isl_set_free(set);
4251 isl_qpolynomial_free(qp);
4252 return -1;
4255 /* Split the domain "set" such that integer division "div"
4256 * has a fixed value (ranging from "min" to "max") on each slice
4257 * and add the results to data->res.
4259 static int split_div(__isl_take isl_set *set,
4260 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4261 struct isl_split_periods_data *data)
4263 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4264 isl_set *set_i = isl_set_copy(set);
4265 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4267 if (set_div(set_i, qp_i, div, min, data) < 0)
4268 goto error;
4270 isl_set_free(set);
4271 isl_qpolynomial_free(qp);
4272 return 0;
4273 error:
4274 isl_set_free(set);
4275 isl_qpolynomial_free(qp);
4276 return -1;
4279 /* If "qp" refers to any integer division
4280 * that can only attain "max_periods" distinct values on "set"
4281 * then split the domain along those distinct values.
4282 * Add the results (or the original if no splitting occurs)
4283 * to data->res.
4285 static int split_periods(__isl_take isl_set *set,
4286 __isl_take isl_qpolynomial *qp, void *user)
4288 int i;
4289 isl_pw_qpolynomial *pwqp;
4290 struct isl_split_periods_data *data;
4291 isl_int min, max;
4292 int total;
4293 int r = 0;
4295 data = (struct isl_split_periods_data *)user;
4297 if (!set || !qp)
4298 goto error;
4300 if (qp->div->n_row == 0) {
4301 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4302 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4303 return 0;
4306 isl_int_init(min);
4307 isl_int_init(max);
4308 total = isl_space_dim(qp->dim, isl_dim_all);
4309 for (i = 0; i < qp->div->n_row; ++i) {
4310 enum isl_lp_result lp_res;
4312 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4313 qp->div->n_row) != -1)
4314 continue;
4316 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4317 set->ctx->one, &min, NULL, NULL);
4318 if (lp_res == isl_lp_error)
4319 goto error2;
4320 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4321 continue;
4322 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4324 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4325 set->ctx->one, &max, NULL, NULL);
4326 if (lp_res == isl_lp_error)
4327 goto error2;
4328 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4329 continue;
4330 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4332 isl_int_sub(max, max, min);
4333 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4334 isl_int_add(max, max, min);
4335 break;
4339 if (i < qp->div->n_row) {
4340 r = split_div(set, qp, i, min, max, data);
4341 } else {
4342 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4343 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4346 isl_int_clear(max);
4347 isl_int_clear(min);
4349 return r;
4350 error2:
4351 isl_int_clear(max);
4352 isl_int_clear(min);
4353 error:
4354 isl_set_free(set);
4355 isl_qpolynomial_free(qp);
4356 return -1;
4359 /* If any quasi-polynomial in pwqp refers to any integer division
4360 * that can only attain "max_periods" distinct values on its domain
4361 * then split the domain along those distinct values.
4363 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4364 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4366 struct isl_split_periods_data data;
4368 data.max_periods = max_periods;
4369 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4371 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4372 goto error;
4374 isl_pw_qpolynomial_free(pwqp);
4376 return data.res;
4377 error:
4378 isl_pw_qpolynomial_free(data.res);
4379 isl_pw_qpolynomial_free(pwqp);
4380 return NULL;
4383 /* Construct a piecewise quasipolynomial that is constant on the given
4384 * domain. In particular, it is
4385 * 0 if cst == 0
4386 * 1 if cst == 1
4387 * infinity if cst == -1
4389 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4390 __isl_take isl_basic_set *bset, int cst)
4392 isl_space *dim;
4393 isl_qpolynomial *qp;
4395 if (!bset)
4396 return NULL;
4398 bset = isl_basic_set_params(bset);
4399 dim = isl_basic_set_get_space(bset);
4400 if (cst < 0)
4401 qp = isl_qpolynomial_infty_on_domain(dim);
4402 else if (cst == 0)
4403 qp = isl_qpolynomial_zero_on_domain(dim);
4404 else
4405 qp = isl_qpolynomial_one_on_domain(dim);
4406 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4409 /* Factor bset, call fn on each of the factors and return the product.
4411 * If no factors can be found, simply call fn on the input.
4412 * Otherwise, construct the factors based on the factorizer,
4413 * call fn on each factor and compute the product.
4415 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4416 __isl_take isl_basic_set *bset,
4417 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4419 int i, n;
4420 isl_space *dim;
4421 isl_set *set;
4422 isl_factorizer *f;
4423 isl_qpolynomial *qp;
4424 isl_pw_qpolynomial *pwqp;
4425 unsigned nparam;
4426 unsigned nvar;
4428 f = isl_basic_set_factorizer(bset);
4429 if (!f)
4430 goto error;
4431 if (f->n_group == 0) {
4432 isl_factorizer_free(f);
4433 return fn(bset);
4436 nparam = isl_basic_set_dim(bset, isl_dim_param);
4437 nvar = isl_basic_set_dim(bset, isl_dim_set);
4439 dim = isl_basic_set_get_space(bset);
4440 dim = isl_space_domain(dim);
4441 set = isl_set_universe(isl_space_copy(dim));
4442 qp = isl_qpolynomial_one_on_domain(dim);
4443 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4445 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4447 for (i = 0, n = 0; i < f->n_group; ++i) {
4448 isl_basic_set *bset_i;
4449 isl_pw_qpolynomial *pwqp_i;
4451 bset_i = isl_basic_set_copy(bset);
4452 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4453 nparam + n + f->len[i], nvar - n - f->len[i]);
4454 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4455 nparam, n);
4456 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4457 n + f->len[i], nvar - n - f->len[i]);
4458 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4460 pwqp_i = fn(bset_i);
4461 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4463 n += f->len[i];
4466 isl_basic_set_free(bset);
4467 isl_factorizer_free(f);
4469 return pwqp;
4470 error:
4471 isl_basic_set_free(bset);
4472 return NULL;
4475 /* Factor bset, call fn on each of the factors and return the product.
4476 * The function is assumed to evaluate to zero on empty domains,
4477 * to one on zero-dimensional domains and to infinity on unbounded domains
4478 * and will not be called explicitly on zero-dimensional or unbounded domains.
4480 * We first check for some special cases and remove all equalities.
4481 * Then we hand over control to compressed_multiplicative_call.
4483 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4484 __isl_take isl_basic_set *bset,
4485 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4487 int bounded;
4488 isl_morph *morph;
4489 isl_pw_qpolynomial *pwqp;
4491 if (!bset)
4492 return NULL;
4494 if (isl_basic_set_plain_is_empty(bset))
4495 return constant_on_domain(bset, 0);
4497 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4498 return constant_on_domain(bset, 1);
4500 bounded = isl_basic_set_is_bounded(bset);
4501 if (bounded < 0)
4502 goto error;
4503 if (!bounded)
4504 return constant_on_domain(bset, -1);
4506 if (bset->n_eq == 0)
4507 return compressed_multiplicative_call(bset, fn);
4509 morph = isl_basic_set_full_compression(bset);
4510 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4512 pwqp = compressed_multiplicative_call(bset, fn);
4514 morph = isl_morph_dom_params(morph);
4515 morph = isl_morph_ran_params(morph);
4516 morph = isl_morph_inverse(morph);
4518 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4520 return pwqp;
4521 error:
4522 isl_basic_set_free(bset);
4523 return NULL;
4526 /* Drop all floors in "qp", turning each integer division [a/m] into
4527 * a rational division a/m. If "down" is set, then the integer division
4528 * is replaced by (a-(m-1))/m instead.
4530 static __isl_give isl_qpolynomial *qp_drop_floors(
4531 __isl_take isl_qpolynomial *qp, int down)
4533 int i;
4534 struct isl_upoly *s;
4536 if (!qp)
4537 return NULL;
4538 if (qp->div->n_row == 0)
4539 return qp;
4541 qp = isl_qpolynomial_cow(qp);
4542 if (!qp)
4543 return NULL;
4545 for (i = qp->div->n_row - 1; i >= 0; --i) {
4546 if (down) {
4547 isl_int_sub(qp->div->row[i][1],
4548 qp->div->row[i][1], qp->div->row[i][0]);
4549 isl_int_add_ui(qp->div->row[i][1],
4550 qp->div->row[i][1], 1);
4552 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4553 qp->div->row[i][0], qp->div->n_col - 1);
4554 qp = substitute_div(qp, i, s);
4555 if (!qp)
4556 return NULL;
4559 return qp;
4562 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4563 * a rational division a/m.
4565 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4566 __isl_take isl_pw_qpolynomial *pwqp)
4568 int i;
4570 if (!pwqp)
4571 return NULL;
4573 if (isl_pw_qpolynomial_is_zero(pwqp))
4574 return pwqp;
4576 pwqp = isl_pw_qpolynomial_cow(pwqp);
4577 if (!pwqp)
4578 return NULL;
4580 for (i = 0; i < pwqp->n; ++i) {
4581 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4582 if (!pwqp->p[i].qp)
4583 goto error;
4586 return pwqp;
4587 error:
4588 isl_pw_qpolynomial_free(pwqp);
4589 return NULL;
4592 /* Adjust all the integer divisions in "qp" such that they are at least
4593 * one over the given orthant (identified by "signs"). This ensures
4594 * that they will still be non-negative even after subtracting (m-1)/m.
4596 * In particular, f is replaced by f' + v, changing f = [a/m]
4597 * to f' = [(a - m v)/m].
4598 * If the constant term k in a is smaller than m,
4599 * the constant term of v is set to floor(k/m) - 1.
4600 * For any other term, if the coefficient c and the variable x have
4601 * the same sign, then no changes are needed.
4602 * Otherwise, if the variable is positive (and c is negative),
4603 * then the coefficient of x in v is set to floor(c/m).
4604 * If the variable is negative (and c is positive),
4605 * then the coefficient of x in v is set to ceil(c/m).
4607 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4608 int *signs)
4610 int i, j;
4611 int total;
4612 isl_vec *v = NULL;
4613 struct isl_upoly *s;
4615 qp = isl_qpolynomial_cow(qp);
4616 if (!qp)
4617 return NULL;
4618 qp->div = isl_mat_cow(qp->div);
4619 if (!qp->div)
4620 goto error;
4622 total = isl_space_dim(qp->dim, isl_dim_all);
4623 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4625 for (i = 0; i < qp->div->n_row; ++i) {
4626 isl_int *row = qp->div->row[i];
4627 v = isl_vec_clr(v);
4628 if (!v)
4629 goto error;
4630 if (isl_int_lt(row[1], row[0])) {
4631 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4632 isl_int_sub_ui(v->el[0], v->el[0], 1);
4633 isl_int_submul(row[1], row[0], v->el[0]);
4635 for (j = 0; j < total; ++j) {
4636 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4637 continue;
4638 if (signs[j] < 0)
4639 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4640 else
4641 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4642 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4644 for (j = 0; j < i; ++j) {
4645 if (isl_int_sgn(row[2 + total + j]) >= 0)
4646 continue;
4647 isl_int_fdiv_q(v->el[1 + total + j],
4648 row[2 + total + j], row[0]);
4649 isl_int_submul(row[2 + total + j],
4650 row[0], v->el[1 + total + j]);
4652 for (j = i + 1; j < qp->div->n_row; ++j) {
4653 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4654 continue;
4655 isl_seq_combine(qp->div->row[j] + 1,
4656 qp->div->ctx->one, qp->div->row[j] + 1,
4657 qp->div->row[j][2 + total + i], v->el, v->size);
4659 isl_int_set_si(v->el[1 + total + i], 1);
4660 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4661 qp->div->ctx->one, v->size);
4662 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4663 isl_upoly_free(s);
4664 if (!qp->upoly)
4665 goto error;
4668 isl_vec_free(v);
4669 return qp;
4670 error:
4671 isl_vec_free(v);
4672 isl_qpolynomial_free(qp);
4673 return NULL;
4676 struct isl_to_poly_data {
4677 int sign;
4678 isl_pw_qpolynomial *res;
4679 isl_qpolynomial *qp;
4682 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4683 * We first make all integer divisions positive and then split the
4684 * quasipolynomials into terms with sign data->sign (the direction
4685 * of the requested approximation) and terms with the opposite sign.
4686 * In the first set of terms, each integer division [a/m] is
4687 * overapproximated by a/m, while in the second it is underapproximated
4688 * by (a-(m-1))/m.
4690 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4691 void *user)
4693 struct isl_to_poly_data *data = user;
4694 isl_pw_qpolynomial *t;
4695 isl_qpolynomial *qp, *up, *down;
4697 qp = isl_qpolynomial_copy(data->qp);
4698 qp = make_divs_pos(qp, signs);
4700 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4701 up = qp_drop_floors(up, 0);
4702 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4703 down = qp_drop_floors(down, 1);
4705 isl_qpolynomial_free(qp);
4706 qp = isl_qpolynomial_add(up, down);
4708 t = isl_pw_qpolynomial_alloc(orthant, qp);
4709 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4711 return 0;
4714 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4715 * the polynomial will be an overapproximation. If "sign" is negative,
4716 * it will be an underapproximation. If "sign" is zero, the approximation
4717 * will lie somewhere in between.
4719 * In particular, is sign == 0, we simply drop the floors, turning
4720 * the integer divisions into rational divisions.
4721 * Otherwise, we split the domains into orthants, make all integer divisions
4722 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4723 * depending on the requested sign and the sign of the term in which
4724 * the integer division appears.
4726 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4727 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4729 int i;
4730 struct isl_to_poly_data data;
4732 if (sign == 0)
4733 return pwqp_drop_floors(pwqp);
4735 if (!pwqp)
4736 return NULL;
4738 data.sign = sign;
4739 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4741 for (i = 0; i < pwqp->n; ++i) {
4742 if (pwqp->p[i].qp->div->n_row == 0) {
4743 isl_pw_qpolynomial *t;
4744 t = isl_pw_qpolynomial_alloc(
4745 isl_set_copy(pwqp->p[i].set),
4746 isl_qpolynomial_copy(pwqp->p[i].qp));
4747 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4748 continue;
4750 data.qp = pwqp->p[i].qp;
4751 if (isl_set_foreach_orthant(pwqp->p[i].set,
4752 &to_polynomial_on_orthant, &data) < 0)
4753 goto error;
4756 isl_pw_qpolynomial_free(pwqp);
4758 return data.res;
4759 error:
4760 isl_pw_qpolynomial_free(pwqp);
4761 isl_pw_qpolynomial_free(data.res);
4762 return NULL;
4765 static int poly_entry(void **entry, void *user)
4767 int *sign = user;
4768 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4770 *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4772 return *pwqp ? 0 : -1;
4775 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4776 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4778 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4779 if (!upwqp)
4780 return NULL;
4782 if (isl_hash_table_foreach(upwqp->space->ctx, &upwqp->table,
4783 &poly_entry, &sign) < 0)
4784 goto error;
4786 return upwqp;
4787 error:
4788 isl_union_pw_qpolynomial_free(upwqp);
4789 return NULL;
4792 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4793 __isl_take isl_qpolynomial *qp)
4795 int i, k;
4796 isl_space *dim;
4797 isl_vec *aff = NULL;
4798 isl_basic_map *bmap = NULL;
4799 unsigned pos;
4800 unsigned n_div;
4802 if (!qp)
4803 return NULL;
4804 if (!isl_upoly_is_affine(qp->upoly))
4805 isl_die(qp->dim->ctx, isl_error_invalid,
4806 "input quasi-polynomial not affine", goto error);
4807 aff = isl_qpolynomial_extract_affine(qp);
4808 if (!aff)
4809 goto error;
4810 dim = isl_qpolynomial_get_space(qp);
4811 pos = 1 + isl_space_offset(dim, isl_dim_out);
4812 n_div = qp->div->n_row;
4813 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4815 for (i = 0; i < n_div; ++i) {
4816 k = isl_basic_map_alloc_div(bmap);
4817 if (k < 0)
4818 goto error;
4819 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4820 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4821 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4822 goto error;
4824 k = isl_basic_map_alloc_equality(bmap);
4825 if (k < 0)
4826 goto error;
4827 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4828 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4829 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4831 isl_vec_free(aff);
4832 isl_qpolynomial_free(qp);
4833 bmap = isl_basic_map_finalize(bmap);
4834 return bmap;
4835 error:
4836 isl_vec_free(aff);
4837 isl_qpolynomial_free(qp);
4838 isl_basic_map_free(bmap);
4839 return NULL;