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_zero(__isl_keep
struct isl_upoly
*up
)
37 struct isl_upoly_cst
*cst
;
41 if (!isl_upoly_is_cst(up
))
44 cst
= isl_upoly_as_cst(up
);
48 return isl_int_is_zero(cst
->n
) && isl_int_is_pos(cst
->d
);
51 int isl_upoly_is_nan(__isl_keep
struct isl_upoly
*up
)
53 struct isl_upoly_cst
*cst
;
57 if (!isl_upoly_is_cst(up
))
60 cst
= isl_upoly_as_cst(up
);
64 return isl_int_is_zero(cst
->n
) && isl_int_is_zero(cst
->d
);
67 int isl_upoly_is_infty(__isl_keep
struct isl_upoly
*up
)
69 struct isl_upoly_cst
*cst
;
73 if (!isl_upoly_is_cst(up
))
76 cst
= isl_upoly_as_cst(up
);
80 return isl_int_is_pos(cst
->n
) && isl_int_is_zero(cst
->d
);
83 int isl_upoly_is_neginfty(__isl_keep
struct isl_upoly
*up
)
85 struct isl_upoly_cst
*cst
;
89 if (!isl_upoly_is_cst(up
))
92 cst
= isl_upoly_as_cst(up
);
96 return isl_int_is_neg(cst
->n
) && isl_int_is_zero(cst
->d
);
99 int isl_upoly_is_one(__isl_keep
struct isl_upoly
*up
)
101 struct isl_upoly_cst
*cst
;
105 if (!isl_upoly_is_cst(up
))
108 cst
= isl_upoly_as_cst(up
);
112 return isl_int_eq(cst
->n
, cst
->d
) && isl_int_is_pos(cst
->d
);
115 int isl_upoly_is_negone(__isl_keep
struct isl_upoly
*up
)
117 struct isl_upoly_cst
*cst
;
121 if (!isl_upoly_is_cst(up
))
124 cst
= isl_upoly_as_cst(up
);
128 return isl_int_is_negone(cst
->n
) && isl_int_is_one(cst
->d
);
131 __isl_give
struct isl_upoly_cst
*isl_upoly_cst_alloc(struct isl_ctx
*ctx
)
133 struct isl_upoly_cst
*cst
;
135 cst
= isl_alloc_type(ctx
, struct isl_upoly_cst
);
144 isl_int_init(cst
->n
);
145 isl_int_init(cst
->d
);
150 __isl_give
struct isl_upoly
*isl_upoly_zero(struct isl_ctx
*ctx
)
152 struct isl_upoly_cst
*cst
;
154 cst
= isl_upoly_cst_alloc(ctx
);
158 isl_int_set_si(cst
->n
, 0);
159 isl_int_set_si(cst
->d
, 1);
164 __isl_give
struct isl_upoly
*isl_upoly_infty(struct isl_ctx
*ctx
)
166 struct isl_upoly_cst
*cst
;
168 cst
= isl_upoly_cst_alloc(ctx
);
172 isl_int_set_si(cst
->n
, 1);
173 isl_int_set_si(cst
->d
, 0);
178 __isl_give
struct isl_upoly
*isl_upoly_nan(struct isl_ctx
*ctx
)
180 struct isl_upoly_cst
*cst
;
182 cst
= isl_upoly_cst_alloc(ctx
);
186 isl_int_set_si(cst
->n
, 0);
187 isl_int_set_si(cst
->d
, 0);
192 __isl_give
struct isl_upoly
*isl_upoly_rat_cst(struct isl_ctx
*ctx
,
193 isl_int n
, isl_int d
)
195 struct isl_upoly_cst
*cst
;
197 cst
= isl_upoly_cst_alloc(ctx
);
201 isl_int_set(cst
->n
, n
);
202 isl_int_set(cst
->d
, d
);
207 __isl_give
struct isl_upoly_rec
*isl_upoly_alloc_rec(struct isl_ctx
*ctx
,
210 struct isl_upoly_rec
*rec
;
212 isl_assert(ctx
, var
>= 0, return NULL
);
213 isl_assert(ctx
, size
>= 0, return NULL
);
214 rec
= isl_calloc(dim
->ctx
, struct isl_upoly_rec
,
215 sizeof(struct isl_upoly_rec
) +
216 (size
- 1) * sizeof(struct isl_upoly
*));
230 isl_upoly_free(&rec
->up
);
234 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial
*qp
)
236 struct isl_upoly_cst
*cst
;
241 return isl_upoly_is_zero(qp
->upoly
);
244 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial
*qp
)
246 struct isl_upoly_cst
*cst
;
251 return isl_upoly_is_one(qp
->upoly
);
254 static void upoly_free_cst(__isl_take
struct isl_upoly_cst
*cst
)
256 isl_int_clear(cst
->n
);
257 isl_int_clear(cst
->d
);
260 static void upoly_free_rec(__isl_take
struct isl_upoly_rec
*rec
)
264 for (i
= 0; i
< rec
->n
; ++i
)
265 isl_upoly_free(rec
->p
[i
]);
268 __isl_give
struct isl_upoly
*isl_upoly_copy(__isl_keep
struct isl_upoly
*up
)
277 __isl_give
struct isl_upoly
*isl_upoly_dup_cst(__isl_keep
struct isl_upoly
*up
)
279 struct isl_upoly_cst
*cst
;
280 struct isl_upoly_cst
*dup
;
282 cst
= isl_upoly_as_cst(up
);
286 dup
= isl_upoly_as_cst(isl_upoly_zero(up
->ctx
));
289 isl_int_set(dup
->n
, cst
->n
);
290 isl_int_set(dup
->d
, cst
->d
);
295 __isl_give
struct isl_upoly
*isl_upoly_dup_rec(__isl_keep
struct isl_upoly
*up
)
298 struct isl_upoly_rec
*rec
;
299 struct isl_upoly_rec
*dup
;
301 rec
= isl_upoly_as_rec(up
);
305 dup
= isl_upoly_alloc_rec(up
->ctx
, up
->var
, rec
->n
);
309 for (i
= 0; i
< rec
->n
; ++i
) {
310 dup
->p
[i
] = isl_upoly_copy(rec
->p
[i
]);
318 isl_upoly_free(&dup
->up
);
322 __isl_give
struct isl_upoly
*isl_upoly_dup(__isl_keep
struct isl_upoly
*up
)
324 struct isl_upoly
*dup
;
329 if (isl_upoly_is_cst(up
))
330 return isl_upoly_dup_cst(up
);
332 return isl_upoly_dup_rec(up
);
335 __isl_give
struct isl_upoly
*isl_upoly_cow(__isl_take
struct isl_upoly
*up
)
343 return isl_upoly_dup(up
);
346 void isl_upoly_free(__isl_take
struct isl_upoly
*up
)
355 upoly_free_cst((struct isl_upoly_cst
*)up
);
357 upoly_free_rec((struct isl_upoly_rec
*)up
);
359 isl_ctx_deref(up
->ctx
);
363 static void isl_upoly_cst_reduce(__isl_keep
struct isl_upoly_cst
*cst
)
368 isl_int_gcd(gcd
, cst
->n
, cst
->d
);
369 if (!isl_int_is_zero(gcd
) && !isl_int_is_one(gcd
)) {
370 isl_int_divexact(cst
->n
, cst
->n
, gcd
);
371 isl_int_divexact(cst
->d
, cst
->d
, gcd
);
376 __isl_give
struct isl_upoly
*isl_upoly_sum_cst(__isl_take
struct isl_upoly
*up1
,
377 __isl_take
struct isl_upoly
*up2
)
379 struct isl_upoly_cst
*cst1
;
380 struct isl_upoly_cst
*cst2
;
382 up1
= isl_upoly_cow(up1
);
386 cst1
= isl_upoly_as_cst(up1
);
387 cst2
= isl_upoly_as_cst(up2
);
389 if (isl_int_eq(cst1
->d
, cst2
->d
))
390 isl_int_add(cst1
->n
, cst1
->n
, cst2
->n
);
392 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->d
);
393 isl_int_addmul(cst1
->n
, cst2
->n
, cst1
->d
);
394 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
397 isl_upoly_cst_reduce(cst1
);
407 static __isl_give
struct isl_upoly
*replace_by_zero(
408 __isl_take
struct isl_upoly
*up
)
416 return isl_upoly_zero(ctx
);
419 __isl_give
struct isl_upoly
*isl_upoly_sum(__isl_take
struct isl_upoly
*up1
,
420 __isl_take
struct isl_upoly
*up2
)
423 struct isl_upoly_rec
*rec1
, *rec2
;
428 if (isl_upoly_is_nan(up1
)) {
433 if (isl_upoly_is_nan(up2
)) {
438 if (isl_upoly_is_zero(up1
)) {
443 if (isl_upoly_is_zero(up2
)) {
448 if (up1
->var
< up2
->var
)
449 return isl_upoly_sum(up2
, up1
);
451 if (up2
->var
< up1
->var
) {
452 struct isl_upoly_rec
*rec
;
453 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
457 up1
= isl_upoly_cow(up1
);
458 rec
= isl_upoly_as_rec(up1
);
461 rec
->p
[0] = isl_upoly_sum(rec
->p
[0], up2
);
462 if (rec
->n
== 1 && isl_upoly_is_zero(rec
->p
[0]))
463 up1
= replace_by_zero(up1
);
467 if (isl_upoly_is_cst(up1
))
468 return isl_upoly_sum_cst(up1
, up2
);
470 rec1
= isl_upoly_as_rec(up1
);
471 rec2
= isl_upoly_as_rec(up2
);
475 if (rec1
->n
< rec2
->n
)
476 return isl_upoly_sum(up2
, up1
);
478 up1
= isl_upoly_cow(up1
);
479 rec1
= isl_upoly_as_rec(up1
);
483 for (i
= rec2
->n
- 1; i
>= 0; --i
) {
484 rec1
->p
[i
] = isl_upoly_sum(rec1
->p
[i
],
485 isl_upoly_copy(rec2
->p
[i
]));
488 if (i
== rec1
->n
- 1 && isl_upoly_is_zero(rec1
->p
[i
])) {
489 isl_upoly_free(rec1
->p
[i
]);
495 up1
= replace_by_zero(up1
);
506 __isl_give
struct isl_upoly
*isl_upoly_neg_cst(__isl_take
struct isl_upoly
*up
)
508 struct isl_upoly_cst
*cst
;
510 if (isl_upoly_is_zero(up
))
513 up
= isl_upoly_cow(up
);
517 cst
= isl_upoly_as_cst(up
);
519 isl_int_neg(cst
->n
, cst
->n
);
524 __isl_give
struct isl_upoly
*isl_upoly_neg(__isl_take
struct isl_upoly
*up
)
527 struct isl_upoly_rec
*rec
;
532 if (isl_upoly_is_cst(up
))
533 return isl_upoly_neg_cst(up
);
535 up
= isl_upoly_cow(up
);
536 rec
= isl_upoly_as_rec(up
);
540 for (i
= 0; i
< rec
->n
; ++i
) {
541 rec
->p
[i
] = isl_upoly_neg(rec
->p
[i
]);
552 __isl_give
struct isl_upoly
*isl_upoly_mul_cst(__isl_take
struct isl_upoly
*up1
,
553 __isl_take
struct isl_upoly
*up2
)
555 struct isl_upoly_cst
*cst1
;
556 struct isl_upoly_cst
*cst2
;
558 up1
= isl_upoly_cow(up1
);
562 cst1
= isl_upoly_as_cst(up1
);
563 cst2
= isl_upoly_as_cst(up2
);
565 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->n
);
566 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
568 isl_upoly_cst_reduce(cst1
);
578 __isl_give
struct isl_upoly
*isl_upoly_mul_rec(__isl_take
struct isl_upoly
*up1
,
579 __isl_take
struct isl_upoly
*up2
)
581 struct isl_upoly_rec
*rec1
;
582 struct isl_upoly_rec
*rec2
;
583 struct isl_upoly_rec
*res
;
587 rec1
= isl_upoly_as_rec(up1
);
588 rec2
= isl_upoly_as_rec(up2
);
591 size
= rec1
->n
+ rec2
->n
- 1;
592 res
= isl_upoly_alloc_rec(up1
->ctx
, up1
->var
, size
);
596 for (i
= 0; i
< rec1
->n
; ++i
) {
597 res
->p
[i
] = isl_upoly_mul(isl_upoly_copy(rec2
->p
[0]),
598 isl_upoly_copy(rec1
->p
[i
]));
603 for (; i
< size
; ++i
) {
604 res
->p
[i
] = isl_upoly_zero(up1
->ctx
);
609 for (i
= 0; i
< rec1
->n
; ++i
) {
610 for (j
= 1; j
< rec2
->n
; ++j
) {
611 struct isl_upoly
*up
;
612 up
= isl_upoly_mul(isl_upoly_copy(rec2
->p
[j
]),
613 isl_upoly_copy(rec1
->p
[i
]));
614 res
->p
[i
+ j
] = isl_upoly_sum(res
->p
[i
+ j
], up
);
627 isl_upoly_free(&res
->up
);
631 __isl_give
struct isl_upoly
*isl_upoly_mul(__isl_take
struct isl_upoly
*up1
,
632 __isl_take
struct isl_upoly
*up2
)
637 if (isl_upoly_is_nan(up1
)) {
642 if (isl_upoly_is_nan(up2
)) {
647 if (isl_upoly_is_zero(up1
)) {
652 if (isl_upoly_is_zero(up2
)) {
657 if (isl_upoly_is_one(up1
)) {
662 if (isl_upoly_is_one(up2
)) {
667 if (up1
->var
< up2
->var
)
668 return isl_upoly_mul(up2
, up1
);
670 if (up2
->var
< up1
->var
) {
672 struct isl_upoly_rec
*rec
;
673 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
674 isl_ctx
*ctx
= up1
->ctx
;
677 return isl_upoly_nan(ctx
);
679 up1
= isl_upoly_cow(up1
);
680 rec
= isl_upoly_as_rec(up1
);
684 for (i
= 0; i
< rec
->n
; ++i
) {
685 rec
->p
[i
] = isl_upoly_mul(rec
->p
[i
],
686 isl_upoly_copy(up2
));
694 if (isl_upoly_is_cst(up1
))
695 return isl_upoly_mul_cst(up1
, up2
);
697 return isl_upoly_mul_rec(up1
, up2
);
704 __isl_give isl_qpolynomial
*isl_qpolynomial_alloc(__isl_take isl_dim
*dim
,
707 struct isl_qpolynomial
*qp
;
713 total
= isl_dim_total(dim
);
715 qp
= isl_calloc_type(dim
->ctx
, struct isl_qpolynomial
);
720 qp
->div
= isl_mat_alloc(dim
->ctx
, n_div
, 1 + 1 + total
+ n_div
);
729 isl_qpolynomial_free(qp
);
733 __isl_give isl_qpolynomial
*isl_qpolynomial_copy(__isl_keep isl_qpolynomial
*qp
)
742 __isl_give isl_qpolynomial
*isl_qpolynomial_dup(__isl_keep isl_qpolynomial
*qp
)
744 struct isl_qpolynomial
*dup
;
749 dup
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
);
752 isl_mat_free(dup
->div
);
753 dup
->div
= isl_mat_copy(qp
->div
);
756 dup
->upoly
= isl_upoly_copy(qp
->upoly
);
762 isl_qpolynomial_free(dup
);
766 __isl_give isl_qpolynomial
*isl_qpolynomial_cow(__isl_take isl_qpolynomial
*qp
)
774 return isl_qpolynomial_dup(qp
);
777 void isl_qpolynomial_free(__isl_take isl_qpolynomial
*qp
)
785 isl_dim_free(qp
->dim
);
786 isl_mat_free(qp
->div
);
787 isl_upoly_free(qp
->upoly
);
792 static int compatible_divs(__isl_keep isl_mat
*div1
, __isl_keep isl_mat
*div2
)
797 isl_assert(div1
->ctx
, div1
->n_row
>= div2
->n_row
&&
798 div1
->n_col
>= div2
->n_col
, return -1);
800 if (div1
->n_row
== div2
->n_row
)
801 return isl_mat_is_equal(div1
, div2
);
805 div1
->n_row
= div2
->n_row
;
806 div1
->n_col
= div2
->n_col
;
808 equal
= isl_mat_is_equal(div1
, div2
);
816 static void expand_row(__isl_keep isl_mat
*dst
, int d
,
817 __isl_keep isl_mat
*src
, int s
, int *exp
)
820 unsigned c
= src
->n_col
- src
->n_row
;
822 isl_seq_cpy(dst
->row
[d
], src
->row
[s
], c
);
823 isl_seq_clr(dst
->row
[d
] + c
, dst
->n_col
- c
);
825 for (i
= 0; i
< s
; ++i
)
826 isl_int_set(dst
->row
[d
][c
+ exp
[i
]], src
->row
[s
][c
+ i
]);
829 static int cmp_row(__isl_keep isl_mat
*div
, int i
, int j
)
833 li
= isl_seq_last_non_zero(div
->row
[i
], div
->n_col
);
834 lj
= isl_seq_last_non_zero(div
->row
[j
], div
->n_col
);
839 return isl_seq_cmp(div
->row
[i
], div
->row
[j
], div
->n_col
);
842 static __isl_give isl_mat
*merge_divs(__isl_keep isl_mat
*div1
,
843 __isl_keep isl_mat
*div2
, int *exp1
, int *exp2
)
847 unsigned d
= div1
->n_col
- div1
->n_row
;
849 div
= isl_mat_alloc(div1
->ctx
, 1 + div1
->n_row
+ div2
->n_row
,
850 d
+ div1
->n_row
+ div2
->n_row
);
854 for (i
= 0, j
= 0, k
= 0; i
< div1
->n_row
&& j
< div2
->n_row
; ++k
) {
857 expand_row(div
, k
, div1
, i
, exp1
);
858 expand_row(div
, k
+ 1, div2
, j
, exp2
);
860 cmp
= cmp_row(div
, k
, k
+ 1);
864 } else if (cmp
< 0) {
868 isl_seq_cpy(div
->row
[k
], div
->row
[k
+ 1], div
->n_col
);
871 for (; i
< div1
->n_row
; ++i
, ++k
) {
872 expand_row(div
, k
, div1
, i
, exp1
);
875 for (; j
< div2
->n_row
; ++j
, ++k
) {
876 expand_row(div
, k
, div2
, j
, exp2
);
886 static __isl_give
struct isl_upoly
*expand(__isl_take
struct isl_upoly
*up
,
890 struct isl_upoly_rec
*rec
;
892 if (isl_upoly_is_cst(up
))
898 if (exp
[up
->var
- first
] == up
->var
- first
)
901 up
= isl_upoly_cow(up
);
905 up
->var
= exp
[up
->var
- first
] + first
;
907 rec
= isl_upoly_as_rec(up
);
911 for (i
= 0; i
< rec
->n
; ++i
) {
912 rec
->p
[i
] = expand(rec
->p
[i
], exp
, first
);
923 static __isl_give isl_qpolynomial
*with_merged_divs(
924 __isl_give isl_qpolynomial
*(*fn
)(__isl_take isl_qpolynomial
*qp1
,
925 __isl_take isl_qpolynomial
*qp2
),
926 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
932 qp1
= isl_qpolynomial_cow(qp1
);
933 qp2
= isl_qpolynomial_cow(qp2
);
938 isl_assert(qp1
->div
->ctx
, qp1
->div
->n_row
>= qp2
->div
->n_row
&&
939 qp1
->div
->n_col
>= qp2
->div
->n_col
, goto error
);
941 exp1
= isl_alloc_array(qp1
->div
->ctx
, int, qp1
->div
->n_row
);
942 exp2
= isl_alloc_array(qp2
->div
->ctx
, int, qp2
->div
->n_row
);
946 div
= merge_divs(qp1
->div
, qp2
->div
, exp1
, exp2
);
950 isl_mat_free(qp1
->div
);
951 qp1
->div
= isl_mat_copy(div
);
952 isl_mat_free(qp2
->div
);
953 qp2
->div
= isl_mat_copy(div
);
955 qp1
->upoly
= expand(qp1
->upoly
, exp1
, div
->n_col
- div
->n_row
- 2);
956 qp2
->upoly
= expand(qp2
->upoly
, exp2
, div
->n_col
- div
->n_row
- 2);
958 if (!qp1
->upoly
|| !qp2
->upoly
)
970 isl_qpolynomial_free(qp1
);
971 isl_qpolynomial_free(qp2
);
975 __isl_give isl_qpolynomial
*isl_qpolynomial_add(__isl_take isl_qpolynomial
*qp1
,
976 __isl_take isl_qpolynomial
*qp2
)
978 qp1
= isl_qpolynomial_cow(qp1
);
983 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
984 return isl_qpolynomial_add(qp2
, qp1
);
986 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
987 if (!compatible_divs(qp1
->div
, qp2
->div
))
988 return with_merged_divs(isl_qpolynomial_add
, qp1
, qp2
);
990 qp1
->upoly
= isl_upoly_sum(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
994 isl_qpolynomial_free(qp2
);
998 isl_qpolynomial_free(qp1
);
999 isl_qpolynomial_free(qp2
);
1003 __isl_give isl_qpolynomial
*isl_qpolynomial_sub(__isl_take isl_qpolynomial
*qp1
,
1004 __isl_take isl_qpolynomial
*qp2
)
1006 return isl_qpolynomial_add(qp1
, isl_qpolynomial_neg(qp2
));
1009 __isl_give isl_qpolynomial
*isl_qpolynomial_neg(__isl_take isl_qpolynomial
*qp
)
1011 qp
= isl_qpolynomial_cow(qp
);
1016 qp
->upoly
= isl_upoly_neg(qp
->upoly
);
1022 isl_qpolynomial_free(qp
);
1026 __isl_give isl_qpolynomial
*isl_qpolynomial_mul(__isl_take isl_qpolynomial
*qp1
,
1027 __isl_take isl_qpolynomial
*qp2
)
1029 qp1
= isl_qpolynomial_cow(qp1
);
1034 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1035 return isl_qpolynomial_mul(qp2
, qp1
);
1037 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1038 if (!compatible_divs(qp1
->div
, qp2
->div
))
1039 return with_merged_divs(isl_qpolynomial_mul
, qp1
, qp2
);
1041 qp1
->upoly
= isl_upoly_mul(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1045 isl_qpolynomial_free(qp2
);
1049 isl_qpolynomial_free(qp1
);
1050 isl_qpolynomial_free(qp2
);
1054 __isl_give isl_qpolynomial
*isl_qpolynomial_zero(__isl_take isl_dim
*dim
)
1056 struct isl_qpolynomial
*qp
;
1057 struct isl_upoly_cst
*cst
;
1059 qp
= isl_qpolynomial_alloc(dim
, 0);
1063 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1069 isl_qpolynomial_free(qp
);
1073 __isl_give isl_qpolynomial
*isl_qpolynomial_infty(__isl_take isl_dim
*dim
)
1075 struct isl_qpolynomial
*qp
;
1076 struct isl_upoly_cst
*cst
;
1078 qp
= isl_qpolynomial_alloc(dim
, 0);
1082 qp
->upoly
= isl_upoly_infty(dim
->ctx
);
1088 isl_qpolynomial_free(qp
);
1092 __isl_give isl_qpolynomial
*isl_qpolynomial_nan(__isl_take isl_dim
*dim
)
1094 struct isl_qpolynomial
*qp
;
1095 struct isl_upoly_cst
*cst
;
1097 qp
= isl_qpolynomial_alloc(dim
, 0);
1101 qp
->upoly
= isl_upoly_nan(dim
->ctx
);
1107 isl_qpolynomial_free(qp
);
1111 __isl_give isl_qpolynomial
*isl_qpolynomial_cst(__isl_take isl_dim
*dim
,
1114 struct isl_qpolynomial
*qp
;
1115 struct isl_upoly_cst
*cst
;
1117 qp
= isl_qpolynomial_alloc(dim
, 0);
1121 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1124 cst
= isl_upoly_as_cst(qp
->upoly
);
1125 isl_int_set(cst
->n
, v
);
1129 isl_qpolynomial_free(qp
);
1133 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial
*qp
,
1134 isl_int
*n
, isl_int
*d
)
1136 struct isl_upoly_cst
*cst
;
1141 if (!isl_upoly_is_cst(qp
->upoly
))
1144 cst
= isl_upoly_as_cst(qp
->upoly
);
1149 isl_int_set(*n
, cst
->n
);
1151 isl_int_set(*d
, cst
->d
);
1156 __isl_give isl_qpolynomial
*isl_qpolynomial_pow(__isl_take isl_dim
*dim
,
1160 struct isl_qpolynomial
*qp
;
1161 struct isl_upoly_rec
*rec
;
1162 struct isl_upoly_cst
*cst
;
1163 struct isl_ctx
*ctx
;
1170 qp
= isl_qpolynomial_alloc(dim
, 0);
1174 qp
->upoly
= &isl_upoly_alloc_rec(ctx
, pos
, 1 + power
)->up
;
1177 rec
= isl_upoly_as_rec(qp
->upoly
);
1178 for (i
= 0; i
< 1 + power
; ++i
) {
1179 rec
->p
[i
] = isl_upoly_zero(ctx
);
1184 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1185 isl_int_set_si(cst
->n
, 1);
1189 isl_qpolynomial_free(qp
);
1193 __isl_give isl_qpolynomial
*isl_qpolynomial_var(__isl_take isl_dim
*dim
,
1194 enum isl_dim_type type
, unsigned pos
)
1199 isl_assert(dim
->ctx
, isl_dim_size(dim
, isl_dim_in
) == 0, goto error
);
1200 isl_assert(dim
->ctx
, pos
< isl_dim_size(dim
, type
), goto error
);
1202 if (type
== isl_dim_set
)
1203 pos
+= isl_dim_size(dim
, isl_dim_param
);
1205 return isl_qpolynomial_pow(dim
, pos
, 1);
1211 __isl_give isl_qpolynomial
*isl_qpolynomial_div_pow(__isl_take isl_div
*div
,
1214 struct isl_qpolynomial
*qp
= NULL
;
1215 struct isl_upoly_rec
*rec
;
1216 struct isl_upoly_cst
*cst
;
1222 isl_assert(div
->ctx
, div
->bmap
->n_div
== 1, goto error
);
1224 qp
= isl_qpolynomial_alloc(isl_basic_map_get_dim(div
->bmap
), 1);
1228 isl_seq_cpy(qp
->div
->row
[0], div
->line
[0], qp
->div
->n_col
- 1);
1229 isl_int_set_si(qp
->div
->row
[0][qp
->div
->n_col
- 1], 0);
1231 pos
= isl_dim_total(qp
->dim
);
1232 qp
->upoly
= &isl_upoly_alloc_rec(div
->ctx
, pos
, 1 + power
)->up
;
1235 rec
= isl_upoly_as_rec(qp
->upoly
);
1236 for (i
= 0; i
< 1 + power
; ++i
) {
1237 rec
->p
[i
] = isl_upoly_zero(div
->ctx
);
1242 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1243 isl_int_set_si(cst
->n
, 1);
1249 isl_qpolynomial_free(qp
);
1254 __isl_give isl_qpolynomial
*isl_qpolynomial_div(__isl_take isl_div
*div
)
1256 return isl_qpolynomial_div_pow(div
, 1);
1259 __isl_give isl_qpolynomial
*isl_qpolynomial_rat_cst(__isl_take isl_dim
*dim
,
1260 const isl_int n
, const isl_int d
)
1262 struct isl_qpolynomial
*qp
;
1263 struct isl_upoly_cst
*cst
;
1265 qp
= isl_qpolynomial_alloc(dim
, 0);
1269 qp
->upoly
= isl_upoly_zero(dim
->ctx
);
1272 cst
= isl_upoly_as_cst(qp
->upoly
);
1273 isl_int_set(cst
->n
, n
);
1274 isl_int_set(cst
->d
, d
);
1278 isl_qpolynomial_free(qp
);
1282 static __isl_give isl_pw_qpolynomial
*pw_qpolynomial_alloc(__isl_take isl_dim
*dim
,
1285 struct isl_pw_qpolynomial
*pwqp
;
1289 isl_assert(dim
->ctx
, n
>= 0, return NULL
);
1290 pwqp
= isl_alloc(dim
->ctx
, struct isl_pw_qpolynomial
,
1291 sizeof(struct isl_pw_qpolynomial
) +
1292 (n
- 1) * sizeof(struct isl_pw_qpolynomial_piece
));
1306 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add_piece(
1307 __isl_take isl_pw_qpolynomial
*pwqp
,
1308 __isl_take isl_set
*set
, __isl_take isl_qpolynomial
*qp
)
1310 if (!pwqp
|| !set
|| !qp
)
1313 if (isl_set_fast_is_empty(set
) || isl_qpolynomial_is_zero(qp
)) {
1315 isl_qpolynomial_free(qp
);
1319 isl_assert(set
->ctx
, isl_dim_equal(pwqp
->dim
, qp
->dim
), goto error
);
1320 isl_assert(set
->ctx
, pwqp
->n
< pwqp
->size
, goto error
);
1322 pwqp
->p
[pwqp
->n
].set
= set
;
1323 pwqp
->p
[pwqp
->n
].qp
= qp
;
1328 isl_pw_qpolynomial_free(pwqp
);
1330 isl_qpolynomial_free(qp
);
1334 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_alloc(__isl_take isl_set
*set
,
1335 __isl_take isl_qpolynomial
*qp
)
1337 struct isl_pw_qpolynomial
*pwqp
;
1342 pwqp
= pw_qpolynomial_alloc(isl_set_get_dim(set
), 1);
1344 return isl_pw_qpolynomial_add_piece(pwqp
, set
, qp
);
1347 isl_qpolynomial_free(qp
);
1351 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_zero(__isl_take isl_dim
*dim
)
1353 return pw_qpolynomial_alloc(dim
, 0);
1356 int isl_pw_qpolynomial_is_zero(__isl_keep isl_pw_qpolynomial
*pwqp
)
1361 return pwqp
->n
== 0;
1364 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial
*pwqp
)
1372 if (!isl_set_fast_is_universe(pwqp
->p
[0].set
))
1375 return isl_qpolynomial_is_one(pwqp
->p
[0].qp
);
1378 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_copy(
1379 __isl_keep isl_pw_qpolynomial
*pwqp
)
1388 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_dup(
1389 __isl_keep isl_pw_qpolynomial
*pwqp
)
1392 struct isl_pw_qpolynomial
*dup
;
1397 dup
= pw_qpolynomial_alloc(isl_dim_copy(pwqp
->dim
), pwqp
->n
);
1401 for (i
= 0; i
< pwqp
->n
; ++i
)
1402 dup
= isl_pw_qpolynomial_add_piece(dup
,
1403 isl_set_copy(pwqp
->p
[i
].set
),
1404 isl_qpolynomial_copy(pwqp
->p
[i
].qp
));
1408 isl_pw_qpolynomial_free(dup
);
1412 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_cow(
1413 __isl_take isl_pw_qpolynomial
*pwqp
)
1421 return isl_pw_qpolynomial_dup(pwqp
);
1424 void isl_pw_qpolynomial_free(__isl_take isl_pw_qpolynomial
*pwqp
)
1430 if (--pwqp
->ref
> 0)
1433 for (i
= 0; i
< pwqp
->n
; ++i
) {
1434 isl_set_free(pwqp
->p
[i
].set
);
1435 isl_qpolynomial_free(pwqp
->p
[i
].qp
);
1437 isl_dim_free(pwqp
->dim
);
1441 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add(
1442 __isl_take isl_pw_qpolynomial
*pwqp1
,
1443 __isl_take isl_pw_qpolynomial
*pwqp2
)
1446 struct isl_pw_qpolynomial
*res
;
1449 if (!pwqp1
|| !pwqp2
)
1452 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1455 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1456 isl_pw_qpolynomial_free(pwqp1
);
1460 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1461 isl_pw_qpolynomial_free(pwqp2
);
1465 n
= (pwqp1
->n
+ 1) * (pwqp2
->n
+ 1);
1466 res
= pw_qpolynomial_alloc(isl_dim_copy(pwqp1
->dim
), n
);
1468 for (i
= 0; i
< pwqp1
->n
; ++i
) {
1469 set
= isl_set_copy(pwqp1
->p
[i
].set
);
1470 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1471 struct isl_set
*common
;
1472 struct isl_qpolynomial
*sum
;
1473 set
= isl_set_subtract(set
,
1474 isl_set_copy(pwqp2
->p
[j
].set
));
1475 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
1476 isl_set_copy(pwqp2
->p
[j
].set
));
1477 if (isl_set_fast_is_empty(common
)) {
1478 isl_set_free(common
);
1482 sum
= isl_qpolynomial_add(
1483 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
1484 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1486 res
= isl_pw_qpolynomial_add_piece(res
, common
, sum
);
1488 res
= isl_pw_qpolynomial_add_piece(res
, set
,
1489 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
));
1492 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1493 set
= isl_set_copy(pwqp2
->p
[j
].set
);
1494 for (i
= 0; i
< pwqp1
->n
; ++i
)
1495 set
= isl_set_subtract(set
,
1496 isl_set_copy(pwqp1
->p
[i
].set
));
1497 res
= isl_pw_qpolynomial_add_piece(res
, set
,
1498 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1501 isl_pw_qpolynomial_free(pwqp1
);
1502 isl_pw_qpolynomial_free(pwqp2
);
1506 isl_pw_qpolynomial_free(pwqp1
);
1507 isl_pw_qpolynomial_free(pwqp2
);
1511 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add_disjoint(
1512 __isl_take isl_pw_qpolynomial
*pwqp1
,
1513 __isl_take isl_pw_qpolynomial
*pwqp2
)
1516 struct isl_pw_qpolynomial
*res
;
1518 if (!pwqp1
|| !pwqp2
)
1521 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1524 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1525 isl_pw_qpolynomial_free(pwqp1
);
1529 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1530 isl_pw_qpolynomial_free(pwqp2
);
1534 res
= pw_qpolynomial_alloc(isl_dim_copy(pwqp1
->dim
), pwqp1
->n
+ pwqp2
->n
);
1536 for (i
= 0; i
< pwqp1
->n
; ++i
)
1537 res
= isl_pw_qpolynomial_add_piece(res
,
1538 isl_set_copy(pwqp1
->p
[i
].set
),
1539 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
));
1541 for (i
= 0; i
< pwqp2
->n
; ++i
)
1542 res
= isl_pw_qpolynomial_add_piece(res
,
1543 isl_set_copy(pwqp2
->p
[i
].set
),
1544 isl_qpolynomial_copy(pwqp2
->p
[i
].qp
));
1546 isl_pw_qpolynomial_free(pwqp1
);
1547 isl_pw_qpolynomial_free(pwqp2
);
1551 isl_pw_qpolynomial_free(pwqp1
);
1552 isl_pw_qpolynomial_free(pwqp2
);
1556 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_mul(
1557 __isl_take isl_pw_qpolynomial
*pwqp1
,
1558 __isl_take isl_pw_qpolynomial
*pwqp2
)
1561 struct isl_pw_qpolynomial
*res
;
1564 if (!pwqp1
|| !pwqp2
)
1567 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1570 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1571 isl_pw_qpolynomial_free(pwqp2
);
1575 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1576 isl_pw_qpolynomial_free(pwqp1
);
1580 if (isl_pw_qpolynomial_is_one(pwqp1
)) {
1581 isl_pw_qpolynomial_free(pwqp1
);
1585 if (isl_pw_qpolynomial_is_one(pwqp2
)) {
1586 isl_pw_qpolynomial_free(pwqp2
);
1590 n
= pwqp1
->n
* pwqp2
->n
;
1591 res
= pw_qpolynomial_alloc(isl_dim_copy(pwqp1
->dim
), n
);
1593 for (i
= 0; i
< pwqp1
->n
; ++i
) {
1594 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1595 struct isl_set
*common
;
1596 struct isl_qpolynomial
*prod
;
1597 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
1598 isl_set_copy(pwqp2
->p
[j
].set
));
1599 if (isl_set_fast_is_empty(common
)) {
1600 isl_set_free(common
);
1604 prod
= isl_qpolynomial_mul(
1605 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
1606 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1608 res
= isl_pw_qpolynomial_add_piece(res
, common
, prod
);
1612 isl_pw_qpolynomial_free(pwqp1
);
1613 isl_pw_qpolynomial_free(pwqp2
);
1617 isl_pw_qpolynomial_free(pwqp1
);
1618 isl_pw_qpolynomial_free(pwqp2
);
1622 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_neg(
1623 __isl_take isl_pw_qpolynomial
*pwqp
)
1626 struct isl_pw_qpolynomial
*res
;
1632 if (isl_pw_qpolynomial_is_zero(pwqp
))
1635 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
1637 for (i
= 0; i
< pwqp
->n
; ++i
) {
1638 pwqp
->p
[i
].qp
= isl_qpolynomial_neg(pwqp
->p
[i
].qp
);
1645 isl_pw_qpolynomial_free(pwqp
);
1649 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_sub(
1650 __isl_take isl_pw_qpolynomial
*pwqp1
,
1651 __isl_take isl_pw_qpolynomial
*pwqp2
)
1653 return isl_pw_qpolynomial_add(pwqp1
, isl_pw_qpolynomial_neg(pwqp2
));
1656 __isl_give
struct isl_upoly
*isl_upoly_eval(
1657 __isl_take
struct isl_upoly
*up
, __isl_take isl_vec
*vec
)
1660 struct isl_upoly_rec
*rec
;
1661 struct isl_upoly
*res
;
1662 struct isl_upoly
*base
;
1664 if (isl_upoly_is_cst(up
)) {
1669 rec
= isl_upoly_as_rec(up
);
1673 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1675 base
= isl_upoly_rat_cst(up
->ctx
, vec
->el
[1 + up
->var
], vec
->el
[0]);
1677 res
= isl_upoly_eval(isl_upoly_copy(rec
->p
[rec
->n
- 1]),
1680 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1681 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1682 res
= isl_upoly_sum(res
,
1683 isl_upoly_eval(isl_upoly_copy(rec
->p
[i
]),
1684 isl_vec_copy(vec
)));
1687 isl_upoly_free(base
);
1697 __isl_give isl_qpolynomial
*isl_qpolynomial_eval(
1698 __isl_take isl_qpolynomial
*qp
, __isl_take isl_point
*pnt
)
1701 struct isl_upoly
*up
;
1706 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, qp
->dim
), goto error
);
1708 if (qp
->div
->n_row
== 0)
1709 ext
= isl_vec_copy(pnt
->vec
);
1712 unsigned dim
= isl_dim_total(qp
->dim
);
1713 ext
= isl_vec_alloc(qp
->dim
->ctx
, 1 + dim
+ qp
->div
->n_row
);
1717 isl_seq_cpy(ext
->el
, pnt
->vec
->el
, pnt
->vec
->size
);
1718 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1719 isl_seq_inner_product(qp
->div
->row
[i
] + 1, ext
->el
,
1720 1 + dim
+ i
, &ext
->el
[1+dim
+i
]);
1721 isl_int_fdiv_q(ext
->el
[1+dim
+i
], ext
->el
[1+dim
+i
],
1722 qp
->div
->row
[i
][0]);
1726 up
= isl_upoly_eval(isl_upoly_copy(qp
->upoly
), ext
);
1730 dim
= isl_dim_copy(qp
->dim
);
1731 isl_qpolynomial_free(qp
);
1732 isl_point_free(pnt
);
1734 qp
= isl_qpolynomial_alloc(dim
, 0);
1742 isl_qpolynomial_free(qp
);
1743 isl_point_free(pnt
);
1747 __isl_give isl_qpolynomial
*isl_pw_qpolynomial_eval(
1748 __isl_take isl_pw_qpolynomial
*pwqp
, __isl_take isl_point
*pnt
)
1752 isl_qpolynomial
*qp
;
1756 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, pwqp
->dim
), goto error
);
1758 for (i
= 0; i
< pwqp
->n
; ++i
) {
1759 found
= isl_set_contains_point(pwqp
->p
[i
].set
, pnt
);
1766 qp
= isl_qpolynomial_eval(isl_qpolynomial_copy(pwqp
->p
[i
].qp
),
1767 isl_point_copy(pnt
));
1769 qp
= isl_qpolynomial_zero(isl_dim_copy(pwqp
->dim
));
1770 isl_pw_qpolynomial_free(pwqp
);
1771 isl_point_free(pnt
);
1774 isl_pw_qpolynomial_free(pwqp
);
1775 isl_point_free(pnt
);
1779 __isl_give isl_qpolynomial
*isl_qpolynomial_move(__isl_take isl_qpolynomial
*qp
,
1780 enum isl_dim_type dst_type
, unsigned dst_pos
,
1781 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
1786 qp
->dim
= isl_dim_move(qp
->dim
, dst_type
, dst_pos
, src_type
, src_pos
, n
);
1790 /* Do something to polynomials when needed; later */
1794 isl_qpolynomial_free(qp
);
1798 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_move(
1799 __isl_take isl_pw_qpolynomial
*pwqp
,
1800 enum isl_dim_type dst_type
, unsigned dst_pos
,
1801 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
1808 pwqp
->dim
= isl_dim_move(pwqp
->dim
,
1809 dst_type
, dst_pos
, src_type
, src_pos
, n
);
1813 for (i
= 0; i
< pwqp
->n
; ++i
) {
1814 pwqp
->p
[i
].set
= isl_set_move(pwqp
->p
[i
].set
, dst_type
, dst_pos
,
1815 src_type
, src_pos
, n
);
1816 if (!pwqp
->p
[i
].set
)
1818 pwqp
->p
[i
].qp
= isl_qpolynomial_move(pwqp
->p
[i
].qp
,
1819 dst_type
, dst_pos
, src_type
, src_pos
, n
);
1826 isl_pw_qpolynomial_free(pwqp
);