2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the GNU LGPLv2.1 license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
13 #include <isl_union_map_private.h>
14 #include <isl_polynomial_private.h>
15 #include <isl_point_private.h>
16 #include <isl_dim_private.h>
17 #include <isl_map_private.h>
19 static unsigned pos(__isl_keep isl_dim
*dim
, enum isl_dim_type type
)
22 case isl_dim_param
: return 0;
23 case isl_dim_in
: return dim
->nparam
;
24 case isl_dim_out
: return dim
->nparam
+ dim
->n_in
;
29 int isl_upoly_is_cst(__isl_keep
struct isl_upoly
*up
)
37 __isl_keep
struct isl_upoly_cst
*isl_upoly_as_cst(__isl_keep
struct isl_upoly
*up
)
42 isl_assert(up
->ctx
, up
->var
< 0, return NULL
);
44 return (struct isl_upoly_cst
*)up
;
47 __isl_keep
struct isl_upoly_rec
*isl_upoly_as_rec(__isl_keep
struct isl_upoly
*up
)
52 isl_assert(up
->ctx
, up
->var
>= 0, return NULL
);
54 return (struct isl_upoly_rec
*)up
;
57 int isl_upoly_is_equal(__isl_keep
struct isl_upoly
*up1
,
58 __isl_keep
struct isl_upoly
*up2
)
61 struct isl_upoly_rec
*rec1
, *rec2
;
67 if (up1
->var
!= up2
->var
)
69 if (isl_upoly_is_cst(up1
)) {
70 struct isl_upoly_cst
*cst1
, *cst2
;
71 cst1
= isl_upoly_as_cst(up1
);
72 cst2
= isl_upoly_as_cst(up2
);
75 return isl_int_eq(cst1
->n
, cst2
->n
) &&
76 isl_int_eq(cst1
->d
, cst2
->d
);
79 rec1
= isl_upoly_as_rec(up1
);
80 rec2
= isl_upoly_as_rec(up2
);
84 if (rec1
->n
!= rec2
->n
)
87 for (i
= 0; i
< rec1
->n
; ++i
) {
88 int eq
= isl_upoly_is_equal(rec1
->p
[i
], rec2
->p
[i
]);
96 int isl_upoly_is_zero(__isl_keep
struct isl_upoly
*up
)
98 struct isl_upoly_cst
*cst
;
102 if (!isl_upoly_is_cst(up
))
105 cst
= isl_upoly_as_cst(up
);
109 return isl_int_is_zero(cst
->n
) && isl_int_is_pos(cst
->d
);
112 int isl_upoly_sgn(__isl_keep
struct isl_upoly
*up
)
114 struct isl_upoly_cst
*cst
;
118 if (!isl_upoly_is_cst(up
))
121 cst
= isl_upoly_as_cst(up
);
125 return isl_int_sgn(cst
->n
);
128 int isl_upoly_is_nan(__isl_keep
struct isl_upoly
*up
)
130 struct isl_upoly_cst
*cst
;
134 if (!isl_upoly_is_cst(up
))
137 cst
= isl_upoly_as_cst(up
);
141 return isl_int_is_zero(cst
->n
) && isl_int_is_zero(cst
->d
);
144 int isl_upoly_is_infty(__isl_keep
struct isl_upoly
*up
)
146 struct isl_upoly_cst
*cst
;
150 if (!isl_upoly_is_cst(up
))
153 cst
= isl_upoly_as_cst(up
);
157 return isl_int_is_pos(cst
->n
) && isl_int_is_zero(cst
->d
);
160 int isl_upoly_is_neginfty(__isl_keep
struct isl_upoly
*up
)
162 struct isl_upoly_cst
*cst
;
166 if (!isl_upoly_is_cst(up
))
169 cst
= isl_upoly_as_cst(up
);
173 return isl_int_is_neg(cst
->n
) && isl_int_is_zero(cst
->d
);
176 int isl_upoly_is_one(__isl_keep
struct isl_upoly
*up
)
178 struct isl_upoly_cst
*cst
;
182 if (!isl_upoly_is_cst(up
))
185 cst
= isl_upoly_as_cst(up
);
189 return isl_int_eq(cst
->n
, cst
->d
) && isl_int_is_pos(cst
->d
);
192 int isl_upoly_is_negone(__isl_keep
struct isl_upoly
*up
)
194 struct isl_upoly_cst
*cst
;
198 if (!isl_upoly_is_cst(up
))
201 cst
= isl_upoly_as_cst(up
);
205 return isl_int_is_negone(cst
->n
) && isl_int_is_one(cst
->d
);
208 __isl_give
struct isl_upoly_cst
*isl_upoly_cst_alloc(struct isl_ctx
*ctx
)
210 struct isl_upoly_cst
*cst
;
212 cst
= isl_alloc_type(ctx
, struct isl_upoly_cst
);
221 isl_int_init(cst
->n
);
222 isl_int_init(cst
->d
);
227 __isl_give
struct isl_upoly
*isl_upoly_zero(struct isl_ctx
*ctx
)
229 struct isl_upoly_cst
*cst
;
231 cst
= isl_upoly_cst_alloc(ctx
);
235 isl_int_set_si(cst
->n
, 0);
236 isl_int_set_si(cst
->d
, 1);
241 __isl_give
struct isl_upoly
*isl_upoly_infty(struct isl_ctx
*ctx
)
243 struct isl_upoly_cst
*cst
;
245 cst
= isl_upoly_cst_alloc(ctx
);
249 isl_int_set_si(cst
->n
, 1);
250 isl_int_set_si(cst
->d
, 0);
255 __isl_give
struct isl_upoly
*isl_upoly_neginfty(struct isl_ctx
*ctx
)
257 struct isl_upoly_cst
*cst
;
259 cst
= isl_upoly_cst_alloc(ctx
);
263 isl_int_set_si(cst
->n
, -1);
264 isl_int_set_si(cst
->d
, 0);
269 __isl_give
struct isl_upoly
*isl_upoly_nan(struct isl_ctx
*ctx
)
271 struct isl_upoly_cst
*cst
;
273 cst
= isl_upoly_cst_alloc(ctx
);
277 isl_int_set_si(cst
->n
, 0);
278 isl_int_set_si(cst
->d
, 0);
283 __isl_give
struct isl_upoly
*isl_upoly_rat_cst(struct isl_ctx
*ctx
,
284 isl_int n
, isl_int d
)
286 struct isl_upoly_cst
*cst
;
288 cst
= isl_upoly_cst_alloc(ctx
);
292 isl_int_set(cst
->n
, n
);
293 isl_int_set(cst
->d
, d
);
298 __isl_give
struct isl_upoly_rec
*isl_upoly_alloc_rec(struct isl_ctx
*ctx
,
301 struct isl_upoly_rec
*rec
;
303 isl_assert(ctx
, var
>= 0, return NULL
);
304 isl_assert(ctx
, size
>= 0, return NULL
);
305 rec
= isl_calloc(ctx
, struct isl_upoly_rec
,
306 sizeof(struct isl_upoly_rec
) +
307 (size
- 1) * sizeof(struct isl_upoly
*));
322 isl_ctx
*isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial
*qp
)
324 return qp
? qp
->dim
->ctx
: NULL
;
327 __isl_give isl_dim
*isl_qpolynomial_get_dim(__isl_keep isl_qpolynomial
*qp
)
329 return qp
? isl_dim_copy(qp
->dim
) : NULL
;
332 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial
*qp
,
333 enum isl_dim_type type
)
335 return qp
? isl_dim_size(qp
->dim
, type
) : 0;
338 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial
*qp
)
340 return qp
? isl_upoly_is_zero(qp
->upoly
) : -1;
343 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial
*qp
)
345 return qp
? isl_upoly_is_one(qp
->upoly
) : -1;
348 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial
*qp
)
350 return qp
? isl_upoly_is_nan(qp
->upoly
) : -1;
353 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial
*qp
)
355 return qp
? isl_upoly_is_infty(qp
->upoly
) : -1;
358 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial
*qp
)
360 return qp
? isl_upoly_is_neginfty(qp
->upoly
) : -1;
363 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial
*qp
)
365 return qp
? isl_upoly_sgn(qp
->upoly
) : 0;
368 static void upoly_free_cst(__isl_take
struct isl_upoly_cst
*cst
)
370 isl_int_clear(cst
->n
);
371 isl_int_clear(cst
->d
);
374 static void upoly_free_rec(__isl_take
struct isl_upoly_rec
*rec
)
378 for (i
= 0; i
< rec
->n
; ++i
)
379 isl_upoly_free(rec
->p
[i
]);
382 __isl_give
struct isl_upoly
*isl_upoly_copy(__isl_keep
struct isl_upoly
*up
)
391 __isl_give
struct isl_upoly
*isl_upoly_dup_cst(__isl_keep
struct isl_upoly
*up
)
393 struct isl_upoly_cst
*cst
;
394 struct isl_upoly_cst
*dup
;
396 cst
= isl_upoly_as_cst(up
);
400 dup
= isl_upoly_as_cst(isl_upoly_zero(up
->ctx
));
403 isl_int_set(dup
->n
, cst
->n
);
404 isl_int_set(dup
->d
, cst
->d
);
409 __isl_give
struct isl_upoly
*isl_upoly_dup_rec(__isl_keep
struct isl_upoly
*up
)
412 struct isl_upoly_rec
*rec
;
413 struct isl_upoly_rec
*dup
;
415 rec
= isl_upoly_as_rec(up
);
419 dup
= isl_upoly_alloc_rec(up
->ctx
, up
->var
, rec
->n
);
423 for (i
= 0; i
< rec
->n
; ++i
) {
424 dup
->p
[i
] = isl_upoly_copy(rec
->p
[i
]);
432 isl_upoly_free(&dup
->up
);
436 __isl_give
struct isl_upoly
*isl_upoly_dup(__isl_keep
struct isl_upoly
*up
)
438 struct isl_upoly
*dup
;
443 if (isl_upoly_is_cst(up
))
444 return isl_upoly_dup_cst(up
);
446 return isl_upoly_dup_rec(up
);
449 __isl_give
struct isl_upoly
*isl_upoly_cow(__isl_take
struct isl_upoly
*up
)
457 return isl_upoly_dup(up
);
460 void isl_upoly_free(__isl_take
struct isl_upoly
*up
)
469 upoly_free_cst((struct isl_upoly_cst
*)up
);
471 upoly_free_rec((struct isl_upoly_rec
*)up
);
473 isl_ctx_deref(up
->ctx
);
477 static void isl_upoly_cst_reduce(__isl_keep
struct isl_upoly_cst
*cst
)
482 isl_int_gcd(gcd
, cst
->n
, cst
->d
);
483 if (!isl_int_is_zero(gcd
) && !isl_int_is_one(gcd
)) {
484 isl_int_divexact(cst
->n
, cst
->n
, gcd
);
485 isl_int_divexact(cst
->d
, cst
->d
, gcd
);
490 __isl_give
struct isl_upoly
*isl_upoly_sum_cst(__isl_take
struct isl_upoly
*up1
,
491 __isl_take
struct isl_upoly
*up2
)
493 struct isl_upoly_cst
*cst1
;
494 struct isl_upoly_cst
*cst2
;
496 up1
= isl_upoly_cow(up1
);
500 cst1
= isl_upoly_as_cst(up1
);
501 cst2
= isl_upoly_as_cst(up2
);
503 if (isl_int_eq(cst1
->d
, cst2
->d
))
504 isl_int_add(cst1
->n
, cst1
->n
, cst2
->n
);
506 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->d
);
507 isl_int_addmul(cst1
->n
, cst2
->n
, cst1
->d
);
508 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
511 isl_upoly_cst_reduce(cst1
);
521 static __isl_give
struct isl_upoly
*replace_by_zero(
522 __isl_take
struct isl_upoly
*up
)
530 return isl_upoly_zero(ctx
);
533 static __isl_give
struct isl_upoly
*replace_by_constant_term(
534 __isl_take
struct isl_upoly
*up
)
536 struct isl_upoly_rec
*rec
;
537 struct isl_upoly
*cst
;
542 rec
= isl_upoly_as_rec(up
);
545 cst
= isl_upoly_copy(rec
->p
[0]);
553 __isl_give
struct isl_upoly
*isl_upoly_sum(__isl_take
struct isl_upoly
*up1
,
554 __isl_take
struct isl_upoly
*up2
)
557 struct isl_upoly_rec
*rec1
, *rec2
;
562 if (isl_upoly_is_nan(up1
)) {
567 if (isl_upoly_is_nan(up2
)) {
572 if (isl_upoly_is_zero(up1
)) {
577 if (isl_upoly_is_zero(up2
)) {
582 if (up1
->var
< up2
->var
)
583 return isl_upoly_sum(up2
, up1
);
585 if (up2
->var
< up1
->var
) {
586 struct isl_upoly_rec
*rec
;
587 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
591 up1
= isl_upoly_cow(up1
);
592 rec
= isl_upoly_as_rec(up1
);
595 rec
->p
[0] = isl_upoly_sum(rec
->p
[0], up2
);
597 up1
= replace_by_constant_term(up1
);
601 if (isl_upoly_is_cst(up1
))
602 return isl_upoly_sum_cst(up1
, up2
);
604 rec1
= isl_upoly_as_rec(up1
);
605 rec2
= isl_upoly_as_rec(up2
);
609 if (rec1
->n
< rec2
->n
)
610 return isl_upoly_sum(up2
, up1
);
612 up1
= isl_upoly_cow(up1
);
613 rec1
= isl_upoly_as_rec(up1
);
617 for (i
= rec2
->n
- 1; i
>= 0; --i
) {
618 rec1
->p
[i
] = isl_upoly_sum(rec1
->p
[i
],
619 isl_upoly_copy(rec2
->p
[i
]));
622 if (i
== rec1
->n
- 1 && isl_upoly_is_zero(rec1
->p
[i
])) {
623 isl_upoly_free(rec1
->p
[i
]);
629 up1
= replace_by_zero(up1
);
630 else if (rec1
->n
== 1)
631 up1
= replace_by_constant_term(up1
);
642 __isl_give
struct isl_upoly
*isl_upoly_neg_cst(__isl_take
struct isl_upoly
*up
)
644 struct isl_upoly_cst
*cst
;
646 if (isl_upoly_is_zero(up
))
649 up
= isl_upoly_cow(up
);
653 cst
= isl_upoly_as_cst(up
);
655 isl_int_neg(cst
->n
, cst
->n
);
660 __isl_give
struct isl_upoly
*isl_upoly_neg(__isl_take
struct isl_upoly
*up
)
663 struct isl_upoly_rec
*rec
;
668 if (isl_upoly_is_cst(up
))
669 return isl_upoly_neg_cst(up
);
671 up
= isl_upoly_cow(up
);
672 rec
= isl_upoly_as_rec(up
);
676 for (i
= 0; i
< rec
->n
; ++i
) {
677 rec
->p
[i
] = isl_upoly_neg(rec
->p
[i
]);
688 __isl_give
struct isl_upoly
*isl_upoly_mul_cst(__isl_take
struct isl_upoly
*up1
,
689 __isl_take
struct isl_upoly
*up2
)
691 struct isl_upoly_cst
*cst1
;
692 struct isl_upoly_cst
*cst2
;
694 up1
= isl_upoly_cow(up1
);
698 cst1
= isl_upoly_as_cst(up1
);
699 cst2
= isl_upoly_as_cst(up2
);
701 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->n
);
702 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
704 isl_upoly_cst_reduce(cst1
);
714 __isl_give
struct isl_upoly
*isl_upoly_mul_rec(__isl_take
struct isl_upoly
*up1
,
715 __isl_take
struct isl_upoly
*up2
)
717 struct isl_upoly_rec
*rec1
;
718 struct isl_upoly_rec
*rec2
;
719 struct isl_upoly_rec
*res
;
723 rec1
= isl_upoly_as_rec(up1
);
724 rec2
= isl_upoly_as_rec(up2
);
727 size
= rec1
->n
+ rec2
->n
- 1;
728 res
= isl_upoly_alloc_rec(up1
->ctx
, up1
->var
, size
);
732 for (i
= 0; i
< rec1
->n
; ++i
) {
733 res
->p
[i
] = isl_upoly_mul(isl_upoly_copy(rec2
->p
[0]),
734 isl_upoly_copy(rec1
->p
[i
]));
739 for (; i
< size
; ++i
) {
740 res
->p
[i
] = isl_upoly_zero(up1
->ctx
);
745 for (i
= 0; i
< rec1
->n
; ++i
) {
746 for (j
= 1; j
< rec2
->n
; ++j
) {
747 struct isl_upoly
*up
;
748 up
= isl_upoly_mul(isl_upoly_copy(rec2
->p
[j
]),
749 isl_upoly_copy(rec1
->p
[i
]));
750 res
->p
[i
+ j
] = isl_upoly_sum(res
->p
[i
+ j
], up
);
763 isl_upoly_free(&res
->up
);
767 __isl_give
struct isl_upoly
*isl_upoly_mul(__isl_take
struct isl_upoly
*up1
,
768 __isl_take
struct isl_upoly
*up2
)
773 if (isl_upoly_is_nan(up1
)) {
778 if (isl_upoly_is_nan(up2
)) {
783 if (isl_upoly_is_zero(up1
)) {
788 if (isl_upoly_is_zero(up2
)) {
793 if (isl_upoly_is_one(up1
)) {
798 if (isl_upoly_is_one(up2
)) {
803 if (up1
->var
< up2
->var
)
804 return isl_upoly_mul(up2
, up1
);
806 if (up2
->var
< up1
->var
) {
808 struct isl_upoly_rec
*rec
;
809 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
810 isl_ctx
*ctx
= up1
->ctx
;
813 return isl_upoly_nan(ctx
);
815 up1
= isl_upoly_cow(up1
);
816 rec
= isl_upoly_as_rec(up1
);
820 for (i
= 0; i
< rec
->n
; ++i
) {
821 rec
->p
[i
] = isl_upoly_mul(rec
->p
[i
],
822 isl_upoly_copy(up2
));
830 if (isl_upoly_is_cst(up1
))
831 return isl_upoly_mul_cst(up1
, up2
);
833 return isl_upoly_mul_rec(up1
, up2
);
840 __isl_give isl_qpolynomial
*isl_qpolynomial_alloc(__isl_take isl_dim
*dim
,
841 unsigned n_div
, __isl_take
struct isl_upoly
*up
)
843 struct isl_qpolynomial
*qp
= NULL
;
849 total
= isl_dim_total(dim
);
851 qp
= isl_calloc_type(dim
->ctx
, struct isl_qpolynomial
);
856 qp
->div
= isl_mat_alloc(dim
->ctx
, n_div
, 1 + 1 + total
+ n_div
);
867 isl_qpolynomial_free(qp
);
871 __isl_give isl_qpolynomial
*isl_qpolynomial_copy(__isl_keep isl_qpolynomial
*qp
)
880 __isl_give isl_qpolynomial
*isl_qpolynomial_dup(__isl_keep isl_qpolynomial
*qp
)
882 struct isl_qpolynomial
*dup
;
887 dup
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
,
888 isl_upoly_copy(qp
->upoly
));
891 isl_mat_free(dup
->div
);
892 dup
->div
= isl_mat_copy(qp
->div
);
898 isl_qpolynomial_free(dup
);
902 __isl_give isl_qpolynomial
*isl_qpolynomial_cow(__isl_take isl_qpolynomial
*qp
)
910 return isl_qpolynomial_dup(qp
);
913 void isl_qpolynomial_free(__isl_take isl_qpolynomial
*qp
)
921 isl_dim_free(qp
->dim
);
922 isl_mat_free(qp
->div
);
923 isl_upoly_free(qp
->upoly
);
928 static int compatible_divs(__isl_keep isl_mat
*div1
, __isl_keep isl_mat
*div2
)
933 isl_assert(div1
->ctx
, div1
->n_row
>= div2
->n_row
&&
934 div1
->n_col
>= div2
->n_col
, return -1);
936 if (div1
->n_row
== div2
->n_row
)
937 return isl_mat_is_equal(div1
, div2
);
941 div1
->n_row
= div2
->n_row
;
942 div1
->n_col
= div2
->n_col
;
944 equal
= isl_mat_is_equal(div1
, div2
);
952 static void expand_row(__isl_keep isl_mat
*dst
, int d
,
953 __isl_keep isl_mat
*src
, int s
, int *exp
)
956 unsigned c
= src
->n_col
- src
->n_row
;
958 isl_seq_cpy(dst
->row
[d
], src
->row
[s
], c
);
959 isl_seq_clr(dst
->row
[d
] + c
, dst
->n_col
- c
);
961 for (i
= 0; i
< s
; ++i
)
962 isl_int_set(dst
->row
[d
][c
+ exp
[i
]], src
->row
[s
][c
+ i
]);
965 static int cmp_row(__isl_keep isl_mat
*div
, int i
, int j
)
969 li
= isl_seq_last_non_zero(div
->row
[i
], div
->n_col
);
970 lj
= isl_seq_last_non_zero(div
->row
[j
], div
->n_col
);
975 return isl_seq_cmp(div
->row
[i
], div
->row
[j
], div
->n_col
);
978 struct isl_div_sort_info
{
983 static int div_sort_cmp(const void *p1
, const void *p2
)
985 const struct isl_div_sort_info
*i1
, *i2
;
986 i1
= (const struct isl_div_sort_info
*) p1
;
987 i2
= (const struct isl_div_sort_info
*) p2
;
989 return cmp_row(i1
->div
, i1
->row
, i2
->row
);
992 static __isl_give isl_mat
*sort_divs(__isl_take isl_mat
*div
)
995 struct isl_div_sort_info
*array
= NULL
;
1000 if (div
->n_row
<= 1)
1003 array
= isl_alloc_array(div
->ctx
, struct isl_div_sort_info
, div
->n_row
);
1004 pos
= isl_alloc_array(div
->ctx
, int, div
->n_row
);
1008 for (i
= 0; i
< div
->n_row
; ++i
) {
1014 qsort(array
, div
->n_row
, sizeof(struct isl_div_sort_info
),
1017 for (i
= 0; i
< div
->n_row
; ++i
) {
1019 if (pos
[array
[i
].row
] == i
)
1021 div
= isl_mat_cow(div
);
1022 div
= isl_mat_swap_rows(div
, i
, pos
[array
[i
].row
]);
1023 t
= pos
[array
[i
].row
];
1024 pos
[array
[i
].row
] = pos
[i
];
1038 static __isl_give isl_mat
*merge_divs(__isl_keep isl_mat
*div1
,
1039 __isl_keep isl_mat
*div2
, int *exp1
, int *exp2
)
1042 isl_mat
*div
= NULL
;
1043 unsigned d
= div1
->n_col
- div1
->n_row
;
1045 div
= isl_mat_alloc(div1
->ctx
, 1 + div1
->n_row
+ div2
->n_row
,
1046 d
+ div1
->n_row
+ div2
->n_row
);
1050 for (i
= 0, j
= 0, k
= 0; i
< div1
->n_row
&& j
< div2
->n_row
; ++k
) {
1053 expand_row(div
, k
, div1
, i
, exp1
);
1054 expand_row(div
, k
+ 1, div2
, j
, exp2
);
1056 cmp
= cmp_row(div
, k
, k
+ 1);
1060 } else if (cmp
< 0) {
1064 isl_seq_cpy(div
->row
[k
], div
->row
[k
+ 1], div
->n_col
);
1067 for (; i
< div1
->n_row
; ++i
, ++k
) {
1068 expand_row(div
, k
, div1
, i
, exp1
);
1071 for (; j
< div2
->n_row
; ++j
, ++k
) {
1072 expand_row(div
, k
, div2
, j
, exp2
);
1082 static __isl_give
struct isl_upoly
*expand(__isl_take
struct isl_upoly
*up
,
1083 int *exp
, int first
)
1086 struct isl_upoly_rec
*rec
;
1088 if (isl_upoly_is_cst(up
))
1091 if (up
->var
< first
)
1094 if (exp
[up
->var
- first
] == up
->var
- first
)
1097 up
= isl_upoly_cow(up
);
1101 up
->var
= exp
[up
->var
- first
] + first
;
1103 rec
= isl_upoly_as_rec(up
);
1107 for (i
= 0; i
< rec
->n
; ++i
) {
1108 rec
->p
[i
] = expand(rec
->p
[i
], exp
, first
);
1119 static __isl_give isl_qpolynomial
*with_merged_divs(
1120 __isl_give isl_qpolynomial
*(*fn
)(__isl_take isl_qpolynomial
*qp1
,
1121 __isl_take isl_qpolynomial
*qp2
),
1122 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1126 isl_mat
*div
= NULL
;
1128 qp1
= isl_qpolynomial_cow(qp1
);
1129 qp2
= isl_qpolynomial_cow(qp2
);
1134 isl_assert(qp1
->div
->ctx
, qp1
->div
->n_row
>= qp2
->div
->n_row
&&
1135 qp1
->div
->n_col
>= qp2
->div
->n_col
, goto error
);
1137 exp1
= isl_alloc_array(qp1
->div
->ctx
, int, qp1
->div
->n_row
);
1138 exp2
= isl_alloc_array(qp2
->div
->ctx
, int, qp2
->div
->n_row
);
1142 div
= merge_divs(qp1
->div
, qp2
->div
, exp1
, exp2
);
1146 isl_mat_free(qp1
->div
);
1147 qp1
->div
= isl_mat_copy(div
);
1148 isl_mat_free(qp2
->div
);
1149 qp2
->div
= isl_mat_copy(div
);
1151 qp1
->upoly
= expand(qp1
->upoly
, exp1
, div
->n_col
- div
->n_row
- 2);
1152 qp2
->upoly
= expand(qp2
->upoly
, exp2
, div
->n_col
- div
->n_row
- 2);
1154 if (!qp1
->upoly
|| !qp2
->upoly
)
1161 return fn(qp1
, qp2
);
1166 isl_qpolynomial_free(qp1
);
1167 isl_qpolynomial_free(qp2
);
1171 __isl_give isl_qpolynomial
*isl_qpolynomial_add(__isl_take isl_qpolynomial
*qp1
,
1172 __isl_take isl_qpolynomial
*qp2
)
1174 qp1
= isl_qpolynomial_cow(qp1
);
1179 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1180 return isl_qpolynomial_add(qp2
, qp1
);
1182 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1183 if (!compatible_divs(qp1
->div
, qp2
->div
))
1184 return with_merged_divs(isl_qpolynomial_add
, qp1
, qp2
);
1186 qp1
->upoly
= isl_upoly_sum(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1190 isl_qpolynomial_free(qp2
);
1194 isl_qpolynomial_free(qp1
);
1195 isl_qpolynomial_free(qp2
);
1199 __isl_give isl_qpolynomial
*isl_qpolynomial_sub(__isl_take isl_qpolynomial
*qp1
,
1200 __isl_take isl_qpolynomial
*qp2
)
1202 return isl_qpolynomial_add(qp1
, isl_qpolynomial_neg(qp2
));
1205 __isl_give isl_qpolynomial
*isl_qpolynomial_neg(__isl_take isl_qpolynomial
*qp
)
1207 qp
= isl_qpolynomial_cow(qp
);
1212 qp
->upoly
= isl_upoly_neg(qp
->upoly
);
1218 isl_qpolynomial_free(qp
);
1222 __isl_give isl_qpolynomial
*isl_qpolynomial_mul(__isl_take isl_qpolynomial
*qp1
,
1223 __isl_take isl_qpolynomial
*qp2
)
1225 qp1
= isl_qpolynomial_cow(qp1
);
1230 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1231 return isl_qpolynomial_mul(qp2
, qp1
);
1233 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1234 if (!compatible_divs(qp1
->div
, qp2
->div
))
1235 return with_merged_divs(isl_qpolynomial_mul
, qp1
, qp2
);
1237 qp1
->upoly
= isl_upoly_mul(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1241 isl_qpolynomial_free(qp2
);
1245 isl_qpolynomial_free(qp1
);
1246 isl_qpolynomial_free(qp2
);
1250 __isl_give isl_qpolynomial
*isl_qpolynomial_zero(__isl_take isl_dim
*dim
)
1252 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1255 __isl_give isl_qpolynomial
*isl_qpolynomial_infty(__isl_take isl_dim
*dim
)
1257 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_infty(dim
->ctx
));
1260 __isl_give isl_qpolynomial
*isl_qpolynomial_neginfty(__isl_take isl_dim
*dim
)
1262 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_neginfty(dim
->ctx
));
1265 __isl_give isl_qpolynomial
*isl_qpolynomial_nan(__isl_take isl_dim
*dim
)
1267 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_nan(dim
->ctx
));
1270 __isl_give isl_qpolynomial
*isl_qpolynomial_cst(__isl_take isl_dim
*dim
,
1273 struct isl_qpolynomial
*qp
;
1274 struct isl_upoly_cst
*cst
;
1276 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1280 cst
= isl_upoly_as_cst(qp
->upoly
);
1281 isl_int_set(cst
->n
, v
);
1286 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial
*qp
,
1287 isl_int
*n
, isl_int
*d
)
1289 struct isl_upoly_cst
*cst
;
1294 if (!isl_upoly_is_cst(qp
->upoly
))
1297 cst
= isl_upoly_as_cst(qp
->upoly
);
1302 isl_int_set(*n
, cst
->n
);
1304 isl_int_set(*d
, cst
->d
);
1309 int isl_upoly_is_affine(__isl_keep
struct isl_upoly
*up
)
1312 struct isl_upoly_rec
*rec
;
1320 rec
= isl_upoly_as_rec(up
);
1327 isl_assert(up
->ctx
, rec
->n
> 1, return -1);
1329 is_cst
= isl_upoly_is_cst(rec
->p
[1]);
1335 return isl_upoly_is_affine(rec
->p
[0]);
1338 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial
*qp
)
1343 if (qp
->div
->n_row
> 0)
1346 return isl_upoly_is_affine(qp
->upoly
);
1349 static void update_coeff(__isl_keep isl_vec
*aff
,
1350 __isl_keep
struct isl_upoly_cst
*cst
, int pos
)
1355 if (isl_int_is_zero(cst
->n
))
1360 isl_int_gcd(gcd
, cst
->d
, aff
->el
[0]);
1361 isl_int_divexact(f
, cst
->d
, gcd
);
1362 isl_int_divexact(gcd
, aff
->el
[0], gcd
);
1363 isl_seq_scale(aff
->el
, aff
->el
, f
, aff
->size
);
1364 isl_int_mul(aff
->el
[1 + pos
], gcd
, cst
->n
);
1369 int isl_upoly_update_affine(__isl_keep
struct isl_upoly
*up
,
1370 __isl_keep isl_vec
*aff
)
1372 struct isl_upoly_cst
*cst
;
1373 struct isl_upoly_rec
*rec
;
1379 struct isl_upoly_cst
*cst
;
1381 cst
= isl_upoly_as_cst(up
);
1384 update_coeff(aff
, cst
, 0);
1388 rec
= isl_upoly_as_rec(up
);
1391 isl_assert(up
->ctx
, rec
->n
== 2, return -1);
1393 cst
= isl_upoly_as_cst(rec
->p
[1]);
1396 update_coeff(aff
, cst
, 1 + up
->var
);
1398 return isl_upoly_update_affine(rec
->p
[0], aff
);
1401 __isl_give isl_vec
*isl_qpolynomial_extract_affine(
1402 __isl_keep isl_qpolynomial
*qp
)
1410 isl_assert(qp
->div
->ctx
, qp
->div
->n_row
== 0, return NULL
);
1411 d
= isl_dim_total(qp
->dim
);
1412 aff
= isl_vec_alloc(qp
->div
->ctx
, 2 + d
);
1416 isl_seq_clr(aff
->el
+ 1, 1 + d
);
1417 isl_int_set_si(aff
->el
[0], 1);
1419 if (isl_upoly_update_affine(qp
->upoly
, aff
) < 0)
1428 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial
*qp1
,
1429 __isl_keep isl_qpolynomial
*qp2
)
1434 return isl_upoly_is_equal(qp1
->upoly
, qp2
->upoly
);
1437 static void upoly_update_den(__isl_keep
struct isl_upoly
*up
, isl_int
*d
)
1440 struct isl_upoly_rec
*rec
;
1442 if (isl_upoly_is_cst(up
)) {
1443 struct isl_upoly_cst
*cst
;
1444 cst
= isl_upoly_as_cst(up
);
1447 isl_int_lcm(*d
, *d
, cst
->d
);
1451 rec
= isl_upoly_as_rec(up
);
1455 for (i
= 0; i
< rec
->n
; ++i
)
1456 upoly_update_den(rec
->p
[i
], d
);
1459 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial
*qp
, isl_int
*d
)
1461 isl_int_set_si(*d
, 1);
1464 upoly_update_den(qp
->upoly
, d
);
1467 __isl_give
struct isl_upoly
*isl_upoly_pow(isl_ctx
*ctx
, int pos
, int power
)
1470 struct isl_upoly
*up
;
1471 struct isl_upoly_rec
*rec
;
1472 struct isl_upoly_cst
*cst
;
1474 rec
= isl_upoly_alloc_rec(ctx
, pos
, 1 + power
);
1477 for (i
= 0; i
< 1 + power
; ++i
) {
1478 rec
->p
[i
] = isl_upoly_zero(ctx
);
1483 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1484 isl_int_set_si(cst
->n
, 1);
1488 isl_upoly_free(&rec
->up
);
1492 __isl_give isl_qpolynomial
*isl_qpolynomial_pow(__isl_take isl_dim
*dim
,
1495 struct isl_ctx
*ctx
;
1502 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_pow(ctx
, pos
, power
));
1505 __isl_give isl_qpolynomial
*isl_qpolynomial_var(__isl_take isl_dim
*dim
,
1506 enum isl_dim_type type
, unsigned pos
)
1511 isl_assert(dim
->ctx
, isl_dim_size(dim
, isl_dim_in
) == 0, goto error
);
1512 isl_assert(dim
->ctx
, pos
< isl_dim_size(dim
, type
), goto error
);
1514 if (type
== isl_dim_set
)
1515 pos
+= isl_dim_size(dim
, isl_dim_param
);
1517 return isl_qpolynomial_pow(dim
, pos
, 1);
1523 __isl_give isl_qpolynomial
*isl_qpolynomial_div_pow(__isl_take isl_div
*div
,
1526 struct isl_qpolynomial
*qp
= NULL
;
1527 struct isl_upoly_rec
*rec
;
1528 struct isl_upoly_cst
*cst
;
1534 isl_assert(div
->ctx
, div
->bmap
->n_div
== 1, goto error
);
1536 pos
= isl_dim_total(div
->bmap
->dim
);
1537 rec
= isl_upoly_alloc_rec(div
->ctx
, pos
, 1 + power
);
1538 qp
= isl_qpolynomial_alloc(isl_basic_map_get_dim(div
->bmap
), 1,
1543 isl_seq_cpy(qp
->div
->row
[0], div
->line
[0], qp
->div
->n_col
- 1);
1544 isl_int_set_si(qp
->div
->row
[0][qp
->div
->n_col
- 1], 0);
1546 for (i
= 0; i
< 1 + power
; ++i
) {
1547 rec
->p
[i
] = isl_upoly_zero(div
->ctx
);
1552 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1553 isl_int_set_si(cst
->n
, 1);
1559 isl_qpolynomial_free(qp
);
1564 __isl_give isl_qpolynomial
*isl_qpolynomial_div(__isl_take isl_div
*div
)
1566 return isl_qpolynomial_div_pow(div
, 1);
1569 __isl_give isl_qpolynomial
*isl_qpolynomial_rat_cst(__isl_take isl_dim
*dim
,
1570 const isl_int n
, const isl_int d
)
1572 struct isl_qpolynomial
*qp
;
1573 struct isl_upoly_cst
*cst
;
1575 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1579 cst
= isl_upoly_as_cst(qp
->upoly
);
1580 isl_int_set(cst
->n
, n
);
1581 isl_int_set(cst
->d
, d
);
1586 static int up_set_active(__isl_keep
struct isl_upoly
*up
, int *active
, int d
)
1588 struct isl_upoly_rec
*rec
;
1594 if (isl_upoly_is_cst(up
))
1598 active
[up
->var
] = 1;
1600 rec
= isl_upoly_as_rec(up
);
1601 for (i
= 0; i
< rec
->n
; ++i
)
1602 if (up_set_active(rec
->p
[i
], active
, d
) < 0)
1608 static int set_active(__isl_keep isl_qpolynomial
*qp
, int *active
)
1611 int d
= isl_dim_total(qp
->dim
);
1616 for (i
= 0; i
< d
; ++i
)
1617 for (j
= 0; j
< qp
->div
->n_row
; ++j
) {
1618 if (isl_int_is_zero(qp
->div
->row
[j
][2 + i
]))
1624 return up_set_active(qp
->upoly
, active
, d
);
1627 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial
*qp
,
1628 enum isl_dim_type type
, unsigned first
, unsigned n
)
1639 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
1641 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
1642 type
== isl_dim_set
, return -1);
1644 active
= isl_calloc_array(set
->ctx
, int, isl_dim_total(qp
->dim
));
1645 if (set_active(qp
, active
) < 0)
1648 if (type
== isl_dim_set
)
1649 first
+= isl_dim_size(qp
->dim
, isl_dim_param
);
1650 for (i
= 0; i
< n
; ++i
)
1651 if (active
[first
+ i
]) {
1664 __isl_give
struct isl_upoly
*isl_upoly_drop(__isl_take
struct isl_upoly
*up
,
1665 unsigned first
, unsigned n
)
1668 struct isl_upoly_rec
*rec
;
1672 if (n
== 0 || up
->var
< 0 || up
->var
< first
)
1674 if (up
->var
< first
+ n
) {
1675 up
= replace_by_constant_term(up
);
1676 return isl_upoly_drop(up
, first
, n
);
1678 up
= isl_upoly_cow(up
);
1682 rec
= isl_upoly_as_rec(up
);
1686 for (i
= 0; i
< rec
->n
; ++i
) {
1687 rec
->p
[i
] = isl_upoly_drop(rec
->p
[i
], first
, n
);
1698 __isl_give isl_qpolynomial
*isl_qpolynomial_drop_dims(
1699 __isl_take isl_qpolynomial
*qp
,
1700 enum isl_dim_type type
, unsigned first
, unsigned n
)
1704 if (n
== 0 && !isl_dim_get_tuple_name(qp
->dim
, type
))
1707 qp
= isl_qpolynomial_cow(qp
);
1711 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
1713 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
1714 type
== isl_dim_set
, goto error
);
1716 qp
->dim
= isl_dim_drop(qp
->dim
, type
, first
, n
);
1720 if (type
== isl_dim_set
)
1721 first
+= isl_dim_size(qp
->dim
, isl_dim_param
);
1723 qp
->div
= isl_mat_drop_cols(qp
->div
, 2 + first
, n
);
1727 qp
->upoly
= isl_upoly_drop(qp
->upoly
, first
, n
);
1733 isl_qpolynomial_free(qp
);
1738 #define PW isl_pw_qpolynomial
1740 #define EL isl_qpolynomial
1742 #define IS_ZERO is_zero
1746 #define ADD(d,e1,e2) isl_qpolynomial_add(e1,e2)
1748 #include <isl_pw_templ.c>
1751 #define UNION isl_union_pw_qpolynomial
1753 #define PART isl_pw_qpolynomial
1755 #define PARTS pw_qpolynomial
1757 #include <isl_union_templ.c>
1759 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial
*pwqp
)
1767 if (!isl_set_fast_is_universe(pwqp
->p
[0].set
))
1770 return isl_qpolynomial_is_one(pwqp
->p
[0].qp
);
1773 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_mul(
1774 __isl_take isl_pw_qpolynomial
*pwqp1
,
1775 __isl_take isl_pw_qpolynomial
*pwqp2
)
1778 struct isl_pw_qpolynomial
*res
;
1781 if (!pwqp1
|| !pwqp2
)
1784 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1787 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1788 isl_pw_qpolynomial_free(pwqp2
);
1792 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1793 isl_pw_qpolynomial_free(pwqp1
);
1797 if (isl_pw_qpolynomial_is_one(pwqp1
)) {
1798 isl_pw_qpolynomial_free(pwqp1
);
1802 if (isl_pw_qpolynomial_is_one(pwqp2
)) {
1803 isl_pw_qpolynomial_free(pwqp2
);
1807 n
= pwqp1
->n
* pwqp2
->n
;
1808 res
= isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1
->dim
), n
);
1810 for (i
= 0; i
< pwqp1
->n
; ++i
) {
1811 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1812 struct isl_set
*common
;
1813 struct isl_qpolynomial
*prod
;
1814 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
1815 isl_set_copy(pwqp2
->p
[j
].set
));
1816 if (isl_set_fast_is_empty(common
)) {
1817 isl_set_free(common
);
1821 prod
= isl_qpolynomial_mul(
1822 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
1823 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1825 res
= isl_pw_qpolynomial_add_piece(res
, common
, prod
);
1829 isl_pw_qpolynomial_free(pwqp1
);
1830 isl_pw_qpolynomial_free(pwqp2
);
1834 isl_pw_qpolynomial_free(pwqp1
);
1835 isl_pw_qpolynomial_free(pwqp2
);
1839 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_neg(
1840 __isl_take isl_pw_qpolynomial
*pwqp
)
1843 struct isl_pw_qpolynomial
*res
;
1849 if (isl_pw_qpolynomial_is_zero(pwqp
))
1852 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
1854 for (i
= 0; i
< pwqp
->n
; ++i
) {
1855 pwqp
->p
[i
].qp
= isl_qpolynomial_neg(pwqp
->p
[i
].qp
);
1862 isl_pw_qpolynomial_free(pwqp
);
1866 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_sub(
1867 __isl_take isl_pw_qpolynomial
*pwqp1
,
1868 __isl_take isl_pw_qpolynomial
*pwqp2
)
1870 return isl_pw_qpolynomial_add(pwqp1
, isl_pw_qpolynomial_neg(pwqp2
));
1873 __isl_give
struct isl_upoly
*isl_upoly_eval(
1874 __isl_take
struct isl_upoly
*up
, __isl_take isl_vec
*vec
)
1877 struct isl_upoly_rec
*rec
;
1878 struct isl_upoly
*res
;
1879 struct isl_upoly
*base
;
1881 if (isl_upoly_is_cst(up
)) {
1886 rec
= isl_upoly_as_rec(up
);
1890 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1892 base
= isl_upoly_rat_cst(up
->ctx
, vec
->el
[1 + up
->var
], vec
->el
[0]);
1894 res
= isl_upoly_eval(isl_upoly_copy(rec
->p
[rec
->n
- 1]),
1897 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1898 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1899 res
= isl_upoly_sum(res
,
1900 isl_upoly_eval(isl_upoly_copy(rec
->p
[i
]),
1901 isl_vec_copy(vec
)));
1904 isl_upoly_free(base
);
1914 __isl_give isl_qpolynomial
*isl_qpolynomial_eval(
1915 __isl_take isl_qpolynomial
*qp
, __isl_take isl_point
*pnt
)
1918 struct isl_upoly
*up
;
1923 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, qp
->dim
), goto error
);
1925 if (qp
->div
->n_row
== 0)
1926 ext
= isl_vec_copy(pnt
->vec
);
1929 unsigned dim
= isl_dim_total(qp
->dim
);
1930 ext
= isl_vec_alloc(qp
->dim
->ctx
, 1 + dim
+ qp
->div
->n_row
);
1934 isl_seq_cpy(ext
->el
, pnt
->vec
->el
, pnt
->vec
->size
);
1935 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1936 isl_seq_inner_product(qp
->div
->row
[i
] + 1, ext
->el
,
1937 1 + dim
+ i
, &ext
->el
[1+dim
+i
]);
1938 isl_int_fdiv_q(ext
->el
[1+dim
+i
], ext
->el
[1+dim
+i
],
1939 qp
->div
->row
[i
][0]);
1943 up
= isl_upoly_eval(isl_upoly_copy(qp
->upoly
), ext
);
1947 dim
= isl_dim_copy(qp
->dim
);
1948 isl_qpolynomial_free(qp
);
1949 isl_point_free(pnt
);
1951 return isl_qpolynomial_alloc(dim
, 0, up
);
1953 isl_qpolynomial_free(qp
);
1954 isl_point_free(pnt
);
1958 int isl_upoly_cmp(__isl_keep
struct isl_upoly_cst
*cst1
,
1959 __isl_keep
struct isl_upoly_cst
*cst2
)
1964 isl_int_mul(t
, cst1
->n
, cst2
->d
);
1965 isl_int_submul(t
, cst2
->n
, cst1
->d
);
1966 cmp
= isl_int_sgn(t
);
1971 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial
*qp1
,
1972 __isl_keep isl_qpolynomial
*qp2
)
1974 struct isl_upoly_cst
*cst1
, *cst2
;
1978 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), return -1);
1979 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), return -1);
1980 if (isl_qpolynomial_is_nan(qp1
))
1982 if (isl_qpolynomial_is_nan(qp2
))
1984 cst1
= isl_upoly_as_cst(qp1
->upoly
);
1985 cst2
= isl_upoly_as_cst(qp2
->upoly
);
1987 return isl_upoly_cmp(cst1
, cst2
) <= 0;
1990 __isl_give isl_qpolynomial
*isl_qpolynomial_min_cst(
1991 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1993 struct isl_upoly_cst
*cst1
, *cst2
;
1998 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
1999 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
2000 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2001 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2002 cmp
= isl_upoly_cmp(cst1
, cst2
);
2005 isl_qpolynomial_free(qp2
);
2007 isl_qpolynomial_free(qp1
);
2012 isl_qpolynomial_free(qp1
);
2013 isl_qpolynomial_free(qp2
);
2017 __isl_give isl_qpolynomial
*isl_qpolynomial_max_cst(
2018 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
2020 struct isl_upoly_cst
*cst1
, *cst2
;
2025 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
2026 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
2027 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2028 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2029 cmp
= isl_upoly_cmp(cst1
, cst2
);
2032 isl_qpolynomial_free(qp2
);
2034 isl_qpolynomial_free(qp1
);
2039 isl_qpolynomial_free(qp1
);
2040 isl_qpolynomial_free(qp2
);
2044 __isl_give isl_qpolynomial
*isl_qpolynomial_insert_dims(
2045 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
,
2046 unsigned first
, unsigned n
)
2055 qp
= isl_qpolynomial_cow(qp
);
2059 isl_assert(qp
->div
->ctx
, first
<= isl_dim_size(qp
->dim
, type
),
2062 g_pos
= pos(qp
->dim
, type
) + first
;
2064 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + g_pos
, n
);
2068 total
= qp
->div
->n_col
- 2;
2069 if (total
> g_pos
) {
2071 exp
= isl_alloc_array(qp
->div
->ctx
, int, total
- g_pos
);
2074 for (i
= 0; i
< total
- g_pos
; ++i
)
2076 qp
->upoly
= expand(qp
->upoly
, exp
, g_pos
);
2082 qp
->dim
= isl_dim_insert(qp
->dim
, type
, first
, n
);
2088 isl_qpolynomial_free(qp
);
2092 __isl_give isl_qpolynomial
*isl_qpolynomial_add_dims(
2093 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
, unsigned n
)
2097 pos
= isl_qpolynomial_dim(qp
, type
);
2099 return isl_qpolynomial_insert_dims(qp
, type
, pos
, n
);
2102 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add_dims(
2103 __isl_take isl_pw_qpolynomial
*pwqp
,
2104 enum isl_dim_type type
, unsigned n
)
2111 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
2115 pwqp
->dim
= isl_dim_add(pwqp
->dim
, type
, n
);
2119 for (i
= 0; i
< pwqp
->n
; ++i
) {
2120 pwqp
->p
[i
].set
= isl_set_add(pwqp
->p
[i
].set
, type
, n
);
2121 if (!pwqp
->p
[i
].set
)
2123 pwqp
->p
[i
].qp
= isl_qpolynomial_add_dims(pwqp
->p
[i
].qp
, type
, n
);
2130 isl_pw_qpolynomial_free(pwqp
);
2134 static int *reordering_move(isl_ctx
*ctx
,
2135 unsigned len
, unsigned dst
, unsigned src
, unsigned n
)
2140 reordering
= isl_alloc_array(ctx
, int, len
);
2145 for (i
= 0; i
< dst
; ++i
)
2147 for (i
= 0; i
< n
; ++i
)
2148 reordering
[src
+ i
] = dst
+ i
;
2149 for (i
= 0; i
< src
- dst
; ++i
)
2150 reordering
[dst
+ i
] = dst
+ n
+ i
;
2151 for (i
= 0; i
< len
- src
- n
; ++i
)
2152 reordering
[src
+ n
+ i
] = src
+ n
+ i
;
2154 for (i
= 0; i
< src
; ++i
)
2156 for (i
= 0; i
< n
; ++i
)
2157 reordering
[src
+ i
] = dst
+ i
;
2158 for (i
= 0; i
< dst
- src
; ++i
)
2159 reordering
[src
+ n
+ i
] = src
+ i
;
2160 for (i
= 0; i
< len
- dst
- n
; ++i
)
2161 reordering
[dst
+ n
+ i
] = dst
+ n
+ i
;
2167 static __isl_give
struct isl_upoly
*reorder(__isl_take
struct isl_upoly
*up
,
2171 struct isl_upoly_rec
*rec
;
2172 struct isl_upoly
*base
;
2173 struct isl_upoly
*res
;
2175 if (isl_upoly_is_cst(up
))
2178 rec
= isl_upoly_as_rec(up
);
2182 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
2184 base
= isl_upoly_pow(up
->ctx
, r
[up
->var
], 1);
2185 res
= reorder(isl_upoly_copy(rec
->p
[rec
->n
- 1]), r
);
2187 for (i
= rec
->n
- 2; i
>= 0; --i
) {
2188 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
2189 res
= isl_upoly_sum(res
, reorder(isl_upoly_copy(rec
->p
[i
]), r
));
2192 isl_upoly_free(base
);
2201 __isl_give isl_qpolynomial
*isl_qpolynomial_move_dims(
2202 __isl_take isl_qpolynomial
*qp
,
2203 enum isl_dim_type dst_type
, unsigned dst_pos
,
2204 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
2210 qp
= isl_qpolynomial_cow(qp
);
2214 isl_assert(qp
->dim
->ctx
, src_pos
+ n
<= isl_dim_size(qp
->dim
, src_type
),
2217 g_dst_pos
= pos(qp
->dim
, dst_type
) + dst_pos
;
2218 g_src_pos
= pos(qp
->dim
, src_type
) + src_pos
;
2219 if (dst_type
> src_type
)
2222 qp
->div
= isl_mat_move_cols(qp
->div
, 2 + g_dst_pos
, 2 + g_src_pos
, n
);
2223 qp
->div
= sort_divs(qp
->div
);
2227 reordering
= reordering_move(qp
->dim
->ctx
,
2228 qp
->div
->n_col
- 2, g_dst_pos
, g_src_pos
, n
);
2232 qp
->upoly
= reorder(qp
->upoly
, reordering
);
2237 qp
->dim
= isl_dim_move(qp
->dim
, dst_type
, dst_pos
, src_type
, src_pos
, n
);
2243 isl_qpolynomial_free(qp
);
2247 __isl_give
struct isl_upoly
*isl_upoly_from_affine(isl_ctx
*ctx
, isl_int
*f
,
2248 isl_int denom
, unsigned len
)
2251 struct isl_upoly
*up
;
2253 isl_assert(ctx
, len
>= 1, return NULL
);
2255 up
= isl_upoly_rat_cst(ctx
, f
[0], denom
);
2256 for (i
= 0; i
< len
- 1; ++i
) {
2257 struct isl_upoly
*t
;
2258 struct isl_upoly
*c
;
2260 if (isl_int_is_zero(f
[1 + i
]))
2263 c
= isl_upoly_rat_cst(ctx
, f
[1 + i
], denom
);
2264 t
= isl_upoly_pow(ctx
, i
, 1);
2265 t
= isl_upoly_mul(c
, t
);
2266 up
= isl_upoly_sum(up
, t
);
2272 __isl_give isl_qpolynomial
*isl_qpolynomial_from_affine(__isl_take isl_dim
*dim
,
2273 isl_int
*f
, isl_int denom
)
2275 struct isl_upoly
*up
;
2280 up
= isl_upoly_from_affine(dim
->ctx
, f
, denom
, 1 + isl_dim_total(dim
));
2282 return isl_qpolynomial_alloc(dim
, 0, up
);
2285 __isl_give isl_qpolynomial
*isl_qpolynomial_from_constraint(
2286 __isl_take isl_constraint
*c
, enum isl_dim_type type
, unsigned pos
)
2290 struct isl_upoly
*up
;
2291 isl_qpolynomial
*qp
;
2297 isl_int_init(denom
);
2299 isl_constraint_get_coefficient(c
, type
, pos
, &denom
);
2300 isl_constraint_set_coefficient(c
, type
, pos
, c
->ctx
->zero
);
2301 sgn
= isl_int_sgn(denom
);
2302 isl_int_abs(denom
, denom
);
2303 up
= isl_upoly_from_affine(c
->ctx
, c
->line
[0], denom
,
2304 1 + isl_constraint_dim(c
, isl_dim_all
));
2306 isl_int_neg(denom
, denom
);
2307 isl_constraint_set_coefficient(c
, type
, pos
, denom
);
2309 dim
= isl_dim_copy(c
->bmap
->dim
);
2311 isl_int_clear(denom
);
2312 isl_constraint_free(c
);
2314 qp
= isl_qpolynomial_alloc(dim
, 0, up
);
2316 qp
= isl_qpolynomial_neg(qp
);
2320 __isl_give
struct isl_upoly
*isl_upoly_subs(__isl_take
struct isl_upoly
*up
,
2321 unsigned first
, unsigned n
, __isl_keep
struct isl_upoly
**subs
)
2324 struct isl_upoly_rec
*rec
;
2325 struct isl_upoly
*base
, *res
;
2330 if (isl_upoly_is_cst(up
))
2333 if (up
->var
< first
)
2336 rec
= isl_upoly_as_rec(up
);
2340 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
2342 if (up
->var
>= first
+ n
)
2343 base
= isl_upoly_pow(up
->ctx
, up
->var
, 1);
2345 base
= isl_upoly_copy(subs
[up
->var
- first
]);
2347 res
= isl_upoly_subs(isl_upoly_copy(rec
->p
[rec
->n
- 1]), first
, n
, subs
);
2348 for (i
= rec
->n
- 2; i
>= 0; --i
) {
2349 struct isl_upoly
*t
;
2350 t
= isl_upoly_subs(isl_upoly_copy(rec
->p
[i
]), first
, n
, subs
);
2351 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
2352 res
= isl_upoly_sum(res
, t
);
2355 isl_upoly_free(base
);
2364 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
2365 * in "qp" by subs[i].
2367 __isl_give isl_qpolynomial
*isl_qpolynomial_substitute(
2368 __isl_take isl_qpolynomial
*qp
,
2369 enum isl_dim_type type
, unsigned first
, unsigned n
,
2370 __isl_keep isl_qpolynomial
**subs
)
2373 struct isl_upoly
**ups
;
2378 qp
= isl_qpolynomial_cow(qp
);
2381 for (i
= 0; i
< n
; ++i
)
2385 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
2388 for (i
= 0; i
< n
; ++i
)
2389 isl_assert(qp
->dim
->ctx
, isl_dim_equal(qp
->dim
, subs
[i
]->dim
),
2392 isl_assert(qp
->dim
->ctx
, qp
->div
->n_row
== 0, goto error
);
2393 for (i
= 0; i
< n
; ++i
)
2394 isl_assert(qp
->dim
->ctx
, subs
[i
]->div
->n_row
== 0, goto error
);
2396 first
+= pos(qp
->dim
, type
);
2398 ups
= isl_alloc_array(qp
->dim
->ctx
, struct isl_upoly
*, n
);
2401 for (i
= 0; i
< n
; ++i
)
2402 ups
[i
] = subs
[i
]->upoly
;
2404 qp
->upoly
= isl_upoly_subs(qp
->upoly
, first
, n
, ups
);
2413 isl_qpolynomial_free(qp
);
2417 __isl_give isl_basic_set
*add_div_constraints(__isl_take isl_basic_set
*bset
,
2418 __isl_take isl_mat
*div
)
2426 bset
= isl_basic_set_extend_constraints(bset
, 0, 2 * div
->n_row
);
2429 total
= isl_basic_set_total_dim(bset
);
2430 for (i
= 0; i
< div
->n_row
; ++i
)
2431 if (isl_basic_set_add_div_constraints_var(bset
,
2432 total
- div
->n_row
+ i
, div
->row
[i
]) < 0)
2439 isl_basic_set_free(bset
);
2443 /* Extend "bset" with extra set dimensions for each integer division
2444 * in "qp" and then call "fn" with the extended bset and the polynomial
2445 * that results from replacing each of the integer divisions by the
2446 * corresponding extra set dimension.
2448 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial
*qp
,
2449 __isl_keep isl_basic_set
*bset
,
2450 int (*fn
)(__isl_take isl_basic_set
*bset
,
2451 __isl_take isl_qpolynomial
*poly
, void *user
), void *user
)
2455 isl_qpolynomial
*poly
;
2459 if (qp
->div
->n_row
== 0)
2460 return fn(isl_basic_set_copy(bset
), isl_qpolynomial_copy(qp
),
2463 div
= isl_mat_copy(qp
->div
);
2464 dim
= isl_dim_copy(qp
->dim
);
2465 dim
= isl_dim_add(dim
, isl_dim_set
, qp
->div
->n_row
);
2466 poly
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_copy(qp
->upoly
));
2467 bset
= isl_basic_set_copy(bset
);
2468 bset
= isl_basic_set_add(bset
, isl_dim_set
, qp
->div
->n_row
);
2469 bset
= add_div_constraints(bset
, div
);
2471 return fn(bset
, poly
, user
);
2476 /* Return total degree in variables first (inclusive) up to last (exclusive).
2478 int isl_upoly_degree(__isl_keep
struct isl_upoly
*up
, int first
, int last
)
2482 struct isl_upoly_rec
*rec
;
2486 if (isl_upoly_is_zero(up
))
2488 if (isl_upoly_is_cst(up
) || up
->var
< first
)
2491 rec
= isl_upoly_as_rec(up
);
2495 for (i
= 0; i
< rec
->n
; ++i
) {
2498 if (isl_upoly_is_zero(rec
->p
[i
]))
2500 d
= isl_upoly_degree(rec
->p
[i
], first
, last
);
2510 /* Return total degree in set variables.
2512 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial
*poly
)
2520 ovar
= isl_dim_offset(poly
->dim
, isl_dim_set
);
2521 nvar
= isl_dim_size(poly
->dim
, isl_dim_set
);
2522 return isl_upoly_degree(poly
->upoly
, ovar
, ovar
+ nvar
);
2525 __isl_give
struct isl_upoly
*isl_upoly_coeff(__isl_keep
struct isl_upoly
*up
,
2526 unsigned pos
, int deg
)
2529 struct isl_upoly_rec
*rec
;
2534 if (isl_upoly_is_cst(up
) || up
->var
< pos
) {
2536 return isl_upoly_copy(up
);
2538 return isl_upoly_zero(up
->ctx
);
2541 rec
= isl_upoly_as_rec(up
);
2545 if (up
->var
== pos
) {
2547 return isl_upoly_copy(rec
->p
[deg
]);
2549 return isl_upoly_zero(up
->ctx
);
2552 up
= isl_upoly_copy(up
);
2553 up
= isl_upoly_cow(up
);
2554 rec
= isl_upoly_as_rec(up
);
2558 for (i
= 0; i
< rec
->n
; ++i
) {
2559 struct isl_upoly
*t
;
2560 t
= isl_upoly_coeff(rec
->p
[i
], pos
, deg
);
2563 isl_upoly_free(rec
->p
[i
]);
2573 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
2575 __isl_give isl_qpolynomial
*isl_qpolynomial_coeff(
2576 __isl_keep isl_qpolynomial
*qp
,
2577 enum isl_dim_type type
, unsigned t_pos
, int deg
)
2580 struct isl_upoly
*up
;
2586 isl_assert(qp
->div
->ctx
, t_pos
< isl_dim_size(qp
->dim
, type
),
2589 g_pos
= pos(qp
->dim
, type
) + t_pos
;
2590 up
= isl_upoly_coeff(qp
->upoly
, g_pos
, deg
);
2592 c
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
, up
);
2595 isl_mat_free(c
->div
);
2596 c
->div
= isl_mat_copy(qp
->div
);
2601 isl_qpolynomial_free(c
);
2605 /* Homogenize the polynomial in the variables first (inclusive) up to
2606 * last (exclusive) by inserting powers of variable first.
2607 * Variable first is assumed not to appear in the input.
2609 __isl_give
struct isl_upoly
*isl_upoly_homogenize(
2610 __isl_take
struct isl_upoly
*up
, int deg
, int target
,
2611 int first
, int last
)
2614 struct isl_upoly_rec
*rec
;
2618 if (isl_upoly_is_zero(up
))
2622 if (isl_upoly_is_cst(up
) || up
->var
< first
) {
2623 struct isl_upoly
*hom
;
2625 hom
= isl_upoly_pow(up
->ctx
, first
, target
- deg
);
2628 rec
= isl_upoly_as_rec(hom
);
2629 rec
->p
[target
- deg
] = isl_upoly_mul(rec
->p
[target
- deg
], up
);
2634 up
= isl_upoly_cow(up
);
2635 rec
= isl_upoly_as_rec(up
);
2639 for (i
= 0; i
< rec
->n
; ++i
) {
2640 if (isl_upoly_is_zero(rec
->p
[i
]))
2642 rec
->p
[i
] = isl_upoly_homogenize(rec
->p
[i
],
2643 up
->var
< last
? deg
+ i
: i
, target
,
2655 /* Homogenize the polynomial in the set variables by introducing
2656 * powers of an extra set variable at position 0.
2658 __isl_give isl_qpolynomial
*isl_qpolynomial_homogenize(
2659 __isl_take isl_qpolynomial
*poly
)
2663 int deg
= isl_qpolynomial_degree(poly
);
2668 poly
= isl_qpolynomial_insert_dims(poly
, isl_dim_set
, 0, 1);
2669 poly
= isl_qpolynomial_cow(poly
);
2673 ovar
= isl_dim_offset(poly
->dim
, isl_dim_set
);
2674 nvar
= isl_dim_size(poly
->dim
, isl_dim_set
);
2675 poly
->upoly
= isl_upoly_homogenize(poly
->upoly
, 0, deg
,
2682 isl_qpolynomial_free(poly
);
2686 __isl_give isl_term
*isl_term_alloc(__isl_take isl_dim
*dim
,
2687 __isl_take isl_mat
*div
)
2695 n
= isl_dim_total(dim
) + div
->n_row
;
2697 term
= isl_calloc(dim
->ctx
, struct isl_term
,
2698 sizeof(struct isl_term
) + (n
- 1) * sizeof(int));
2705 isl_int_init(term
->n
);
2706 isl_int_init(term
->d
);
2715 __isl_give isl_term
*isl_term_copy(__isl_keep isl_term
*term
)
2724 __isl_give isl_term
*isl_term_dup(__isl_keep isl_term
*term
)
2733 total
= isl_dim_total(term
->dim
) + term
->div
->n_row
;
2735 dup
= isl_term_alloc(isl_dim_copy(term
->dim
), isl_mat_copy(term
->div
));
2739 isl_int_set(dup
->n
, term
->n
);
2740 isl_int_set(dup
->d
, term
->d
);
2742 for (i
= 0; i
< total
; ++i
)
2743 dup
->pow
[i
] = term
->pow
[i
];
2748 __isl_give isl_term
*isl_term_cow(__isl_take isl_term
*term
)
2756 return isl_term_dup(term
);
2759 void isl_term_free(__isl_take isl_term
*term
)
2764 if (--term
->ref
> 0)
2767 isl_dim_free(term
->dim
);
2768 isl_mat_free(term
->div
);
2769 isl_int_clear(term
->n
);
2770 isl_int_clear(term
->d
);
2774 unsigned isl_term_dim(__isl_keep isl_term
*term
, enum isl_dim_type type
)
2782 case isl_dim_out
: return isl_dim_size(term
->dim
, type
);
2783 case isl_dim_div
: return term
->div
->n_row
;
2784 case isl_dim_all
: return isl_dim_total(term
->dim
) + term
->div
->n_row
;
2789 isl_ctx
*isl_term_get_ctx(__isl_keep isl_term
*term
)
2791 return term
? term
->dim
->ctx
: NULL
;
2794 void isl_term_get_num(__isl_keep isl_term
*term
, isl_int
*n
)
2798 isl_int_set(*n
, term
->n
);
2801 void isl_term_get_den(__isl_keep isl_term
*term
, isl_int
*d
)
2805 isl_int_set(*d
, term
->d
);
2808 int isl_term_get_exp(__isl_keep isl_term
*term
,
2809 enum isl_dim_type type
, unsigned pos
)
2814 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, type
), return -1);
2816 if (type
>= isl_dim_set
)
2817 pos
+= isl_dim_size(term
->dim
, isl_dim_param
);
2818 if (type
>= isl_dim_div
)
2819 pos
+= isl_dim_size(term
->dim
, isl_dim_set
);
2821 return term
->pow
[pos
];
2824 __isl_give isl_div
*isl_term_get_div(__isl_keep isl_term
*term
, unsigned pos
)
2826 isl_basic_map
*bmap
;
2833 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, isl_dim_div
),
2836 total
= term
->div
->n_col
- term
->div
->n_row
- 2;
2837 /* No nested divs for now */
2838 isl_assert(term
->dim
->ctx
,
2839 isl_seq_first_non_zero(term
->div
->row
[pos
] + 2 + total
,
2840 term
->div
->n_row
) == -1,
2843 bmap
= isl_basic_map_alloc_dim(isl_dim_copy(term
->dim
), 1, 0, 0);
2844 if ((k
= isl_basic_map_alloc_div(bmap
)) < 0)
2847 isl_seq_cpy(bmap
->div
[k
], term
->div
->row
[pos
], 2 + total
);
2849 return isl_basic_map_div(bmap
, k
);
2851 isl_basic_map_free(bmap
);
2855 __isl_give isl_term
*isl_upoly_foreach_term(__isl_keep
struct isl_upoly
*up
,
2856 int (*fn
)(__isl_take isl_term
*term
, void *user
),
2857 __isl_take isl_term
*term
, void *user
)
2860 struct isl_upoly_rec
*rec
;
2865 if (isl_upoly_is_zero(up
))
2868 isl_assert(up
->ctx
, !isl_upoly_is_nan(up
), goto error
);
2869 isl_assert(up
->ctx
, !isl_upoly_is_infty(up
), goto error
);
2870 isl_assert(up
->ctx
, !isl_upoly_is_neginfty(up
), goto error
);
2872 if (isl_upoly_is_cst(up
)) {
2873 struct isl_upoly_cst
*cst
;
2874 cst
= isl_upoly_as_cst(up
);
2877 term
= isl_term_cow(term
);
2880 isl_int_set(term
->n
, cst
->n
);
2881 isl_int_set(term
->d
, cst
->d
);
2882 if (fn(isl_term_copy(term
), user
) < 0)
2887 rec
= isl_upoly_as_rec(up
);
2891 for (i
= 0; i
< rec
->n
; ++i
) {
2892 term
= isl_term_cow(term
);
2895 term
->pow
[up
->var
] = i
;
2896 term
= isl_upoly_foreach_term(rec
->p
[i
], fn
, term
, user
);
2900 term
->pow
[up
->var
] = 0;
2904 isl_term_free(term
);
2908 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial
*qp
,
2909 int (*fn
)(__isl_take isl_term
*term
, void *user
), void *user
)
2916 term
= isl_term_alloc(isl_dim_copy(qp
->dim
), isl_mat_copy(qp
->div
));
2920 term
= isl_upoly_foreach_term(qp
->upoly
, fn
, term
, user
);
2922 isl_term_free(term
);
2924 return term
? 0 : -1;
2927 __isl_give isl_qpolynomial
*isl_qpolynomial_from_term(__isl_take isl_term
*term
)
2929 struct isl_upoly
*up
;
2930 isl_qpolynomial
*qp
;
2936 n
= isl_dim_total(term
->dim
) + term
->div
->n_row
;
2938 up
= isl_upoly_rat_cst(term
->dim
->ctx
, term
->n
, term
->d
);
2939 for (i
= 0; i
< n
; ++i
) {
2942 up
= isl_upoly_mul(up
,
2943 isl_upoly_pow(term
->dim
->ctx
, i
, term
->pow
[i
]));
2946 qp
= isl_qpolynomial_alloc(isl_dim_copy(term
->dim
), term
->div
->n_row
, up
);
2949 isl_mat_free(qp
->div
);
2950 qp
->div
= isl_mat_copy(term
->div
);
2954 isl_term_free(term
);
2957 isl_qpolynomial_free(qp
);
2958 isl_term_free(term
);
2962 __isl_give isl_qpolynomial
*isl_qpolynomial_lift(__isl_take isl_qpolynomial
*qp
,
2963 __isl_take isl_dim
*dim
)
2972 if (isl_dim_equal(qp
->dim
, dim
)) {
2977 qp
= isl_qpolynomial_cow(qp
);
2981 extra
= isl_dim_size(dim
, isl_dim_set
) -
2982 isl_dim_size(qp
->dim
, isl_dim_set
);
2983 total
= isl_dim_total(qp
->dim
);
2984 if (qp
->div
->n_row
) {
2987 exp
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
2990 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
2992 qp
->upoly
= expand(qp
->upoly
, exp
, total
);
2997 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + total
, extra
);
3000 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
3001 isl_seq_clr(qp
->div
->row
[i
] + 2 + total
, extra
);
3003 isl_dim_free(qp
->dim
);
3009 isl_qpolynomial_free(qp
);
3013 /* For each parameter or variable that does not appear in qp,
3014 * first eliminate the variable from all constraints and then set it to zero.
3016 static __isl_give isl_set
*fix_inactive(__isl_take isl_set
*set
,
3017 __isl_keep isl_qpolynomial
*qp
)
3028 d
= isl_dim_total(set
->dim
);
3029 active
= isl_calloc_array(set
->ctx
, int, d
);
3030 if (set_active(qp
, active
) < 0)
3033 for (i
= 0; i
< d
; ++i
)
3042 nparam
= isl_dim_size(set
->dim
, isl_dim_param
);
3043 nvar
= isl_dim_size(set
->dim
, isl_dim_set
);
3044 for (i
= 0; i
< nparam
; ++i
) {
3047 set
= isl_set_eliminate(set
, isl_dim_param
, i
, 1);
3048 set
= isl_set_fix_si(set
, isl_dim_param
, i
, 0);
3050 for (i
= 0; i
< nvar
; ++i
) {
3051 if (active
[nparam
+ i
])
3053 set
= isl_set_eliminate(set
, isl_dim_set
, i
, 1);
3054 set
= isl_set_fix_si(set
, isl_dim_set
, i
, 0);
3066 struct isl_opt_data
{
3067 isl_qpolynomial
*qp
;
3069 isl_qpolynomial
*opt
;
3073 static int opt_fn(__isl_take isl_point
*pnt
, void *user
)
3075 struct isl_opt_data
*data
= (struct isl_opt_data
*)user
;
3076 isl_qpolynomial
*val
;
3078 val
= isl_qpolynomial_eval(isl_qpolynomial_copy(data
->qp
), pnt
);
3082 } else if (data
->max
) {
3083 data
->opt
= isl_qpolynomial_max_cst(data
->opt
, val
);
3085 data
->opt
= isl_qpolynomial_min_cst(data
->opt
, val
);
3091 __isl_give isl_qpolynomial
*isl_qpolynomial_opt_on_domain(
3092 __isl_take isl_qpolynomial
*qp
, __isl_take isl_set
*set
, int max
)
3094 struct isl_opt_data data
= { NULL
, 1, NULL
, max
};
3099 if (isl_upoly_is_cst(qp
->upoly
)) {
3104 set
= fix_inactive(set
, qp
);
3107 if (isl_set_foreach_point(set
, opt_fn
, &data
) < 0)
3111 data
.opt
= isl_qpolynomial_zero(isl_qpolynomial_get_dim(qp
));
3114 isl_qpolynomial_free(qp
);
3118 isl_qpolynomial_free(qp
);
3119 isl_qpolynomial_free(data
.opt
);
3123 __isl_give isl_qpolynomial
*isl_qpolynomial_morph(__isl_take isl_qpolynomial
*qp
,
3124 __isl_take isl_morph
*morph
)
3128 struct isl_upoly
*up
;
3130 struct isl_upoly
**subs
;
3133 qp
= isl_qpolynomial_cow(qp
);
3138 isl_assert(ctx
, isl_dim_equal(qp
->dim
, morph
->dom
->dim
), goto error
);
3140 subs
= isl_calloc_array(ctx
, struct isl_upoly
*, morph
->inv
->n_row
- 1);
3144 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
)
3145 subs
[i
] = isl_upoly_from_affine(ctx
, morph
->inv
->row
[1 + i
],
3146 morph
->inv
->row
[0][0], morph
->inv
->n_col
);
3148 qp
->upoly
= isl_upoly_subs(qp
->upoly
, 0, morph
->inv
->n_row
- 1, subs
);
3150 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
)
3151 isl_upoly_free(subs
[i
]);
3154 mat
= isl_mat_diagonal(isl_mat_identity(ctx
, 1), isl_mat_copy(morph
->inv
));
3155 mat
= isl_mat_diagonal(mat
, isl_mat_identity(ctx
, qp
->div
->n_row
));
3156 qp
->div
= isl_mat_product(qp
->div
, mat
);
3157 isl_dim_free(qp
->dim
);
3158 qp
->dim
= isl_dim_copy(morph
->ran
->dim
);
3160 if (!qp
->upoly
|| !qp
->div
|| !qp
->dim
)
3163 isl_morph_free(morph
);
3167 isl_qpolynomial_free(qp
);
3168 isl_morph_free(morph
);
3172 static int neg_entry(void **entry
, void *user
)
3174 isl_pw_qpolynomial
**pwqp
= (isl_pw_qpolynomial
**)entry
;
3176 *pwqp
= isl_pw_qpolynomial_neg(*pwqp
);
3178 return *pwqp
? 0 : -1;
3181 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_neg(
3182 __isl_take isl_union_pw_qpolynomial
*upwqp
)
3184 upwqp
= isl_union_pw_qpolynomial_cow(upwqp
);
3188 if (isl_hash_table_foreach(upwqp
->dim
->ctx
, &upwqp
->table
,
3189 &neg_entry
, NULL
) < 0)
3194 isl_union_pw_qpolynomial_free(upwqp
);
3198 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_sub(
3199 __isl_take isl_union_pw_qpolynomial
*upwqp1
,
3200 __isl_take isl_union_pw_qpolynomial
*upwqp2
)
3202 return isl_union_pw_qpolynomial_add(upwqp1
,
3203 isl_union_pw_qpolynomial_neg(upwqp2
));
3206 static int mul_entry(void **entry
, void *user
)
3208 struct isl_union_pw_qpolynomial_match_bin_data
*data
= user
;
3210 struct isl_hash_table_entry
*entry2
;
3211 isl_pw_qpolynomial
*pwpq
= *entry
;
3214 hash
= isl_dim_get_hash(pwpq
->dim
);
3215 entry2
= isl_hash_table_find(data
->upwqp2
->dim
->ctx
,
3216 &data
->upwqp2
->table
,
3217 hash
, &has_dim
, pwpq
->dim
, 0);
3221 pwpq
= isl_pw_qpolynomial_copy(pwpq
);
3222 pwpq
= isl_pw_qpolynomial_mul(pwpq
,
3223 isl_pw_qpolynomial_copy(entry2
->data
));
3225 empty
= isl_pw_qpolynomial_is_zero(pwpq
);
3227 isl_pw_qpolynomial_free(pwpq
);
3231 isl_pw_qpolynomial_free(pwpq
);
3235 data
->res
= isl_union_pw_qpolynomial_add_pw_qpolynomial(data
->res
, pwpq
);
3240 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_mul(
3241 __isl_take isl_union_pw_qpolynomial
*upwqp1
,
3242 __isl_take isl_union_pw_qpolynomial
*upwqp2
)
3244 return match_bin_op(upwqp1
, upwqp2
, &mul_entry
);