2 #include <isl_polynomial_private.h>
3 #include <isl_point_private.h>
4 #include <isl_dim_private.h>
5 #include <isl_map_private.h>
7 int isl_upoly_is_cst(__isl_keep
struct isl_upoly
*up
)
15 __isl_keep
struct isl_upoly_cst
*isl_upoly_as_cst(__isl_keep
struct isl_upoly
*up
)
20 isl_assert(up
->ctx
, up
->var
< 0, return NULL
);
22 return (struct isl_upoly_cst
*)up
;
25 __isl_keep
struct isl_upoly_rec
*isl_upoly_as_rec(__isl_keep
struct isl_upoly
*up
)
30 isl_assert(up
->ctx
, up
->var
>= 0, return NULL
);
32 return (struct isl_upoly_rec
*)up
;
35 int isl_upoly_is_equal(__isl_keep
struct isl_upoly
*up1
,
36 __isl_keep
struct isl_upoly
*up2
)
39 struct isl_upoly_rec
*rec1
, *rec2
;
45 if (up1
->var
!= up2
->var
)
47 if (isl_upoly_is_cst(up1
)) {
48 struct isl_upoly_cst
*cst1
, *cst2
;
49 cst1
= isl_upoly_as_cst(up1
);
50 cst2
= isl_upoly_as_cst(up2
);
53 return isl_int_eq(cst1
->n
, cst2
->n
) &&
54 isl_int_eq(cst1
->d
, cst2
->d
);
57 rec1
= isl_upoly_as_rec(up1
);
58 rec2
= isl_upoly_as_rec(up2
);
62 if (rec1
->n
!= rec2
->n
)
65 for (i
= 0; i
< rec1
->n
; ++i
) {
66 int eq
= isl_upoly_is_equal(rec1
->p
[i
], rec2
->p
[i
]);
74 int isl_upoly_is_zero(__isl_keep
struct isl_upoly
*up
)
76 struct isl_upoly_cst
*cst
;
80 if (!isl_upoly_is_cst(up
))
83 cst
= isl_upoly_as_cst(up
);
87 return isl_int_is_zero(cst
->n
) && isl_int_is_pos(cst
->d
);
90 int isl_upoly_is_nan(__isl_keep
struct isl_upoly
*up
)
92 struct isl_upoly_cst
*cst
;
96 if (!isl_upoly_is_cst(up
))
99 cst
= isl_upoly_as_cst(up
);
103 return isl_int_is_zero(cst
->n
) && isl_int_is_zero(cst
->d
);
106 int isl_upoly_is_infty(__isl_keep
struct isl_upoly
*up
)
108 struct isl_upoly_cst
*cst
;
112 if (!isl_upoly_is_cst(up
))
115 cst
= isl_upoly_as_cst(up
);
119 return isl_int_is_pos(cst
->n
) && isl_int_is_zero(cst
->d
);
122 int isl_upoly_is_neginfty(__isl_keep
struct isl_upoly
*up
)
124 struct isl_upoly_cst
*cst
;
128 if (!isl_upoly_is_cst(up
))
131 cst
= isl_upoly_as_cst(up
);
135 return isl_int_is_neg(cst
->n
) && isl_int_is_zero(cst
->d
);
138 int isl_upoly_is_one(__isl_keep
struct isl_upoly
*up
)
140 struct isl_upoly_cst
*cst
;
144 if (!isl_upoly_is_cst(up
))
147 cst
= isl_upoly_as_cst(up
);
151 return isl_int_eq(cst
->n
, cst
->d
) && isl_int_is_pos(cst
->d
);
154 int isl_upoly_is_negone(__isl_keep
struct isl_upoly
*up
)
156 struct isl_upoly_cst
*cst
;
160 if (!isl_upoly_is_cst(up
))
163 cst
= isl_upoly_as_cst(up
);
167 return isl_int_is_negone(cst
->n
) && isl_int_is_one(cst
->d
);
170 __isl_give
struct isl_upoly_cst
*isl_upoly_cst_alloc(struct isl_ctx
*ctx
)
172 struct isl_upoly_cst
*cst
;
174 cst
= isl_alloc_type(ctx
, struct isl_upoly_cst
);
183 isl_int_init(cst
->n
);
184 isl_int_init(cst
->d
);
189 __isl_give
struct isl_upoly
*isl_upoly_zero(struct isl_ctx
*ctx
)
191 struct isl_upoly_cst
*cst
;
193 cst
= isl_upoly_cst_alloc(ctx
);
197 isl_int_set_si(cst
->n
, 0);
198 isl_int_set_si(cst
->d
, 1);
203 __isl_give
struct isl_upoly
*isl_upoly_infty(struct isl_ctx
*ctx
)
205 struct isl_upoly_cst
*cst
;
207 cst
= isl_upoly_cst_alloc(ctx
);
211 isl_int_set_si(cst
->n
, 1);
212 isl_int_set_si(cst
->d
, 0);
217 __isl_give
struct isl_upoly
*isl_upoly_nan(struct isl_ctx
*ctx
)
219 struct isl_upoly_cst
*cst
;
221 cst
= isl_upoly_cst_alloc(ctx
);
225 isl_int_set_si(cst
->n
, 0);
226 isl_int_set_si(cst
->d
, 0);
231 __isl_give
struct isl_upoly
*isl_upoly_rat_cst(struct isl_ctx
*ctx
,
232 isl_int n
, isl_int d
)
234 struct isl_upoly_cst
*cst
;
236 cst
= isl_upoly_cst_alloc(ctx
);
240 isl_int_set(cst
->n
, n
);
241 isl_int_set(cst
->d
, d
);
246 __isl_give
struct isl_upoly_rec
*isl_upoly_alloc_rec(struct isl_ctx
*ctx
,
249 struct isl_upoly_rec
*rec
;
251 isl_assert(ctx
, var
>= 0, return NULL
);
252 isl_assert(ctx
, size
>= 0, return NULL
);
253 rec
= isl_calloc(dim
->ctx
, struct isl_upoly_rec
,
254 sizeof(struct isl_upoly_rec
) +
255 (size
- 1) * sizeof(struct isl_upoly
*));
269 isl_upoly_free(&rec
->up
);
273 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial
*qp
)
275 struct isl_upoly_cst
*cst
;
280 return isl_upoly_is_zero(qp
->upoly
);
283 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial
*qp
)
285 struct isl_upoly_cst
*cst
;
290 return isl_upoly_is_one(qp
->upoly
);
293 static void upoly_free_cst(__isl_take
struct isl_upoly_cst
*cst
)
295 isl_int_clear(cst
->n
);
296 isl_int_clear(cst
->d
);
299 static void upoly_free_rec(__isl_take
struct isl_upoly_rec
*rec
)
303 for (i
= 0; i
< rec
->n
; ++i
)
304 isl_upoly_free(rec
->p
[i
]);
307 __isl_give
struct isl_upoly
*isl_upoly_copy(__isl_keep
struct isl_upoly
*up
)
316 __isl_give
struct isl_upoly
*isl_upoly_dup_cst(__isl_keep
struct isl_upoly
*up
)
318 struct isl_upoly_cst
*cst
;
319 struct isl_upoly_cst
*dup
;
321 cst
= isl_upoly_as_cst(up
);
325 dup
= isl_upoly_as_cst(isl_upoly_zero(up
->ctx
));
328 isl_int_set(dup
->n
, cst
->n
);
329 isl_int_set(dup
->d
, cst
->d
);
334 __isl_give
struct isl_upoly
*isl_upoly_dup_rec(__isl_keep
struct isl_upoly
*up
)
337 struct isl_upoly_rec
*rec
;
338 struct isl_upoly_rec
*dup
;
340 rec
= isl_upoly_as_rec(up
);
344 dup
= isl_upoly_alloc_rec(up
->ctx
, up
->var
, rec
->n
);
348 for (i
= 0; i
< rec
->n
; ++i
) {
349 dup
->p
[i
] = isl_upoly_copy(rec
->p
[i
]);
357 isl_upoly_free(&dup
->up
);
361 __isl_give
struct isl_upoly
*isl_upoly_dup(__isl_keep
struct isl_upoly
*up
)
363 struct isl_upoly
*dup
;
368 if (isl_upoly_is_cst(up
))
369 return isl_upoly_dup_cst(up
);
371 return isl_upoly_dup_rec(up
);
374 __isl_give
struct isl_upoly
*isl_upoly_cow(__isl_take
struct isl_upoly
*up
)
382 return isl_upoly_dup(up
);
385 void isl_upoly_free(__isl_take
struct isl_upoly
*up
)
394 upoly_free_cst((struct isl_upoly_cst
*)up
);
396 upoly_free_rec((struct isl_upoly_rec
*)up
);
398 isl_ctx_deref(up
->ctx
);
402 static void isl_upoly_cst_reduce(__isl_keep
struct isl_upoly_cst
*cst
)
407 isl_int_gcd(gcd
, cst
->n
, cst
->d
);
408 if (!isl_int_is_zero(gcd
) && !isl_int_is_one(gcd
)) {
409 isl_int_divexact(cst
->n
, cst
->n
, gcd
);
410 isl_int_divexact(cst
->d
, cst
->d
, gcd
);
415 __isl_give
struct isl_upoly
*isl_upoly_sum_cst(__isl_take
struct isl_upoly
*up1
,
416 __isl_take
struct isl_upoly
*up2
)
418 struct isl_upoly_cst
*cst1
;
419 struct isl_upoly_cst
*cst2
;
421 up1
= isl_upoly_cow(up1
);
425 cst1
= isl_upoly_as_cst(up1
);
426 cst2
= isl_upoly_as_cst(up2
);
428 if (isl_int_eq(cst1
->d
, cst2
->d
))
429 isl_int_add(cst1
->n
, cst1
->n
, cst2
->n
);
431 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->d
);
432 isl_int_addmul(cst1
->n
, cst2
->n
, cst1
->d
);
433 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
436 isl_upoly_cst_reduce(cst1
);
446 static __isl_give
struct isl_upoly
*replace_by_zero(
447 __isl_take
struct isl_upoly
*up
)
455 return isl_upoly_zero(ctx
);
458 static __isl_give
struct isl_upoly
*replace_by_constant_term(
459 __isl_take
struct isl_upoly
*up
)
461 struct isl_upoly_rec
*rec
;
462 struct isl_upoly
*cst
;
467 rec
= isl_upoly_as_rec(up
);
470 isl_assert(up
->ctx
, rec
->n
== 1, goto error
);
471 cst
= isl_upoly_copy(rec
->p
[0]);
479 __isl_give
struct isl_upoly
*isl_upoly_sum(__isl_take
struct isl_upoly
*up1
,
480 __isl_take
struct isl_upoly
*up2
)
483 struct isl_upoly_rec
*rec1
, *rec2
;
488 if (isl_upoly_is_nan(up1
)) {
493 if (isl_upoly_is_nan(up2
)) {
498 if (isl_upoly_is_zero(up1
)) {
503 if (isl_upoly_is_zero(up2
)) {
508 if (up1
->var
< up2
->var
)
509 return isl_upoly_sum(up2
, up1
);
511 if (up2
->var
< up1
->var
) {
512 struct isl_upoly_rec
*rec
;
513 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
517 up1
= isl_upoly_cow(up1
);
518 rec
= isl_upoly_as_rec(up1
);
521 rec
->p
[0] = isl_upoly_sum(rec
->p
[0], up2
);
523 up1
= replace_by_constant_term(up1
);
527 if (isl_upoly_is_cst(up1
))
528 return isl_upoly_sum_cst(up1
, up2
);
530 rec1
= isl_upoly_as_rec(up1
);
531 rec2
= isl_upoly_as_rec(up2
);
535 if (rec1
->n
< rec2
->n
)
536 return isl_upoly_sum(up2
, up1
);
538 up1
= isl_upoly_cow(up1
);
539 rec1
= isl_upoly_as_rec(up1
);
543 for (i
= rec2
->n
- 1; i
>= 0; --i
) {
544 rec1
->p
[i
] = isl_upoly_sum(rec1
->p
[i
],
545 isl_upoly_copy(rec2
->p
[i
]));
548 if (i
== rec1
->n
- 1 && isl_upoly_is_zero(rec1
->p
[i
])) {
549 isl_upoly_free(rec1
->p
[i
]);
555 up1
= replace_by_zero(up1
);
556 else if (rec1
->n
== 1)
557 up1
= replace_by_constant_term(up1
);
568 __isl_give
struct isl_upoly
*isl_upoly_neg_cst(__isl_take
struct isl_upoly
*up
)
570 struct isl_upoly_cst
*cst
;
572 if (isl_upoly_is_zero(up
))
575 up
= isl_upoly_cow(up
);
579 cst
= isl_upoly_as_cst(up
);
581 isl_int_neg(cst
->n
, cst
->n
);
586 __isl_give
struct isl_upoly
*isl_upoly_neg(__isl_take
struct isl_upoly
*up
)
589 struct isl_upoly_rec
*rec
;
594 if (isl_upoly_is_cst(up
))
595 return isl_upoly_neg_cst(up
);
597 up
= isl_upoly_cow(up
);
598 rec
= isl_upoly_as_rec(up
);
602 for (i
= 0; i
< rec
->n
; ++i
) {
603 rec
->p
[i
] = isl_upoly_neg(rec
->p
[i
]);
614 __isl_give
struct isl_upoly
*isl_upoly_mul_cst(__isl_take
struct isl_upoly
*up1
,
615 __isl_take
struct isl_upoly
*up2
)
617 struct isl_upoly_cst
*cst1
;
618 struct isl_upoly_cst
*cst2
;
620 up1
= isl_upoly_cow(up1
);
624 cst1
= isl_upoly_as_cst(up1
);
625 cst2
= isl_upoly_as_cst(up2
);
627 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->n
);
628 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
630 isl_upoly_cst_reduce(cst1
);
640 __isl_give
struct isl_upoly
*isl_upoly_mul_rec(__isl_take
struct isl_upoly
*up1
,
641 __isl_take
struct isl_upoly
*up2
)
643 struct isl_upoly_rec
*rec1
;
644 struct isl_upoly_rec
*rec2
;
645 struct isl_upoly_rec
*res
;
649 rec1
= isl_upoly_as_rec(up1
);
650 rec2
= isl_upoly_as_rec(up2
);
653 size
= rec1
->n
+ rec2
->n
- 1;
654 res
= isl_upoly_alloc_rec(up1
->ctx
, up1
->var
, size
);
658 for (i
= 0; i
< rec1
->n
; ++i
) {
659 res
->p
[i
] = isl_upoly_mul(isl_upoly_copy(rec2
->p
[0]),
660 isl_upoly_copy(rec1
->p
[i
]));
665 for (; i
< size
; ++i
) {
666 res
->p
[i
] = isl_upoly_zero(up1
->ctx
);
671 for (i
= 0; i
< rec1
->n
; ++i
) {
672 for (j
= 1; j
< rec2
->n
; ++j
) {
673 struct isl_upoly
*up
;
674 up
= isl_upoly_mul(isl_upoly_copy(rec2
->p
[j
]),
675 isl_upoly_copy(rec1
->p
[i
]));
676 res
->p
[i
+ j
] = isl_upoly_sum(res
->p
[i
+ j
], up
);
689 isl_upoly_free(&res
->up
);
693 __isl_give
struct isl_upoly
*isl_upoly_mul(__isl_take
struct isl_upoly
*up1
,
694 __isl_take
struct isl_upoly
*up2
)
699 if (isl_upoly_is_nan(up1
)) {
704 if (isl_upoly_is_nan(up2
)) {
709 if (isl_upoly_is_zero(up1
)) {
714 if (isl_upoly_is_zero(up2
)) {
719 if (isl_upoly_is_one(up1
)) {
724 if (isl_upoly_is_one(up2
)) {
729 if (up1
->var
< up2
->var
)
730 return isl_upoly_mul(up2
, up1
);
732 if (up2
->var
< up1
->var
) {
734 struct isl_upoly_rec
*rec
;
735 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
736 isl_ctx
*ctx
= up1
->ctx
;
739 return isl_upoly_nan(ctx
);
741 up1
= isl_upoly_cow(up1
);
742 rec
= isl_upoly_as_rec(up1
);
746 for (i
= 0; i
< rec
->n
; ++i
) {
747 rec
->p
[i
] = isl_upoly_mul(rec
->p
[i
],
748 isl_upoly_copy(up2
));
756 if (isl_upoly_is_cst(up1
))
757 return isl_upoly_mul_cst(up1
, up2
);
759 return isl_upoly_mul_rec(up1
, up2
);
766 __isl_give isl_qpolynomial
*isl_qpolynomial_alloc(__isl_take isl_dim
*dim
,
769 struct isl_qpolynomial
*qp
;
775 total
= isl_dim_total(dim
);
777 qp
= isl_calloc_type(dim
->ctx
, struct isl_qpolynomial
);
782 qp
->div
= isl_mat_alloc(dim
->ctx
, n_div
, 1 + 1 + total
+ n_div
);
791 isl_qpolynomial_free(qp
);
795 __isl_give isl_qpolynomial
*isl_qpolynomial_copy(__isl_keep isl_qpolynomial
*qp
)
804 __isl_give isl_qpolynomial
*isl_qpolynomial_dup(__isl_keep isl_qpolynomial
*qp
)
806 struct isl_qpolynomial
*dup
;
811 dup
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
);
814 isl_mat_free(dup
->div
);
815 dup
->div
= isl_mat_copy(qp
->div
);
818 dup
->upoly
= isl_upoly_copy(qp
->upoly
);
824 isl_qpolynomial_free(dup
);
828 __isl_give isl_qpolynomial
*isl_qpolynomial_cow(__isl_take isl_qpolynomial
*qp
)
836 return isl_qpolynomial_dup(qp
);
839 void isl_qpolynomial_free(__isl_take isl_qpolynomial
*qp
)
847 isl_dim_free(qp
->dim
);
848 isl_mat_free(qp
->div
);
849 isl_upoly_free(qp
->upoly
);
854 static int compatible_divs(__isl_keep isl_mat
*div1
, __isl_keep isl_mat
*div2
)
859 isl_assert(div1
->ctx
, div1
->n_row
>= div2
->n_row
&&
860 div1
->n_col
>= div2
->n_col
, return -1);
862 if (div1
->n_row
== div2
->n_row
)
863 return isl_mat_is_equal(div1
, div2
);
867 div1
->n_row
= div2
->n_row
;
868 div1
->n_col
= div2
->n_col
;
870 equal
= isl_mat_is_equal(div1
, div2
);
878 static void expand_row(__isl_keep isl_mat
*dst
, int d
,
879 __isl_keep isl_mat
*src
, int s
, int *exp
)
882 unsigned c
= src
->n_col
- src
->n_row
;
884 isl_seq_cpy(dst
->row
[d
], src
->row
[s
], c
);
885 isl_seq_clr(dst
->row
[d
] + c
, dst
->n_col
- c
);
887 for (i
= 0; i
< s
; ++i
)
888 isl_int_set(dst
->row
[d
][c
+ exp
[i
]], src
->row
[s
][c
+ i
]);
891 static int cmp_row(__isl_keep isl_mat
*div
, int i
, int j
)
895 li
= isl_seq_last_non_zero(div
->row
[i
], div
->n_col
);
896 lj
= isl_seq_last_non_zero(div
->row
[j
], div
->n_col
);
901 return isl_seq_cmp(div
->row
[i
], div
->row
[j
], div
->n_col
);
904 static __isl_give isl_mat
*merge_divs(__isl_keep isl_mat
*div1
,
905 __isl_keep isl_mat
*div2
, int *exp1
, int *exp2
)
909 unsigned d
= div1
->n_col
- div1
->n_row
;
911 div
= isl_mat_alloc(div1
->ctx
, 1 + div1
->n_row
+ div2
->n_row
,
912 d
+ div1
->n_row
+ div2
->n_row
);
916 for (i
= 0, j
= 0, k
= 0; i
< div1
->n_row
&& j
< div2
->n_row
; ++k
) {
919 expand_row(div
, k
, div1
, i
, exp1
);
920 expand_row(div
, k
+ 1, div2
, j
, exp2
);
922 cmp
= cmp_row(div
, k
, k
+ 1);
926 } else if (cmp
< 0) {
930 isl_seq_cpy(div
->row
[k
], div
->row
[k
+ 1], div
->n_col
);
933 for (; i
< div1
->n_row
; ++i
, ++k
) {
934 expand_row(div
, k
, div1
, i
, exp1
);
937 for (; j
< div2
->n_row
; ++j
, ++k
) {
938 expand_row(div
, k
, div2
, j
, exp2
);
948 static __isl_give
struct isl_upoly
*expand(__isl_take
struct isl_upoly
*up
,
952 struct isl_upoly_rec
*rec
;
954 if (isl_upoly_is_cst(up
))
960 if (exp
[up
->var
- first
] == up
->var
- first
)
963 up
= isl_upoly_cow(up
);
967 up
->var
= exp
[up
->var
- first
] + first
;
969 rec
= isl_upoly_as_rec(up
);
973 for (i
= 0; i
< rec
->n
; ++i
) {
974 rec
->p
[i
] = expand(rec
->p
[i
], exp
, first
);
985 static __isl_give isl_qpolynomial
*with_merged_divs(
986 __isl_give isl_qpolynomial
*(*fn
)(__isl_take isl_qpolynomial
*qp1
,
987 __isl_take isl_qpolynomial
*qp2
),
988 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
994 qp1
= isl_qpolynomial_cow(qp1
);
995 qp2
= isl_qpolynomial_cow(qp2
);
1000 isl_assert(qp1
->div
->ctx
, qp1
->div
->n_row
>= qp2
->div
->n_row
&&
1001 qp1
->div
->n_col
>= qp2
->div
->n_col
, goto error
);
1003 exp1
= isl_alloc_array(qp1
->div
->ctx
, int, qp1
->div
->n_row
);
1004 exp2
= isl_alloc_array(qp2
->div
->ctx
, int, qp2
->div
->n_row
);
1008 div
= merge_divs(qp1
->div
, qp2
->div
, exp1
, exp2
);
1012 isl_mat_free(qp1
->div
);
1013 qp1
->div
= isl_mat_copy(div
);
1014 isl_mat_free(qp2
->div
);
1015 qp2
->div
= isl_mat_copy(div
);
1017 qp1
->upoly
= expand(qp1
->upoly
, exp1
, div
->n_col
- div
->n_row
- 2);
1018 qp2
->upoly
= expand(qp2
->upoly
, exp2
, div
->n_col
- div
->n_row
- 2);
1020 if (!qp1
->upoly
|| !qp2
->upoly
)
1027 return fn(qp1
, qp2
);
1032 isl_qpolynomial_free(qp1
);
1033 isl_qpolynomial_free(qp2
);
1037 __isl_give isl_qpolynomial
*isl_qpolynomial_add(__isl_take isl_qpolynomial
*qp1
,
1038 __isl_take isl_qpolynomial
*qp2
)
1040 qp1
= isl_qpolynomial_cow(qp1
);
1045 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1046 return isl_qpolynomial_add(qp2
, qp1
);
1048 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1049 if (!compatible_divs(qp1
->div
, qp2
->div
))
1050 return with_merged_divs(isl_qpolynomial_add
, qp1
, qp2
);
1052 qp1
->upoly
= isl_upoly_sum(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1056 isl_qpolynomial_free(qp2
);
1060 isl_qpolynomial_free(qp1
);
1061 isl_qpolynomial_free(qp2
);
1065 __isl_give isl_qpolynomial
*isl_qpolynomial_sub(__isl_take isl_qpolynomial
*qp1
,
1066 __isl_take isl_qpolynomial
*qp2
)
1068 return isl_qpolynomial_add(qp1
, isl_qpolynomial_neg(qp2
));
1071 __isl_give isl_qpolynomial
*isl_qpolynomial_neg(__isl_take isl_qpolynomial
*qp
)
1073 qp
= isl_qpolynomial_cow(qp
);
1078 qp
->upoly
= isl_upoly_neg(qp
->upoly
);
1084 isl_qpolynomial_free(qp
);
1088 __isl_give isl_qpolynomial
*isl_qpolynomial_mul(__isl_take isl_qpolynomial
*qp1
,
1089 __isl_take isl_qpolynomial
*qp2
)
1091 qp1
= isl_qpolynomial_cow(qp1
);
1096 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1097 return isl_qpolynomial_mul(qp2
, qp1
);
1099 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1100 if (!compatible_divs(qp1
->div
, qp2
->div
))
1101 return with_merged_divs(isl_qpolynomial_mul
, qp1
, qp2
);
1103 qp1
->upoly
= isl_upoly_mul(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1107 isl_qpolynomial_free(qp2
);
1111 isl_qpolynomial_free(qp1
);
1112 isl_qpolynomial_free(qp2
);
1116 __isl_give isl_qpolynomial
*isl_qpolynomial_zero(__isl_take isl_dim
*dim
)
1118 struct isl_qpolynomial
*qp
;
1119 struct isl_upoly_cst
*cst
;
1121 qp
= isl_qpolynomial_alloc(dim
, 0);
1125 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1131 isl_qpolynomial_free(qp
);
1135 __isl_give isl_qpolynomial
*isl_qpolynomial_infty(__isl_take isl_dim
*dim
)
1137 struct isl_qpolynomial
*qp
;
1138 struct isl_upoly_cst
*cst
;
1140 qp
= isl_qpolynomial_alloc(dim
, 0);
1144 qp
->upoly
= isl_upoly_infty(dim
->ctx
);
1150 isl_qpolynomial_free(qp
);
1154 __isl_give isl_qpolynomial
*isl_qpolynomial_nan(__isl_take isl_dim
*dim
)
1156 struct isl_qpolynomial
*qp
;
1157 struct isl_upoly_cst
*cst
;
1159 qp
= isl_qpolynomial_alloc(dim
, 0);
1163 qp
->upoly
= isl_upoly_nan(dim
->ctx
);
1169 isl_qpolynomial_free(qp
);
1173 __isl_give isl_qpolynomial
*isl_qpolynomial_cst(__isl_take isl_dim
*dim
,
1176 struct isl_qpolynomial
*qp
;
1177 struct isl_upoly_cst
*cst
;
1179 qp
= isl_qpolynomial_alloc(dim
, 0);
1183 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1186 cst
= isl_upoly_as_cst(qp
->upoly
);
1187 isl_int_set(cst
->n
, v
);
1191 isl_qpolynomial_free(qp
);
1195 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial
*qp
,
1196 isl_int
*n
, isl_int
*d
)
1198 struct isl_upoly_cst
*cst
;
1203 if (!isl_upoly_is_cst(qp
->upoly
))
1206 cst
= isl_upoly_as_cst(qp
->upoly
);
1211 isl_int_set(*n
, cst
->n
);
1213 isl_int_set(*d
, cst
->d
);
1218 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial
*qp1
,
1219 __isl_keep isl_qpolynomial
*qp2
)
1224 return isl_upoly_is_equal(qp1
->upoly
, qp2
->upoly
);
1227 static void upoly_update_den(__isl_keep
struct isl_upoly
*up
, isl_int
*d
)
1230 struct isl_upoly_rec
*rec
;
1232 if (isl_upoly_is_cst(up
)) {
1233 struct isl_upoly_cst
*cst
;
1234 cst
= isl_upoly_as_cst(up
);
1237 isl_int_lcm(*d
, *d
, cst
->d
);
1241 rec
= isl_upoly_as_rec(up
);
1245 for (i
= 0; i
< rec
->n
; ++i
)
1246 upoly_update_den(rec
->p
[i
], d
);
1249 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial
*qp
, isl_int
*d
)
1251 isl_int_set_si(*d
, 1);
1254 upoly_update_den(qp
->upoly
, d
);
1257 __isl_give isl_qpolynomial
*isl_qpolynomial_pow(__isl_take isl_dim
*dim
,
1261 struct isl_qpolynomial
*qp
;
1262 struct isl_upoly_rec
*rec
;
1263 struct isl_upoly_cst
*cst
;
1264 struct isl_ctx
*ctx
;
1271 qp
= isl_qpolynomial_alloc(dim
, 0);
1275 qp
->upoly
= &isl_upoly_alloc_rec(ctx
, pos
, 1 + power
)->up
;
1278 rec
= isl_upoly_as_rec(qp
->upoly
);
1279 for (i
= 0; i
< 1 + power
; ++i
) {
1280 rec
->p
[i
] = isl_upoly_zero(ctx
);
1285 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1286 isl_int_set_si(cst
->n
, 1);
1290 isl_qpolynomial_free(qp
);
1294 __isl_give isl_qpolynomial
*isl_qpolynomial_var(__isl_take isl_dim
*dim
,
1295 enum isl_dim_type type
, unsigned pos
)
1300 isl_assert(dim
->ctx
, isl_dim_size(dim
, isl_dim_in
) == 0, goto error
);
1301 isl_assert(dim
->ctx
, pos
< isl_dim_size(dim
, type
), goto error
);
1303 if (type
== isl_dim_set
)
1304 pos
+= isl_dim_size(dim
, isl_dim_param
);
1306 return isl_qpolynomial_pow(dim
, pos
, 1);
1312 __isl_give isl_qpolynomial
*isl_qpolynomial_div_pow(__isl_take isl_div
*div
,
1315 struct isl_qpolynomial
*qp
= NULL
;
1316 struct isl_upoly_rec
*rec
;
1317 struct isl_upoly_cst
*cst
;
1323 isl_assert(div
->ctx
, div
->bmap
->n_div
== 1, goto error
);
1325 qp
= isl_qpolynomial_alloc(isl_basic_map_get_dim(div
->bmap
), 1);
1329 isl_seq_cpy(qp
->div
->row
[0], div
->line
[0], qp
->div
->n_col
- 1);
1330 isl_int_set_si(qp
->div
->row
[0][qp
->div
->n_col
- 1], 0);
1332 pos
= isl_dim_total(qp
->dim
);
1333 qp
->upoly
= &isl_upoly_alloc_rec(div
->ctx
, pos
, 1 + power
)->up
;
1336 rec
= isl_upoly_as_rec(qp
->upoly
);
1337 for (i
= 0; i
< 1 + power
; ++i
) {
1338 rec
->p
[i
] = isl_upoly_zero(div
->ctx
);
1343 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1344 isl_int_set_si(cst
->n
, 1);
1350 isl_qpolynomial_free(qp
);
1355 __isl_give isl_qpolynomial
*isl_qpolynomial_div(__isl_take isl_div
*div
)
1357 return isl_qpolynomial_div_pow(div
, 1);
1360 __isl_give isl_qpolynomial
*isl_qpolynomial_rat_cst(__isl_take isl_dim
*dim
,
1361 const isl_int n
, const isl_int d
)
1363 struct isl_qpolynomial
*qp
;
1364 struct isl_upoly_cst
*cst
;
1366 qp
= isl_qpolynomial_alloc(dim
, 0);
1370 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1373 cst
= isl_upoly_as_cst(qp
->upoly
);
1374 isl_int_set(cst
->n
, n
);
1375 isl_int_set(cst
->d
, d
);
1379 isl_qpolynomial_free(qp
);
1383 static int up_set_active(__isl_keep
struct isl_upoly
*up
, int *active
, int d
)
1385 struct isl_upoly_rec
*rec
;
1391 if (isl_upoly_is_cst(up
))
1395 active
[up
->var
] = 1;
1397 rec
= isl_upoly_as_rec(up
);
1398 for (i
= 0; i
< rec
->n
; ++i
)
1399 if (up_set_active(rec
->p
[i
], active
, d
) < 0)
1405 static int set_active(__isl_keep isl_qpolynomial
*qp
, int *active
)
1408 int d
= isl_dim_total(qp
->dim
);
1413 for (i
= 0; i
< d
; ++i
)
1414 for (j
= 0; j
< qp
->div
->n_row
; ++j
) {
1415 if (isl_int_is_zero(qp
->div
->row
[j
][2 + i
]))
1421 return up_set_active(qp
->upoly
, active
, d
);
1424 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial
*qp
,
1425 enum isl_dim_type type
, unsigned first
, unsigned n
)
1436 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
1438 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
1439 type
== isl_dim_set
, return -1);
1441 active
= isl_calloc_array(set
->ctx
, int, isl_dim_total(qp
->dim
));
1442 if (set_active(qp
, active
) < 0)
1445 if (type
== isl_dim_set
)
1446 first
+= isl_dim_size(qp
->dim
, isl_dim_param
);
1447 for (i
= 0; i
< n
; ++i
)
1448 if (active
[first
+ i
]) {
1461 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold
*fold
,
1462 enum isl_dim_type type
, unsigned first
, unsigned n
)
1468 if (fold
->n
== 0 || n
== 0)
1471 for (i
= 0; i
< fold
->n
; ++i
) {
1472 int involves
= isl_qpolynomial_involves_dims(fold
->qp
[i
],
1474 if (involves
< 0 || involves
)
1481 #define PW isl_pw_qpolynomial
1483 #define EL isl_qpolynomial
1485 #define IS_ZERO is_zero
1491 #include <isl_pw_templ.c>
1494 #define PW isl_pw_qpolynomial_fold
1496 #define EL isl_qpolynomial_fold
1498 #define IS_ZERO is_empty
1504 #include <isl_pw_templ.c>
1506 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial
*pwqp
)
1514 if (!isl_set_fast_is_universe(pwqp
->p
[0].set
))
1517 return isl_qpolynomial_is_one(pwqp
->p
[0].qp
);
1520 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_mul(
1521 __isl_take isl_pw_qpolynomial
*pwqp1
,
1522 __isl_take isl_pw_qpolynomial
*pwqp2
)
1525 struct isl_pw_qpolynomial
*res
;
1528 if (!pwqp1
|| !pwqp2
)
1531 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1534 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1535 isl_pw_qpolynomial_free(pwqp2
);
1539 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1540 isl_pw_qpolynomial_free(pwqp1
);
1544 if (isl_pw_qpolynomial_is_one(pwqp1
)) {
1545 isl_pw_qpolynomial_free(pwqp1
);
1549 if (isl_pw_qpolynomial_is_one(pwqp2
)) {
1550 isl_pw_qpolynomial_free(pwqp2
);
1554 n
= pwqp1
->n
* pwqp2
->n
;
1555 res
= isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1
->dim
), n
);
1557 for (i
= 0; i
< pwqp1
->n
; ++i
) {
1558 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1559 struct isl_set
*common
;
1560 struct isl_qpolynomial
*prod
;
1561 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
1562 isl_set_copy(pwqp2
->p
[j
].set
));
1563 if (isl_set_fast_is_empty(common
)) {
1564 isl_set_free(common
);
1568 prod
= isl_qpolynomial_mul(
1569 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
1570 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1572 res
= isl_pw_qpolynomial_add_piece(res
, common
, prod
);
1576 isl_pw_qpolynomial_free(pwqp1
);
1577 isl_pw_qpolynomial_free(pwqp2
);
1581 isl_pw_qpolynomial_free(pwqp1
);
1582 isl_pw_qpolynomial_free(pwqp2
);
1586 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_neg(
1587 __isl_take isl_pw_qpolynomial
*pwqp
)
1590 struct isl_pw_qpolynomial
*res
;
1596 if (isl_pw_qpolynomial_is_zero(pwqp
))
1599 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
1601 for (i
= 0; i
< pwqp
->n
; ++i
) {
1602 pwqp
->p
[i
].qp
= isl_qpolynomial_neg(pwqp
->p
[i
].qp
);
1609 isl_pw_qpolynomial_free(pwqp
);
1613 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_sub(
1614 __isl_take isl_pw_qpolynomial
*pwqp1
,
1615 __isl_take isl_pw_qpolynomial
*pwqp2
)
1617 return isl_pw_qpolynomial_add(pwqp1
, isl_pw_qpolynomial_neg(pwqp2
));
1620 __isl_give
struct isl_upoly
*isl_upoly_eval(
1621 __isl_take
struct isl_upoly
*up
, __isl_take isl_vec
*vec
)
1624 struct isl_upoly_rec
*rec
;
1625 struct isl_upoly
*res
;
1626 struct isl_upoly
*base
;
1628 if (isl_upoly_is_cst(up
)) {
1633 rec
= isl_upoly_as_rec(up
);
1637 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1639 base
= isl_upoly_rat_cst(up
->ctx
, vec
->el
[1 + up
->var
], vec
->el
[0]);
1641 res
= isl_upoly_eval(isl_upoly_copy(rec
->p
[rec
->n
- 1]),
1644 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1645 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1646 res
= isl_upoly_sum(res
,
1647 isl_upoly_eval(isl_upoly_copy(rec
->p
[i
]),
1648 isl_vec_copy(vec
)));
1651 isl_upoly_free(base
);
1661 __isl_give isl_qpolynomial
*isl_qpolynomial_eval(
1662 __isl_take isl_qpolynomial
*qp
, __isl_take isl_point
*pnt
)
1665 struct isl_upoly
*up
;
1670 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, qp
->dim
), goto error
);
1672 if (qp
->div
->n_row
== 0)
1673 ext
= isl_vec_copy(pnt
->vec
);
1676 unsigned dim
= isl_dim_total(qp
->dim
);
1677 ext
= isl_vec_alloc(qp
->dim
->ctx
, 1 + dim
+ qp
->div
->n_row
);
1681 isl_seq_cpy(ext
->el
, pnt
->vec
->el
, pnt
->vec
->size
);
1682 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1683 isl_seq_inner_product(qp
->div
->row
[i
] + 1, ext
->el
,
1684 1 + dim
+ i
, &ext
->el
[1+dim
+i
]);
1685 isl_int_fdiv_q(ext
->el
[1+dim
+i
], ext
->el
[1+dim
+i
],
1686 qp
->div
->row
[i
][0]);
1690 up
= isl_upoly_eval(isl_upoly_copy(qp
->upoly
), ext
);
1694 dim
= isl_dim_copy(qp
->dim
);
1695 isl_qpolynomial_free(qp
);
1696 isl_point_free(pnt
);
1698 qp
= isl_qpolynomial_alloc(dim
, 0);
1706 isl_qpolynomial_free(qp
);
1707 isl_point_free(pnt
);
1711 int isl_upoly_cmp(__isl_keep
struct isl_upoly_cst
*cst1
,
1712 __isl_keep
struct isl_upoly_cst
*cst2
)
1717 isl_int_mul(t
, cst1
->n
, cst2
->d
);
1718 isl_int_submul(t
, cst2
->n
, cst1
->d
);
1719 cmp
= isl_int_sgn(t
);
1724 static __isl_give isl_qpolynomial
*qpolynomial_min(
1725 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1727 struct isl_upoly_cst
*cst1
, *cst2
;
1732 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
1733 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
1734 cst1
= isl_upoly_as_cst(qp1
->upoly
);
1735 cst2
= isl_upoly_as_cst(qp2
->upoly
);
1736 cmp
= isl_upoly_cmp(cst1
, cst2
);
1739 isl_qpolynomial_free(qp2
);
1741 isl_qpolynomial_free(qp1
);
1746 isl_qpolynomial_free(qp1
);
1747 isl_qpolynomial_free(qp2
);
1751 static __isl_give isl_qpolynomial
*qpolynomial_max(
1752 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1754 struct isl_upoly_cst
*cst1
, *cst2
;
1759 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
1760 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
1761 cst1
= isl_upoly_as_cst(qp1
->upoly
);
1762 cst2
= isl_upoly_as_cst(qp2
->upoly
);
1763 cmp
= isl_upoly_cmp(cst1
, cst2
);
1766 isl_qpolynomial_free(qp2
);
1768 isl_qpolynomial_free(qp1
);
1773 isl_qpolynomial_free(qp1
);
1774 isl_qpolynomial_free(qp2
);
1778 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold
*fold1
,
1779 __isl_keep isl_qpolynomial_fold
*fold2
)
1783 if (!fold1
|| !fold2
)
1786 if (fold1
->n
!= fold2
->n
)
1789 /* We probably want to sort the qps first... */
1790 for (i
= 0; i
< fold1
->n
; ++i
) {
1791 int eq
= isl_qpolynomial_is_equal(fold1
->qp
[i
], fold2
->qp
[i
]);
1799 __isl_give isl_qpolynomial
*isl_qpolynomial_fold_eval(
1800 __isl_take isl_qpolynomial_fold
*fold
, __isl_take isl_point
*pnt
)
1802 isl_qpolynomial
*qp
;
1806 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, fold
->dim
), goto error
);
1807 isl_assert(pnt
->dim
->ctx
,
1808 fold
->type
== isl_fold_max
|| fold
->type
== isl_fold_min
,
1812 qp
= isl_qpolynomial_zero(isl_dim_copy(fold
->dim
));
1815 qp
= isl_qpolynomial_eval(isl_qpolynomial_copy(fold
->qp
[0]),
1816 isl_point_copy(pnt
));
1817 for (i
= 1; i
< fold
->n
; ++i
) {
1818 isl_qpolynomial
*qp_i
;
1819 qp_i
= isl_qpolynomial_eval(
1820 isl_qpolynomial_copy(fold
->qp
[i
]),
1821 isl_point_copy(pnt
));
1822 if (fold
->type
== isl_fold_max
)
1823 qp
= qpolynomial_max(qp
, qp_i
);
1825 qp
= qpolynomial_min(qp
, qp_i
);
1828 isl_qpolynomial_fold_free(fold
);
1829 isl_point_free(pnt
);
1833 isl_qpolynomial_fold_free(fold
);
1834 isl_point_free(pnt
);
1838 __isl_give isl_qpolynomial
*isl_qpolynomial_move(__isl_take isl_qpolynomial
*qp
,
1839 enum isl_dim_type dst_type
, unsigned dst_pos
,
1840 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
1845 qp
->dim
= isl_dim_move(qp
->dim
, dst_type
, dst_pos
, src_type
, src_pos
, n
);
1849 /* Do something to polynomials when needed; later */
1853 isl_qpolynomial_free(qp
);
1857 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_move(
1858 __isl_take isl_pw_qpolynomial
*pwqp
,
1859 enum isl_dim_type dst_type
, unsigned dst_pos
,
1860 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
1867 pwqp
->dim
= isl_dim_move(pwqp
->dim
,
1868 dst_type
, dst_pos
, src_type
, src_pos
, n
);
1872 for (i
= 0; i
< pwqp
->n
; ++i
) {
1873 pwqp
->p
[i
].set
= isl_set_move(pwqp
->p
[i
].set
, dst_type
, dst_pos
,
1874 src_type
, src_pos
, n
);
1875 if (!pwqp
->p
[i
].set
)
1877 pwqp
->p
[i
].qp
= isl_qpolynomial_move(pwqp
->p
[i
].qp
,
1878 dst_type
, dst_pos
, src_type
, src_pos
, n
);
1885 isl_pw_qpolynomial_free(pwqp
);
1889 __isl_give isl_dim
*isl_pw_qpolynomial_get_dim(
1890 __isl_keep isl_pw_qpolynomial
*pwqp
)
1895 return isl_dim_copy(pwqp
->dim
);
1898 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial
*pwqp
,
1899 enum isl_dim_type type
)
1901 return pwqp
? isl_dim_size(pwqp
->dim
, type
) : 0;
1904 __isl_give isl_term
*isl_term_alloc(__isl_take isl_dim
*dim
,
1905 __isl_take isl_mat
*div
)
1913 n
= isl_dim_total(dim
) + div
->n_row
;
1915 term
= isl_calloc(dim
->ctx
, struct isl_term
,
1916 sizeof(struct isl_term
) + (n
- 1) * sizeof(int));
1923 isl_int_init(term
->n
);
1924 isl_int_init(term
->d
);
1933 __isl_give isl_term
*isl_term_copy(__isl_keep isl_term
*term
)
1942 __isl_give isl_term
*isl_term_dup(__isl_keep isl_term
*term
)
1951 total
= isl_dim_total(term
->dim
) + term
->div
->n_row
;
1953 dup
= isl_term_alloc(isl_dim_copy(term
->dim
), isl_mat_copy(term
->div
));
1957 isl_int_set(dup
->n
, term
->n
);
1958 isl_int_set(dup
->d
, term
->d
);
1960 for (i
= 0; i
< total
; ++i
)
1961 dup
->pow
[i
] = term
->pow
[i
];
1966 __isl_give isl_term
*isl_term_cow(__isl_take isl_term
*term
)
1974 return isl_term_dup(term
);
1977 void isl_term_free(__isl_take isl_term
*term
)
1982 if (--term
->ref
> 0)
1985 isl_dim_free(term
->dim
);
1986 isl_mat_free(term
->div
);
1987 isl_int_clear(term
->n
);
1988 isl_int_clear(term
->d
);
1992 unsigned isl_term_dim(__isl_keep isl_term
*term
, enum isl_dim_type type
)
2000 case isl_dim_out
: return isl_dim_size(term
->dim
, type
);
2001 case isl_dim_div
: return term
->div
->n_row
;
2002 case isl_dim_all
: return isl_dim_total(term
->dim
) + term
->div
->n_row
;
2006 isl_ctx
*isl_term_get_ctx(__isl_keep isl_term
*term
)
2008 return term
? term
->dim
->ctx
: NULL
;
2011 void isl_term_get_num(__isl_keep isl_term
*term
, isl_int
*n
)
2015 isl_int_set(*n
, term
->n
);
2018 void isl_term_get_den(__isl_keep isl_term
*term
, isl_int
*d
)
2022 isl_int_set(*d
, term
->d
);
2025 int isl_term_get_exp(__isl_keep isl_term
*term
,
2026 enum isl_dim_type type
, unsigned pos
)
2031 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, type
), return -1);
2033 if (type
>= isl_dim_set
)
2034 pos
+= isl_dim_size(term
->dim
, isl_dim_param
);
2035 if (type
>= isl_dim_div
)
2036 pos
+= isl_dim_size(term
->dim
, isl_dim_set
);
2038 return term
->pow
[pos
];
2041 __isl_give isl_div
*isl_term_get_div(__isl_keep isl_term
*term
, unsigned pos
)
2043 isl_basic_map
*bmap
;
2050 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, isl_dim_div
),
2053 total
= term
->div
->n_col
- term
->div
->n_row
- 2;
2054 /* No nested divs for now */
2055 isl_assert(term
->dim
->ctx
,
2056 isl_seq_first_non_zero(term
->div
->row
[pos
] + 2 + total
,
2057 term
->div
->n_row
) == -1,
2060 bmap
= isl_basic_map_alloc_dim(isl_dim_copy(term
->dim
), 1, 0, 0);
2061 if ((k
= isl_basic_map_alloc_div(bmap
)) < 0)
2064 isl_seq_cpy(bmap
->div
[k
], term
->div
->row
[pos
], 2 + total
);
2066 return isl_basic_map_div(bmap
, k
);
2068 isl_basic_map_free(bmap
);
2072 __isl_give isl_term
*isl_upoly_foreach_term(__isl_keep
struct isl_upoly
*up
,
2073 int (*fn
)(__isl_take isl_term
*term
, void *user
),
2074 __isl_take isl_term
*term
, void *user
)
2077 struct isl_upoly_rec
*rec
;
2082 if (isl_upoly_is_zero(up
))
2085 isl_assert(up
->ctx
, !isl_upoly_is_nan(up
), goto error
);
2086 isl_assert(up
->ctx
, !isl_upoly_is_infty(up
), goto error
);
2087 isl_assert(up
->ctx
, !isl_upoly_is_neginfty(up
), goto error
);
2089 if (isl_upoly_is_cst(up
)) {
2090 struct isl_upoly_cst
*cst
;
2091 cst
= isl_upoly_as_cst(up
);
2094 term
= isl_term_cow(term
);
2097 isl_int_set(term
->n
, cst
->n
);
2098 isl_int_set(term
->d
, cst
->d
);
2099 if (fn(isl_term_copy(term
), user
) < 0)
2104 rec
= isl_upoly_as_rec(up
);
2108 for (i
= 0; i
< rec
->n
; ++i
) {
2109 term
= isl_term_cow(term
);
2112 term
->pow
[up
->var
] = i
;
2113 term
= isl_upoly_foreach_term(rec
->p
[i
], fn
, term
, user
);
2117 term
->pow
[up
->var
] = 0;
2121 isl_term_free(term
);
2125 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial
*qp
,
2126 int (*fn
)(__isl_take isl_term
*term
, void *user
), void *user
)
2133 term
= isl_term_alloc(isl_dim_copy(qp
->dim
), isl_mat_copy(qp
->div
));
2137 term
= isl_upoly_foreach_term(qp
->upoly
, fn
, term
, user
);
2139 isl_term_free(term
);
2141 return term
? 0 : -1;
2144 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial
*pwqp
,
2145 int (*fn
)(__isl_take isl_set
*set
, __isl_take isl_qpolynomial
*qp
,
2146 void *user
), void *user
)
2153 for (i
= 0; i
< pwqp
->n
; ++i
)
2154 if (fn(isl_set_copy(pwqp
->p
[i
].set
),
2155 isl_qpolynomial_copy(pwqp
->p
[i
].qp
), user
) < 0)
2161 static int any_divs(__isl_keep isl_set
*set
)
2168 for (i
= 0; i
< set
->n
; ++i
)
2169 if (set
->p
[i
]->n_div
> 0)
2175 __isl_give isl_qpolynomial
*isl_qpolynomial_lift(__isl_take isl_qpolynomial
*qp
,
2176 __isl_take isl_dim
*dim
)
2181 if (isl_dim_equal(qp
->dim
, dim
)) {
2186 qp
= isl_qpolynomial_cow(qp
);
2190 if (qp
->div
->n_row
) {
2196 extra
= isl_dim_size(dim
, isl_dim_set
) -
2197 isl_dim_size(qp
->dim
, isl_dim_set
);
2198 total
= isl_dim_total(qp
->dim
);
2199 exp
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
2202 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
2204 qp
->upoly
= expand(qp
->upoly
, exp
, total
);
2208 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + total
, extra
);
2211 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
2212 isl_seq_clr(qp
->div
->row
[i
] + 2 + total
, extra
);
2215 isl_dim_free(qp
->dim
);
2221 isl_qpolynomial_free(qp
);
2225 static int foreach_lifted_subset(__isl_take isl_set
*set
,
2226 __isl_take isl_qpolynomial
*qp
,
2227 int (*fn
)(__isl_take isl_set
*set
, __isl_take isl_qpolynomial
*qp
,
2228 void *user
), void *user
)
2235 for (i
= 0; i
< set
->n
; ++i
) {
2237 isl_qpolynomial
*copy
;
2239 lift
= isl_set_from_basic_set(isl_basic_set_copy(set
->p
[i
]));
2240 lift
= isl_set_lift(lift
);
2242 copy
= isl_qpolynomial_copy(qp
);
2243 copy
= isl_qpolynomial_lift(copy
, isl_set_get_dim(lift
));
2245 if (fn(lift
, copy
, user
) < 0)
2250 isl_qpolynomial_free(qp
);
2255 isl_qpolynomial_free(qp
);
2259 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial
*pwqp
,
2260 int (*fn
)(__isl_take isl_set
*set
, __isl_take isl_qpolynomial
*qp
,
2261 void *user
), void *user
)
2268 for (i
= 0; i
< pwqp
->n
; ++i
) {
2270 isl_qpolynomial
*qp
;
2272 set
= isl_set_copy(pwqp
->p
[i
].set
);
2273 qp
= isl_qpolynomial_copy(pwqp
->p
[i
].qp
);
2274 if (!any_divs(set
)) {
2275 if (fn(set
, qp
, user
) < 0)
2279 if (foreach_lifted_subset(set
, qp
, fn
, user
) < 0)
2286 static __isl_give isl_qpolynomial_fold
*qpolynomial_fold_alloc(
2287 enum isl_fold type
, __isl_take isl_dim
*dim
, int n
)
2289 isl_qpolynomial_fold
*fold
;
2294 isl_assert(dim
->ctx
, n
>= 0, goto error
);
2295 fold
= isl_alloc(dim
->ctx
, struct isl_qpolynomial_fold
,
2296 sizeof(struct isl_qpolynomial_fold
) +
2297 (n
- 1) * sizeof(struct isl_qpolynomial
*));
2313 __isl_give isl_qpolynomial_fold
*isl_qpolynomial_fold_empty(enum isl_fold type
,
2314 __isl_take isl_dim
*dim
)
2316 return qpolynomial_fold_alloc(type
, dim
, 0);
2319 __isl_give isl_qpolynomial_fold
*isl_qpolynomial_fold_alloc(
2320 enum isl_fold type
, __isl_take isl_qpolynomial
*qp
)
2322 isl_qpolynomial_fold
*fold
;
2327 fold
= qpolynomial_fold_alloc(type
, isl_dim_copy(qp
->dim
), 1);
2336 isl_qpolynomial_fold_free(fold
);
2337 isl_qpolynomial_free(qp
);
2341 __isl_give isl_qpolynomial_fold
*isl_qpolynomial_fold_copy(
2342 __isl_keep isl_qpolynomial_fold
*fold
)
2351 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold
*fold
)
2357 if (--fold
->ref
> 0)
2360 for (i
= 0; i
< fold
->n
; ++i
)
2361 isl_qpolynomial_free(fold
->qp
[i
]);
2362 isl_dim_free(fold
->dim
);
2366 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold
*fold
)
2371 return fold
->n
== 0;
2374 __isl_give isl_qpolynomial_fold
*isl_qpolynomial_fold_fold(
2375 __isl_take isl_qpolynomial_fold
*fold1
,
2376 __isl_take isl_qpolynomial_fold
*fold2
)
2379 struct isl_qpolynomial_fold
*res
= NULL
;
2381 if (!fold1
|| !fold2
)
2384 isl_assert(fold1
->dim
->ctx
, fold1
->type
== fold2
->type
, goto error
);
2385 isl_assert(fold1
->dim
->ctx
, isl_dim_equal(fold1
->dim
, fold2
->dim
),
2388 if (isl_qpolynomial_fold_is_empty(fold1
)) {
2389 isl_qpolynomial_fold_free(fold1
);
2393 if (isl_qpolynomial_fold_is_empty(fold2
)) {
2394 isl_qpolynomial_fold_free(fold2
);
2398 res
= qpolynomial_fold_alloc(fold1
->type
, isl_dim_copy(fold1
->dim
),
2399 fold1
->n
+ fold2
->n
);
2403 for (i
= 0; i
< fold1
->n
; ++i
) {
2404 res
->qp
[res
->n
] = isl_qpolynomial_copy(fold1
->qp
[i
]);
2405 if (!res
->qp
[res
->n
])
2410 for (i
= 0; i
< fold2
->n
; ++i
) {
2411 res
->qp
[res
->n
] = isl_qpolynomial_copy(fold2
->qp
[i
]);
2412 if (!res
->qp
[res
->n
])
2417 isl_qpolynomial_fold_free(fold1
);
2418 isl_qpolynomial_fold_free(fold2
);
2422 isl_qpolynomial_fold_free(res
);
2423 isl_qpolynomial_fold_free(fold1
);
2424 isl_qpolynomial_fold_free(fold2
);
2428 __isl_give isl_pw_qpolynomial_fold
*isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2429 enum isl_fold type
, __isl_take isl_pw_qpolynomial
*pwqp
)
2432 isl_pw_qpolynomial_fold
*pwf
;
2437 pwf
= isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp
->dim
), pwqp
->n
);
2439 for (i
= 0; i
< pwqp
->n
; ++i
)
2440 pwf
= isl_pw_qpolynomial_fold_add_piece(pwf
,
2441 isl_set_copy(pwqp
->p
[i
].set
),
2442 isl_qpolynomial_fold_alloc(type
,
2443 isl_qpolynomial_copy(pwqp
->p
[i
].qp
)));
2445 isl_pw_qpolynomial_free(pwqp
);
2450 /* For each parameter or variable that does not appear in qp,
2451 * first eliminate the variable from all constraints and then set it to zero.
2453 static __isl_give isl_set
*fix_inactive(__isl_take isl_set
*set
,
2454 __isl_keep isl_qpolynomial
*qp
)
2465 d
= isl_dim_total(set
->dim
);
2466 active
= isl_calloc_array(set
->ctx
, int, d
);
2467 if (set_active(qp
, active
) < 0)
2470 for (i
= 0; i
< d
; ++i
)
2479 nparam
= isl_dim_size(set
->dim
, isl_dim_param
);
2480 nvar
= isl_dim_size(set
->dim
, isl_dim_set
);
2481 for (i
= 0; i
< nparam
; ++i
) {
2484 set
= isl_set_eliminate(set
, isl_dim_param
, i
, 1);
2485 set
= isl_set_fix_si(set
, isl_dim_param
, i
, 0);
2487 for (i
= 0; i
< nvar
; ++i
) {
2490 set
= isl_set_eliminate(set
, isl_dim_set
, i
, 1);
2491 set
= isl_set_fix_si(set
, isl_dim_set
, i
, 0);
2503 struct isl_max_data
{
2504 isl_qpolynomial
*qp
;
2506 isl_qpolynomial
*max
;
2509 static int max_fn(__isl_take isl_point
*pnt
, void *user
)
2511 struct isl_max_data
*data
= (struct isl_max_data
*)user
;
2512 isl_qpolynomial
*val
;
2514 val
= isl_qpolynomial_eval(isl_qpolynomial_copy(data
->qp
), pnt
);
2519 data
->max
= qpolynomial_max(data
->max
, val
);
2525 static __isl_give isl_qpolynomial
*guarded_qpolynomial_max(
2526 __isl_take isl_set
*set
, __isl_take isl_qpolynomial
*qp
)
2528 struct isl_max_data data
= { NULL
, 1, NULL
};
2533 if (isl_upoly_is_cst(qp
->upoly
)) {
2538 set
= fix_inactive(set
, qp
);
2541 if (isl_set_foreach_point(set
, max_fn
, &data
) < 0)
2545 isl_qpolynomial_free(qp
);
2549 isl_qpolynomial_free(qp
);
2550 isl_qpolynomial_free(data
.max
);
2554 /* Compute the maximal value attained by the piecewise quasipolynomial
2555 * on its domain or zero if the domain is empty.
2556 * In the worst case, the domain is scanned completely,
2557 * so the domain is assumed to be bounded.
2559 __isl_give isl_qpolynomial
*isl_pw_qpolynomial_max(
2560 __isl_take isl_pw_qpolynomial
*pwqp
)
2563 isl_qpolynomial
*max
;
2569 isl_dim
*dim
= isl_dim_copy(pwqp
->dim
);
2570 isl_pw_qpolynomial_free(pwqp
);
2571 return isl_qpolynomial_zero(dim
);
2574 max
= guarded_qpolynomial_max(isl_set_copy(pwqp
->p
[0].set
),
2575 isl_qpolynomial_copy(pwqp
->p
[0].qp
));
2576 for (i
= 1; i
< pwqp
->n
; ++i
) {
2577 isl_qpolynomial
*max_i
;
2578 max_i
= guarded_qpolynomial_max(isl_set_copy(pwqp
->p
[i
].set
),
2579 isl_qpolynomial_copy(pwqp
->p
[i
].qp
));
2580 max
= qpolynomial_max(max
, max_i
);
2583 isl_pw_qpolynomial_free(pwqp
);