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_give isl_qpolynomial
*isl_qpolynomial_reset_dim(
323 __isl_take isl_qpolynomial
*qp
, __isl_take isl_dim
*dim
)
328 isl_dim_free(qp
->dim
);
333 isl_qpolynomial_free(qp
);
338 isl_ctx
*isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial
*qp
)
340 return qp
? qp
->dim
->ctx
: NULL
;
343 __isl_give isl_dim
*isl_qpolynomial_get_dim(__isl_keep isl_qpolynomial
*qp
)
345 return qp
? isl_dim_copy(qp
->dim
) : NULL
;
348 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial
*qp
,
349 enum isl_dim_type type
)
351 return qp
? isl_dim_size(qp
->dim
, type
) : 0;
354 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial
*qp
)
356 return qp
? isl_upoly_is_zero(qp
->upoly
) : -1;
359 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial
*qp
)
361 return qp
? isl_upoly_is_one(qp
->upoly
) : -1;
364 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial
*qp
)
366 return qp
? isl_upoly_is_nan(qp
->upoly
) : -1;
369 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial
*qp
)
371 return qp
? isl_upoly_is_infty(qp
->upoly
) : -1;
374 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial
*qp
)
376 return qp
? isl_upoly_is_neginfty(qp
->upoly
) : -1;
379 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial
*qp
)
381 return qp
? isl_upoly_sgn(qp
->upoly
) : 0;
384 static void upoly_free_cst(__isl_take
struct isl_upoly_cst
*cst
)
386 isl_int_clear(cst
->n
);
387 isl_int_clear(cst
->d
);
390 static void upoly_free_rec(__isl_take
struct isl_upoly_rec
*rec
)
394 for (i
= 0; i
< rec
->n
; ++i
)
395 isl_upoly_free(rec
->p
[i
]);
398 __isl_give
struct isl_upoly
*isl_upoly_copy(__isl_keep
struct isl_upoly
*up
)
407 __isl_give
struct isl_upoly
*isl_upoly_dup_cst(__isl_keep
struct isl_upoly
*up
)
409 struct isl_upoly_cst
*cst
;
410 struct isl_upoly_cst
*dup
;
412 cst
= isl_upoly_as_cst(up
);
416 dup
= isl_upoly_as_cst(isl_upoly_zero(up
->ctx
));
419 isl_int_set(dup
->n
, cst
->n
);
420 isl_int_set(dup
->d
, cst
->d
);
425 __isl_give
struct isl_upoly
*isl_upoly_dup_rec(__isl_keep
struct isl_upoly
*up
)
428 struct isl_upoly_rec
*rec
;
429 struct isl_upoly_rec
*dup
;
431 rec
= isl_upoly_as_rec(up
);
435 dup
= isl_upoly_alloc_rec(up
->ctx
, up
->var
, rec
->n
);
439 for (i
= 0; i
< rec
->n
; ++i
) {
440 dup
->p
[i
] = isl_upoly_copy(rec
->p
[i
]);
448 isl_upoly_free(&dup
->up
);
452 __isl_give
struct isl_upoly
*isl_upoly_dup(__isl_keep
struct isl_upoly
*up
)
454 struct isl_upoly
*dup
;
459 if (isl_upoly_is_cst(up
))
460 return isl_upoly_dup_cst(up
);
462 return isl_upoly_dup_rec(up
);
465 __isl_give
struct isl_upoly
*isl_upoly_cow(__isl_take
struct isl_upoly
*up
)
473 return isl_upoly_dup(up
);
476 void isl_upoly_free(__isl_take
struct isl_upoly
*up
)
485 upoly_free_cst((struct isl_upoly_cst
*)up
);
487 upoly_free_rec((struct isl_upoly_rec
*)up
);
489 isl_ctx_deref(up
->ctx
);
493 static void isl_upoly_cst_reduce(__isl_keep
struct isl_upoly_cst
*cst
)
498 isl_int_gcd(gcd
, cst
->n
, cst
->d
);
499 if (!isl_int_is_zero(gcd
) && !isl_int_is_one(gcd
)) {
500 isl_int_divexact(cst
->n
, cst
->n
, gcd
);
501 isl_int_divexact(cst
->d
, cst
->d
, gcd
);
506 __isl_give
struct isl_upoly
*isl_upoly_sum_cst(__isl_take
struct isl_upoly
*up1
,
507 __isl_take
struct isl_upoly
*up2
)
509 struct isl_upoly_cst
*cst1
;
510 struct isl_upoly_cst
*cst2
;
512 up1
= isl_upoly_cow(up1
);
516 cst1
= isl_upoly_as_cst(up1
);
517 cst2
= isl_upoly_as_cst(up2
);
519 if (isl_int_eq(cst1
->d
, cst2
->d
))
520 isl_int_add(cst1
->n
, cst1
->n
, cst2
->n
);
522 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->d
);
523 isl_int_addmul(cst1
->n
, cst2
->n
, cst1
->d
);
524 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
527 isl_upoly_cst_reduce(cst1
);
537 static __isl_give
struct isl_upoly
*replace_by_zero(
538 __isl_take
struct isl_upoly
*up
)
546 return isl_upoly_zero(ctx
);
549 static __isl_give
struct isl_upoly
*replace_by_constant_term(
550 __isl_take
struct isl_upoly
*up
)
552 struct isl_upoly_rec
*rec
;
553 struct isl_upoly
*cst
;
558 rec
= isl_upoly_as_rec(up
);
561 cst
= isl_upoly_copy(rec
->p
[0]);
569 __isl_give
struct isl_upoly
*isl_upoly_sum(__isl_take
struct isl_upoly
*up1
,
570 __isl_take
struct isl_upoly
*up2
)
573 struct isl_upoly_rec
*rec1
, *rec2
;
578 if (isl_upoly_is_nan(up1
)) {
583 if (isl_upoly_is_nan(up2
)) {
588 if (isl_upoly_is_zero(up1
)) {
593 if (isl_upoly_is_zero(up2
)) {
598 if (up1
->var
< up2
->var
)
599 return isl_upoly_sum(up2
, up1
);
601 if (up2
->var
< up1
->var
) {
602 struct isl_upoly_rec
*rec
;
603 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
607 up1
= isl_upoly_cow(up1
);
608 rec
= isl_upoly_as_rec(up1
);
611 rec
->p
[0] = isl_upoly_sum(rec
->p
[0], up2
);
613 up1
= replace_by_constant_term(up1
);
617 if (isl_upoly_is_cst(up1
))
618 return isl_upoly_sum_cst(up1
, up2
);
620 rec1
= isl_upoly_as_rec(up1
);
621 rec2
= isl_upoly_as_rec(up2
);
625 if (rec1
->n
< rec2
->n
)
626 return isl_upoly_sum(up2
, up1
);
628 up1
= isl_upoly_cow(up1
);
629 rec1
= isl_upoly_as_rec(up1
);
633 for (i
= rec2
->n
- 1; i
>= 0; --i
) {
634 rec1
->p
[i
] = isl_upoly_sum(rec1
->p
[i
],
635 isl_upoly_copy(rec2
->p
[i
]));
638 if (i
== rec1
->n
- 1 && isl_upoly_is_zero(rec1
->p
[i
])) {
639 isl_upoly_free(rec1
->p
[i
]);
645 up1
= replace_by_zero(up1
);
646 else if (rec1
->n
== 1)
647 up1
= replace_by_constant_term(up1
);
658 __isl_give
struct isl_upoly
*isl_upoly_neg_cst(__isl_take
struct isl_upoly
*up
)
660 struct isl_upoly_cst
*cst
;
662 if (isl_upoly_is_zero(up
))
665 up
= isl_upoly_cow(up
);
669 cst
= isl_upoly_as_cst(up
);
671 isl_int_neg(cst
->n
, cst
->n
);
676 __isl_give
struct isl_upoly
*isl_upoly_neg(__isl_take
struct isl_upoly
*up
)
679 struct isl_upoly_rec
*rec
;
684 if (isl_upoly_is_cst(up
))
685 return isl_upoly_neg_cst(up
);
687 up
= isl_upoly_cow(up
);
688 rec
= isl_upoly_as_rec(up
);
692 for (i
= 0; i
< rec
->n
; ++i
) {
693 rec
->p
[i
] = isl_upoly_neg(rec
->p
[i
]);
704 __isl_give
struct isl_upoly
*isl_upoly_mul_cst(__isl_take
struct isl_upoly
*up1
,
705 __isl_take
struct isl_upoly
*up2
)
707 struct isl_upoly_cst
*cst1
;
708 struct isl_upoly_cst
*cst2
;
710 up1
= isl_upoly_cow(up1
);
714 cst1
= isl_upoly_as_cst(up1
);
715 cst2
= isl_upoly_as_cst(up2
);
717 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->n
);
718 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
720 isl_upoly_cst_reduce(cst1
);
730 __isl_give
struct isl_upoly
*isl_upoly_mul_rec(__isl_take
struct isl_upoly
*up1
,
731 __isl_take
struct isl_upoly
*up2
)
733 struct isl_upoly_rec
*rec1
;
734 struct isl_upoly_rec
*rec2
;
735 struct isl_upoly_rec
*res
;
739 rec1
= isl_upoly_as_rec(up1
);
740 rec2
= isl_upoly_as_rec(up2
);
743 size
= rec1
->n
+ rec2
->n
- 1;
744 res
= isl_upoly_alloc_rec(up1
->ctx
, up1
->var
, size
);
748 for (i
= 0; i
< rec1
->n
; ++i
) {
749 res
->p
[i
] = isl_upoly_mul(isl_upoly_copy(rec2
->p
[0]),
750 isl_upoly_copy(rec1
->p
[i
]));
755 for (; i
< size
; ++i
) {
756 res
->p
[i
] = isl_upoly_zero(up1
->ctx
);
761 for (i
= 0; i
< rec1
->n
; ++i
) {
762 for (j
= 1; j
< rec2
->n
; ++j
) {
763 struct isl_upoly
*up
;
764 up
= isl_upoly_mul(isl_upoly_copy(rec2
->p
[j
]),
765 isl_upoly_copy(rec1
->p
[i
]));
766 res
->p
[i
+ j
] = isl_upoly_sum(res
->p
[i
+ j
], up
);
779 isl_upoly_free(&res
->up
);
783 __isl_give
struct isl_upoly
*isl_upoly_mul(__isl_take
struct isl_upoly
*up1
,
784 __isl_take
struct isl_upoly
*up2
)
789 if (isl_upoly_is_nan(up1
)) {
794 if (isl_upoly_is_nan(up2
)) {
799 if (isl_upoly_is_zero(up1
)) {
804 if (isl_upoly_is_zero(up2
)) {
809 if (isl_upoly_is_one(up1
)) {
814 if (isl_upoly_is_one(up2
)) {
819 if (up1
->var
< up2
->var
)
820 return isl_upoly_mul(up2
, up1
);
822 if (up2
->var
< up1
->var
) {
824 struct isl_upoly_rec
*rec
;
825 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
826 isl_ctx
*ctx
= up1
->ctx
;
829 return isl_upoly_nan(ctx
);
831 up1
= isl_upoly_cow(up1
);
832 rec
= isl_upoly_as_rec(up1
);
836 for (i
= 0; i
< rec
->n
; ++i
) {
837 rec
->p
[i
] = isl_upoly_mul(rec
->p
[i
],
838 isl_upoly_copy(up2
));
846 if (isl_upoly_is_cst(up1
))
847 return isl_upoly_mul_cst(up1
, up2
);
849 return isl_upoly_mul_rec(up1
, up2
);
856 __isl_give isl_qpolynomial
*isl_qpolynomial_alloc(__isl_take isl_dim
*dim
,
857 unsigned n_div
, __isl_take
struct isl_upoly
*up
)
859 struct isl_qpolynomial
*qp
= NULL
;
865 total
= isl_dim_total(dim
);
867 qp
= isl_calloc_type(dim
->ctx
, struct isl_qpolynomial
);
872 qp
->div
= isl_mat_alloc(dim
->ctx
, n_div
, 1 + 1 + total
+ n_div
);
883 isl_qpolynomial_free(qp
);
887 __isl_give isl_qpolynomial
*isl_qpolynomial_copy(__isl_keep isl_qpolynomial
*qp
)
896 __isl_give isl_qpolynomial
*isl_qpolynomial_dup(__isl_keep isl_qpolynomial
*qp
)
898 struct isl_qpolynomial
*dup
;
903 dup
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
,
904 isl_upoly_copy(qp
->upoly
));
907 isl_mat_free(dup
->div
);
908 dup
->div
= isl_mat_copy(qp
->div
);
914 isl_qpolynomial_free(dup
);
918 __isl_give isl_qpolynomial
*isl_qpolynomial_cow(__isl_take isl_qpolynomial
*qp
)
926 return isl_qpolynomial_dup(qp
);
929 void isl_qpolynomial_free(__isl_take isl_qpolynomial
*qp
)
937 isl_dim_free(qp
->dim
);
938 isl_mat_free(qp
->div
);
939 isl_upoly_free(qp
->upoly
);
944 __isl_give
struct isl_upoly
*isl_upoly_pow(isl_ctx
*ctx
, int pos
, int power
)
947 struct isl_upoly
*up
;
948 struct isl_upoly_rec
*rec
;
949 struct isl_upoly_cst
*cst
;
951 rec
= isl_upoly_alloc_rec(ctx
, pos
, 1 + power
);
954 for (i
= 0; i
< 1 + power
; ++i
) {
955 rec
->p
[i
] = isl_upoly_zero(ctx
);
960 cst
= isl_upoly_as_cst(rec
->p
[power
]);
961 isl_int_set_si(cst
->n
, 1);
965 isl_upoly_free(&rec
->up
);
969 /* r array maps original positions to new positions.
971 static __isl_give
struct isl_upoly
*reorder(__isl_take
struct isl_upoly
*up
,
975 struct isl_upoly_rec
*rec
;
976 struct isl_upoly
*base
;
977 struct isl_upoly
*res
;
979 if (isl_upoly_is_cst(up
))
982 rec
= isl_upoly_as_rec(up
);
986 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
988 base
= isl_upoly_pow(up
->ctx
, r
[up
->var
], 1);
989 res
= reorder(isl_upoly_copy(rec
->p
[rec
->n
- 1]), r
);
991 for (i
= rec
->n
- 2; i
>= 0; --i
) {
992 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
993 res
= isl_upoly_sum(res
, reorder(isl_upoly_copy(rec
->p
[i
]), r
));
996 isl_upoly_free(base
);
1005 static int compatible_divs(__isl_keep isl_mat
*div1
, __isl_keep isl_mat
*div2
)
1010 isl_assert(div1
->ctx
, div1
->n_row
>= div2
->n_row
&&
1011 div1
->n_col
>= div2
->n_col
, return -1);
1013 if (div1
->n_row
== div2
->n_row
)
1014 return isl_mat_is_equal(div1
, div2
);
1016 n_row
= div1
->n_row
;
1017 n_col
= div1
->n_col
;
1018 div1
->n_row
= div2
->n_row
;
1019 div1
->n_col
= div2
->n_col
;
1021 equal
= isl_mat_is_equal(div1
, div2
);
1023 div1
->n_row
= n_row
;
1024 div1
->n_col
= n_col
;
1029 static void expand_row(__isl_keep isl_mat
*dst
, int d
,
1030 __isl_keep isl_mat
*src
, int s
, int *exp
)
1033 unsigned c
= src
->n_col
- src
->n_row
;
1035 isl_seq_cpy(dst
->row
[d
], src
->row
[s
], c
);
1036 isl_seq_clr(dst
->row
[d
] + c
, dst
->n_col
- c
);
1038 for (i
= 0; i
< s
; ++i
)
1039 isl_int_set(dst
->row
[d
][c
+ exp
[i
]], src
->row
[s
][c
+ i
]);
1042 static int cmp_row(__isl_keep isl_mat
*div
, int i
, int j
)
1046 li
= isl_seq_last_non_zero(div
->row
[i
], div
->n_col
);
1047 lj
= isl_seq_last_non_zero(div
->row
[j
], div
->n_col
);
1052 return isl_seq_cmp(div
->row
[i
], div
->row
[j
], div
->n_col
);
1055 struct isl_div_sort_info
{
1060 static int div_sort_cmp(const void *p1
, const void *p2
)
1062 const struct isl_div_sort_info
*i1
, *i2
;
1063 i1
= (const struct isl_div_sort_info
*) p1
;
1064 i2
= (const struct isl_div_sort_info
*) p2
;
1066 return cmp_row(i1
->div
, i1
->row
, i2
->row
);
1069 static __isl_give isl_qpolynomial
*sort_divs(__isl_take isl_qpolynomial
*qp
)
1072 struct isl_div_sort_info
*array
= NULL
;
1074 int *reordering
= NULL
;
1079 if (qp
->div
->n_row
<= 1)
1082 div_pos
= isl_dim_total(qp
->dim
);
1084 array
= isl_alloc_array(qp
->div
->ctx
, struct isl_div_sort_info
,
1086 pos
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
1087 reordering
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_col
- 2);
1088 if (!array
|| !pos
|| !reordering
)
1091 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1092 array
[i
].div
= qp
->div
;
1097 qsort(array
, qp
->div
->n_row
, sizeof(struct isl_div_sort_info
),
1100 for (i
= 0; i
< div_pos
; ++i
)
1103 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
1104 reordering
[div_pos
+ array
[i
].row
] = div_pos
+ i
;
1106 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1108 if (pos
[array
[i
].row
] == i
)
1110 qp
->div
= isl_mat_cow(qp
->div
);
1111 qp
->div
= isl_mat_swap_rows(qp
->div
, i
, pos
[array
[i
].row
]);
1112 t
= pos
[array
[i
].row
];
1113 pos
[array
[i
].row
] = pos
[i
];
1117 qp
->upoly
= reorder(qp
->upoly
, reordering
);
1119 if (!qp
->upoly
|| !qp
->div
)
1131 isl_qpolynomial_free(qp
);
1135 static __isl_give isl_mat
*merge_divs(__isl_keep isl_mat
*div1
,
1136 __isl_keep isl_mat
*div2
, int *exp1
, int *exp2
)
1139 isl_mat
*div
= NULL
;
1140 unsigned d
= div1
->n_col
- div1
->n_row
;
1142 div
= isl_mat_alloc(div1
->ctx
, 1 + div1
->n_row
+ div2
->n_row
,
1143 d
+ div1
->n_row
+ div2
->n_row
);
1147 for (i
= 0, j
= 0, k
= 0; i
< div1
->n_row
&& j
< div2
->n_row
; ++k
) {
1150 expand_row(div
, k
, div1
, i
, exp1
);
1151 expand_row(div
, k
+ 1, div2
, j
, exp2
);
1153 cmp
= cmp_row(div
, k
, k
+ 1);
1157 } else if (cmp
< 0) {
1161 isl_seq_cpy(div
->row
[k
], div
->row
[k
+ 1], div
->n_col
);
1164 for (; i
< div1
->n_row
; ++i
, ++k
) {
1165 expand_row(div
, k
, div1
, i
, exp1
);
1168 for (; j
< div2
->n_row
; ++j
, ++k
) {
1169 expand_row(div
, k
, div2
, j
, exp2
);
1179 static __isl_give
struct isl_upoly
*expand(__isl_take
struct isl_upoly
*up
,
1180 int *exp
, int first
)
1183 struct isl_upoly_rec
*rec
;
1185 if (isl_upoly_is_cst(up
))
1188 if (up
->var
< first
)
1191 if (exp
[up
->var
- first
] == up
->var
- first
)
1194 up
= isl_upoly_cow(up
);
1198 up
->var
= exp
[up
->var
- first
] + first
;
1200 rec
= isl_upoly_as_rec(up
);
1204 for (i
= 0; i
< rec
->n
; ++i
) {
1205 rec
->p
[i
] = expand(rec
->p
[i
], exp
, first
);
1216 static __isl_give isl_qpolynomial
*with_merged_divs(
1217 __isl_give isl_qpolynomial
*(*fn
)(__isl_take isl_qpolynomial
*qp1
,
1218 __isl_take isl_qpolynomial
*qp2
),
1219 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1223 isl_mat
*div
= NULL
;
1225 qp1
= isl_qpolynomial_cow(qp1
);
1226 qp2
= isl_qpolynomial_cow(qp2
);
1231 isl_assert(qp1
->div
->ctx
, qp1
->div
->n_row
>= qp2
->div
->n_row
&&
1232 qp1
->div
->n_col
>= qp2
->div
->n_col
, goto error
);
1234 exp1
= isl_alloc_array(qp1
->div
->ctx
, int, qp1
->div
->n_row
);
1235 exp2
= isl_alloc_array(qp2
->div
->ctx
, int, qp2
->div
->n_row
);
1239 div
= merge_divs(qp1
->div
, qp2
->div
, exp1
, exp2
);
1243 isl_mat_free(qp1
->div
);
1244 qp1
->div
= isl_mat_copy(div
);
1245 isl_mat_free(qp2
->div
);
1246 qp2
->div
= isl_mat_copy(div
);
1248 qp1
->upoly
= expand(qp1
->upoly
, exp1
, div
->n_col
- div
->n_row
- 2);
1249 qp2
->upoly
= expand(qp2
->upoly
, exp2
, div
->n_col
- div
->n_row
- 2);
1251 if (!qp1
->upoly
|| !qp2
->upoly
)
1258 return fn(qp1
, qp2
);
1263 isl_qpolynomial_free(qp1
);
1264 isl_qpolynomial_free(qp2
);
1268 __isl_give isl_qpolynomial
*isl_qpolynomial_add(__isl_take isl_qpolynomial
*qp1
,
1269 __isl_take isl_qpolynomial
*qp2
)
1271 qp1
= isl_qpolynomial_cow(qp1
);
1276 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1277 return isl_qpolynomial_add(qp2
, qp1
);
1279 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1280 if (!compatible_divs(qp1
->div
, qp2
->div
))
1281 return with_merged_divs(isl_qpolynomial_add
, qp1
, qp2
);
1283 qp1
->upoly
= isl_upoly_sum(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1287 isl_qpolynomial_free(qp2
);
1291 isl_qpolynomial_free(qp1
);
1292 isl_qpolynomial_free(qp2
);
1296 __isl_give isl_qpolynomial
*isl_qpolynomial_sub(__isl_take isl_qpolynomial
*qp1
,
1297 __isl_take isl_qpolynomial
*qp2
)
1299 return isl_qpolynomial_add(qp1
, isl_qpolynomial_neg(qp2
));
1302 __isl_give isl_qpolynomial
*isl_qpolynomial_neg(__isl_take isl_qpolynomial
*qp
)
1304 qp
= isl_qpolynomial_cow(qp
);
1309 qp
->upoly
= isl_upoly_neg(qp
->upoly
);
1315 isl_qpolynomial_free(qp
);
1319 __isl_give isl_qpolynomial
*isl_qpolynomial_mul(__isl_take isl_qpolynomial
*qp1
,
1320 __isl_take isl_qpolynomial
*qp2
)
1322 qp1
= isl_qpolynomial_cow(qp1
);
1327 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1328 return isl_qpolynomial_mul(qp2
, qp1
);
1330 isl_assert(qp1
->dim
->ctx
, isl_dim_equal(qp1
->dim
, qp2
->dim
), goto error
);
1331 if (!compatible_divs(qp1
->div
, qp2
->div
))
1332 return with_merged_divs(isl_qpolynomial_mul
, qp1
, qp2
);
1334 qp1
->upoly
= isl_upoly_mul(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1338 isl_qpolynomial_free(qp2
);
1342 isl_qpolynomial_free(qp1
);
1343 isl_qpolynomial_free(qp2
);
1347 __isl_give isl_qpolynomial
*isl_qpolynomial_zero(__isl_take isl_dim
*dim
)
1349 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1352 __isl_give isl_qpolynomial
*isl_qpolynomial_infty(__isl_take isl_dim
*dim
)
1354 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_infty(dim
->ctx
));
1357 __isl_give isl_qpolynomial
*isl_qpolynomial_neginfty(__isl_take isl_dim
*dim
)
1359 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_neginfty(dim
->ctx
));
1362 __isl_give isl_qpolynomial
*isl_qpolynomial_nan(__isl_take isl_dim
*dim
)
1364 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_nan(dim
->ctx
));
1367 __isl_give isl_qpolynomial
*isl_qpolynomial_cst(__isl_take isl_dim
*dim
,
1370 struct isl_qpolynomial
*qp
;
1371 struct isl_upoly_cst
*cst
;
1373 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1377 cst
= isl_upoly_as_cst(qp
->upoly
);
1378 isl_int_set(cst
->n
, v
);
1383 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial
*qp
,
1384 isl_int
*n
, isl_int
*d
)
1386 struct isl_upoly_cst
*cst
;
1391 if (!isl_upoly_is_cst(qp
->upoly
))
1394 cst
= isl_upoly_as_cst(qp
->upoly
);
1399 isl_int_set(*n
, cst
->n
);
1401 isl_int_set(*d
, cst
->d
);
1406 int isl_upoly_is_affine(__isl_keep
struct isl_upoly
*up
)
1409 struct isl_upoly_rec
*rec
;
1417 rec
= isl_upoly_as_rec(up
);
1424 isl_assert(up
->ctx
, rec
->n
> 1, return -1);
1426 is_cst
= isl_upoly_is_cst(rec
->p
[1]);
1432 return isl_upoly_is_affine(rec
->p
[0]);
1435 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial
*qp
)
1440 if (qp
->div
->n_row
> 0)
1443 return isl_upoly_is_affine(qp
->upoly
);
1446 static void update_coeff(__isl_keep isl_vec
*aff
,
1447 __isl_keep
struct isl_upoly_cst
*cst
, int pos
)
1452 if (isl_int_is_zero(cst
->n
))
1457 isl_int_gcd(gcd
, cst
->d
, aff
->el
[0]);
1458 isl_int_divexact(f
, cst
->d
, gcd
);
1459 isl_int_divexact(gcd
, aff
->el
[0], gcd
);
1460 isl_seq_scale(aff
->el
, aff
->el
, f
, aff
->size
);
1461 isl_int_mul(aff
->el
[1 + pos
], gcd
, cst
->n
);
1466 int isl_upoly_update_affine(__isl_keep
struct isl_upoly
*up
,
1467 __isl_keep isl_vec
*aff
)
1469 struct isl_upoly_cst
*cst
;
1470 struct isl_upoly_rec
*rec
;
1476 struct isl_upoly_cst
*cst
;
1478 cst
= isl_upoly_as_cst(up
);
1481 update_coeff(aff
, cst
, 0);
1485 rec
= isl_upoly_as_rec(up
);
1488 isl_assert(up
->ctx
, rec
->n
== 2, return -1);
1490 cst
= isl_upoly_as_cst(rec
->p
[1]);
1493 update_coeff(aff
, cst
, 1 + up
->var
);
1495 return isl_upoly_update_affine(rec
->p
[0], aff
);
1498 __isl_give isl_vec
*isl_qpolynomial_extract_affine(
1499 __isl_keep isl_qpolynomial
*qp
)
1507 isl_assert(qp
->div
->ctx
, qp
->div
->n_row
== 0, return NULL
);
1508 d
= isl_dim_total(qp
->dim
);
1509 aff
= isl_vec_alloc(qp
->div
->ctx
, 2 + d
);
1513 isl_seq_clr(aff
->el
+ 1, 1 + d
);
1514 isl_int_set_si(aff
->el
[0], 1);
1516 if (isl_upoly_update_affine(qp
->upoly
, aff
) < 0)
1525 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial
*qp1
,
1526 __isl_keep isl_qpolynomial
*qp2
)
1531 return isl_upoly_is_equal(qp1
->upoly
, qp2
->upoly
);
1534 static void upoly_update_den(__isl_keep
struct isl_upoly
*up
, isl_int
*d
)
1537 struct isl_upoly_rec
*rec
;
1539 if (isl_upoly_is_cst(up
)) {
1540 struct isl_upoly_cst
*cst
;
1541 cst
= isl_upoly_as_cst(up
);
1544 isl_int_lcm(*d
, *d
, cst
->d
);
1548 rec
= isl_upoly_as_rec(up
);
1552 for (i
= 0; i
< rec
->n
; ++i
)
1553 upoly_update_den(rec
->p
[i
], d
);
1556 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial
*qp
, isl_int
*d
)
1558 isl_int_set_si(*d
, 1);
1561 upoly_update_den(qp
->upoly
, d
);
1564 __isl_give isl_qpolynomial
*isl_qpolynomial_pow(__isl_take isl_dim
*dim
,
1567 struct isl_ctx
*ctx
;
1574 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_pow(ctx
, pos
, power
));
1577 __isl_give isl_qpolynomial
*isl_qpolynomial_var(__isl_take isl_dim
*dim
,
1578 enum isl_dim_type type
, unsigned pos
)
1583 isl_assert(dim
->ctx
, isl_dim_size(dim
, isl_dim_in
) == 0, goto error
);
1584 isl_assert(dim
->ctx
, pos
< isl_dim_size(dim
, type
), goto error
);
1586 if (type
== isl_dim_set
)
1587 pos
+= isl_dim_size(dim
, isl_dim_param
);
1589 return isl_qpolynomial_pow(dim
, pos
, 1);
1595 __isl_give isl_qpolynomial
*isl_qpolynomial_div_pow(__isl_take isl_div
*div
,
1598 struct isl_qpolynomial
*qp
= NULL
;
1599 struct isl_upoly_rec
*rec
;
1600 struct isl_upoly_cst
*cst
;
1606 isl_assert(div
->ctx
, div
->bmap
->n_div
== 1, goto error
);
1608 pos
= isl_dim_total(div
->bmap
->dim
);
1609 rec
= isl_upoly_alloc_rec(div
->ctx
, pos
, 1 + power
);
1610 qp
= isl_qpolynomial_alloc(isl_basic_map_get_dim(div
->bmap
), 1,
1615 isl_seq_cpy(qp
->div
->row
[0], div
->line
[0], qp
->div
->n_col
- 1);
1616 isl_int_set_si(qp
->div
->row
[0][qp
->div
->n_col
- 1], 0);
1618 for (i
= 0; i
< 1 + power
; ++i
) {
1619 rec
->p
[i
] = isl_upoly_zero(div
->ctx
);
1624 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1625 isl_int_set_si(cst
->n
, 1);
1631 isl_qpolynomial_free(qp
);
1636 __isl_give isl_qpolynomial
*isl_qpolynomial_div(__isl_take isl_div
*div
)
1638 return isl_qpolynomial_div_pow(div
, 1);
1641 __isl_give isl_qpolynomial
*isl_qpolynomial_rat_cst(__isl_take isl_dim
*dim
,
1642 const isl_int n
, const isl_int d
)
1644 struct isl_qpolynomial
*qp
;
1645 struct isl_upoly_cst
*cst
;
1647 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1651 cst
= isl_upoly_as_cst(qp
->upoly
);
1652 isl_int_set(cst
->n
, n
);
1653 isl_int_set(cst
->d
, d
);
1658 static int up_set_active(__isl_keep
struct isl_upoly
*up
, int *active
, int d
)
1660 struct isl_upoly_rec
*rec
;
1666 if (isl_upoly_is_cst(up
))
1670 active
[up
->var
] = 1;
1672 rec
= isl_upoly_as_rec(up
);
1673 for (i
= 0; i
< rec
->n
; ++i
)
1674 if (up_set_active(rec
->p
[i
], active
, d
) < 0)
1680 static int set_active(__isl_keep isl_qpolynomial
*qp
, int *active
)
1683 int d
= isl_dim_total(qp
->dim
);
1688 for (i
= 0; i
< d
; ++i
)
1689 for (j
= 0; j
< qp
->div
->n_row
; ++j
) {
1690 if (isl_int_is_zero(qp
->div
->row
[j
][2 + i
]))
1696 return up_set_active(qp
->upoly
, active
, d
);
1699 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial
*qp
,
1700 enum isl_dim_type type
, unsigned first
, unsigned n
)
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
, return -1);
1716 active
= isl_calloc_array(set
->ctx
, int, isl_dim_total(qp
->dim
));
1717 if (set_active(qp
, active
) < 0)
1720 if (type
== isl_dim_set
)
1721 first
+= isl_dim_size(qp
->dim
, isl_dim_param
);
1722 for (i
= 0; i
< n
; ++i
)
1723 if (active
[first
+ i
]) {
1736 __isl_give
struct isl_upoly
*isl_upoly_drop(__isl_take
struct isl_upoly
*up
,
1737 unsigned first
, unsigned n
)
1740 struct isl_upoly_rec
*rec
;
1744 if (n
== 0 || up
->var
< 0 || up
->var
< first
)
1746 if (up
->var
< first
+ n
) {
1747 up
= replace_by_constant_term(up
);
1748 return isl_upoly_drop(up
, first
, n
);
1750 up
= isl_upoly_cow(up
);
1754 rec
= isl_upoly_as_rec(up
);
1758 for (i
= 0; i
< rec
->n
; ++i
) {
1759 rec
->p
[i
] = isl_upoly_drop(rec
->p
[i
], first
, n
);
1770 __isl_give isl_qpolynomial
*isl_qpolynomial_drop_dims(
1771 __isl_take isl_qpolynomial
*qp
,
1772 enum isl_dim_type type
, unsigned first
, unsigned n
)
1776 if (n
== 0 && !isl_dim_get_tuple_name(qp
->dim
, type
))
1779 qp
= isl_qpolynomial_cow(qp
);
1783 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
1785 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
1786 type
== isl_dim_set
, goto error
);
1788 qp
->dim
= isl_dim_drop(qp
->dim
, type
, first
, n
);
1792 if (type
== isl_dim_set
)
1793 first
+= isl_dim_size(qp
->dim
, isl_dim_param
);
1795 qp
->div
= isl_mat_drop_cols(qp
->div
, 2 + first
, n
);
1799 qp
->upoly
= isl_upoly_drop(qp
->upoly
, first
, n
);
1805 isl_qpolynomial_free(qp
);
1810 #define PW isl_pw_qpolynomial
1812 #define EL isl_qpolynomial
1814 #define IS_ZERO is_zero
1818 #define ADD(d,e1,e2) isl_qpolynomial_add(e1,e2)
1820 #include <isl_pw_templ.c>
1823 #define UNION isl_union_pw_qpolynomial
1825 #define PART isl_pw_qpolynomial
1827 #define PARTS pw_qpolynomial
1829 #include <isl_union_templ.c>
1831 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial
*pwqp
)
1839 if (!isl_set_fast_is_universe(pwqp
->p
[0].set
))
1842 return isl_qpolynomial_is_one(pwqp
->p
[0].qp
);
1845 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_mul(
1846 __isl_take isl_pw_qpolynomial
*pwqp1
,
1847 __isl_take isl_pw_qpolynomial
*pwqp2
)
1850 struct isl_pw_qpolynomial
*res
;
1853 if (!pwqp1
|| !pwqp2
)
1856 isl_assert(pwqp1
->dim
->ctx
, isl_dim_equal(pwqp1
->dim
, pwqp2
->dim
),
1859 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
1860 isl_pw_qpolynomial_free(pwqp2
);
1864 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
1865 isl_pw_qpolynomial_free(pwqp1
);
1869 if (isl_pw_qpolynomial_is_one(pwqp1
)) {
1870 isl_pw_qpolynomial_free(pwqp1
);
1874 if (isl_pw_qpolynomial_is_one(pwqp2
)) {
1875 isl_pw_qpolynomial_free(pwqp2
);
1879 n
= pwqp1
->n
* pwqp2
->n
;
1880 res
= isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1
->dim
), n
);
1882 for (i
= 0; i
< pwqp1
->n
; ++i
) {
1883 for (j
= 0; j
< pwqp2
->n
; ++j
) {
1884 struct isl_set
*common
;
1885 struct isl_qpolynomial
*prod
;
1886 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
1887 isl_set_copy(pwqp2
->p
[j
].set
));
1888 if (isl_set_fast_is_empty(common
)) {
1889 isl_set_free(common
);
1893 prod
= isl_qpolynomial_mul(
1894 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
1895 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
1897 res
= isl_pw_qpolynomial_add_piece(res
, common
, prod
);
1901 isl_pw_qpolynomial_free(pwqp1
);
1902 isl_pw_qpolynomial_free(pwqp2
);
1906 isl_pw_qpolynomial_free(pwqp1
);
1907 isl_pw_qpolynomial_free(pwqp2
);
1911 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_neg(
1912 __isl_take isl_pw_qpolynomial
*pwqp
)
1915 struct isl_pw_qpolynomial
*res
;
1921 if (isl_pw_qpolynomial_is_zero(pwqp
))
1924 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
1926 for (i
= 0; i
< pwqp
->n
; ++i
) {
1927 pwqp
->p
[i
].qp
= isl_qpolynomial_neg(pwqp
->p
[i
].qp
);
1934 isl_pw_qpolynomial_free(pwqp
);
1938 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_sub(
1939 __isl_take isl_pw_qpolynomial
*pwqp1
,
1940 __isl_take isl_pw_qpolynomial
*pwqp2
)
1942 return isl_pw_qpolynomial_add(pwqp1
, isl_pw_qpolynomial_neg(pwqp2
));
1945 __isl_give
struct isl_upoly
*isl_upoly_eval(
1946 __isl_take
struct isl_upoly
*up
, __isl_take isl_vec
*vec
)
1949 struct isl_upoly_rec
*rec
;
1950 struct isl_upoly
*res
;
1951 struct isl_upoly
*base
;
1953 if (isl_upoly_is_cst(up
)) {
1958 rec
= isl_upoly_as_rec(up
);
1962 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1964 base
= isl_upoly_rat_cst(up
->ctx
, vec
->el
[1 + up
->var
], vec
->el
[0]);
1966 res
= isl_upoly_eval(isl_upoly_copy(rec
->p
[rec
->n
- 1]),
1969 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1970 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1971 res
= isl_upoly_sum(res
,
1972 isl_upoly_eval(isl_upoly_copy(rec
->p
[i
]),
1973 isl_vec_copy(vec
)));
1976 isl_upoly_free(base
);
1986 __isl_give isl_qpolynomial
*isl_qpolynomial_eval(
1987 __isl_take isl_qpolynomial
*qp
, __isl_take isl_point
*pnt
)
1990 struct isl_upoly
*up
;
1995 isl_assert(pnt
->dim
->ctx
, isl_dim_equal(pnt
->dim
, qp
->dim
), goto error
);
1997 if (qp
->div
->n_row
== 0)
1998 ext
= isl_vec_copy(pnt
->vec
);
2001 unsigned dim
= isl_dim_total(qp
->dim
);
2002 ext
= isl_vec_alloc(qp
->dim
->ctx
, 1 + dim
+ qp
->div
->n_row
);
2006 isl_seq_cpy(ext
->el
, pnt
->vec
->el
, pnt
->vec
->size
);
2007 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
2008 isl_seq_inner_product(qp
->div
->row
[i
] + 1, ext
->el
,
2009 1 + dim
+ i
, &ext
->el
[1+dim
+i
]);
2010 isl_int_fdiv_q(ext
->el
[1+dim
+i
], ext
->el
[1+dim
+i
],
2011 qp
->div
->row
[i
][0]);
2015 up
= isl_upoly_eval(isl_upoly_copy(qp
->upoly
), ext
);
2019 dim
= isl_dim_copy(qp
->dim
);
2020 isl_qpolynomial_free(qp
);
2021 isl_point_free(pnt
);
2023 return isl_qpolynomial_alloc(dim
, 0, up
);
2025 isl_qpolynomial_free(qp
);
2026 isl_point_free(pnt
);
2030 int isl_upoly_cmp(__isl_keep
struct isl_upoly_cst
*cst1
,
2031 __isl_keep
struct isl_upoly_cst
*cst2
)
2036 isl_int_mul(t
, cst1
->n
, cst2
->d
);
2037 isl_int_submul(t
, cst2
->n
, cst1
->d
);
2038 cmp
= isl_int_sgn(t
);
2043 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial
*qp1
,
2044 __isl_keep isl_qpolynomial
*qp2
)
2046 struct isl_upoly_cst
*cst1
, *cst2
;
2050 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), return -1);
2051 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), return -1);
2052 if (isl_qpolynomial_is_nan(qp1
))
2054 if (isl_qpolynomial_is_nan(qp2
))
2056 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2057 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2059 return isl_upoly_cmp(cst1
, cst2
) <= 0;
2062 __isl_give isl_qpolynomial
*isl_qpolynomial_min_cst(
2063 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
2065 struct isl_upoly_cst
*cst1
, *cst2
;
2070 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
2071 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
2072 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2073 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2074 cmp
= isl_upoly_cmp(cst1
, cst2
);
2077 isl_qpolynomial_free(qp2
);
2079 isl_qpolynomial_free(qp1
);
2084 isl_qpolynomial_free(qp1
);
2085 isl_qpolynomial_free(qp2
);
2089 __isl_give isl_qpolynomial
*isl_qpolynomial_max_cst(
2090 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
2092 struct isl_upoly_cst
*cst1
, *cst2
;
2097 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
2098 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
2099 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2100 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2101 cmp
= isl_upoly_cmp(cst1
, cst2
);
2104 isl_qpolynomial_free(qp2
);
2106 isl_qpolynomial_free(qp1
);
2111 isl_qpolynomial_free(qp1
);
2112 isl_qpolynomial_free(qp2
);
2116 __isl_give isl_qpolynomial
*isl_qpolynomial_insert_dims(
2117 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
,
2118 unsigned first
, unsigned n
)
2127 qp
= isl_qpolynomial_cow(qp
);
2131 isl_assert(qp
->div
->ctx
, first
<= isl_dim_size(qp
->dim
, type
),
2134 g_pos
= pos(qp
->dim
, type
) + first
;
2136 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + g_pos
, n
);
2140 total
= qp
->div
->n_col
- 2;
2141 if (total
> g_pos
) {
2143 exp
= isl_alloc_array(qp
->div
->ctx
, int, total
- g_pos
);
2146 for (i
= 0; i
< total
- g_pos
; ++i
)
2148 qp
->upoly
= expand(qp
->upoly
, exp
, g_pos
);
2154 qp
->dim
= isl_dim_insert(qp
->dim
, type
, first
, n
);
2160 isl_qpolynomial_free(qp
);
2164 __isl_give isl_qpolynomial
*isl_qpolynomial_add_dims(
2165 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
, unsigned n
)
2169 pos
= isl_qpolynomial_dim(qp
, type
);
2171 return isl_qpolynomial_insert_dims(qp
, type
, pos
, n
);
2174 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_insert_dims(
2175 __isl_take isl_pw_qpolynomial
*pwqp
, enum isl_dim_type type
,
2176 unsigned first
, unsigned n
)
2183 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
2187 pwqp
->dim
= isl_dim_insert(pwqp
->dim
, type
, first
, n
);
2191 for (i
= 0; i
< pwqp
->n
; ++i
) {
2192 pwqp
->p
[i
].set
= isl_set_insert(pwqp
->p
[i
].set
, type
, first
, n
);
2193 if (!pwqp
->p
[i
].set
)
2195 pwqp
->p
[i
].qp
= isl_qpolynomial_insert_dims(pwqp
->p
[i
].qp
,
2203 isl_pw_qpolynomial_free(pwqp
);
2207 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add_dims(
2208 __isl_take isl_pw_qpolynomial
*pwqp
,
2209 enum isl_dim_type type
, unsigned n
)
2213 pos
= isl_pw_qpolynomial_dim(pwqp
, type
);
2215 return isl_pw_qpolynomial_insert_dims(pwqp
, type
, pos
, n
);
2218 static int *reordering_move(isl_ctx
*ctx
,
2219 unsigned len
, unsigned dst
, unsigned src
, unsigned n
)
2224 reordering
= isl_alloc_array(ctx
, int, len
);
2229 for (i
= 0; i
< dst
; ++i
)
2231 for (i
= 0; i
< n
; ++i
)
2232 reordering
[src
+ i
] = dst
+ i
;
2233 for (i
= 0; i
< src
- dst
; ++i
)
2234 reordering
[dst
+ i
] = dst
+ n
+ i
;
2235 for (i
= 0; i
< len
- src
- n
; ++i
)
2236 reordering
[src
+ n
+ i
] = src
+ n
+ i
;
2238 for (i
= 0; i
< src
; ++i
)
2240 for (i
= 0; i
< n
; ++i
)
2241 reordering
[src
+ i
] = dst
+ i
;
2242 for (i
= 0; i
< dst
- src
; ++i
)
2243 reordering
[src
+ n
+ i
] = src
+ i
;
2244 for (i
= 0; i
< len
- dst
- n
; ++i
)
2245 reordering
[dst
+ n
+ i
] = dst
+ n
+ i
;
2251 __isl_give isl_qpolynomial
*isl_qpolynomial_move_dims(
2252 __isl_take isl_qpolynomial
*qp
,
2253 enum isl_dim_type dst_type
, unsigned dst_pos
,
2254 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
2260 qp
= isl_qpolynomial_cow(qp
);
2264 isl_assert(qp
->dim
->ctx
, src_pos
+ n
<= isl_dim_size(qp
->dim
, src_type
),
2267 g_dst_pos
= pos(qp
->dim
, dst_type
) + dst_pos
;
2268 g_src_pos
= pos(qp
->dim
, src_type
) + src_pos
;
2269 if (dst_type
> src_type
)
2272 qp
->div
= isl_mat_move_cols(qp
->div
, 2 + g_dst_pos
, 2 + g_src_pos
, n
);
2279 reordering
= reordering_move(qp
->dim
->ctx
,
2280 qp
->div
->n_col
- 2, g_dst_pos
, g_src_pos
, n
);
2284 qp
->upoly
= reorder(qp
->upoly
, reordering
);
2289 qp
->dim
= isl_dim_move(qp
->dim
, dst_type
, dst_pos
, src_type
, src_pos
, n
);
2295 isl_qpolynomial_free(qp
);
2299 __isl_give
struct isl_upoly
*isl_upoly_from_affine(isl_ctx
*ctx
, isl_int
*f
,
2300 isl_int denom
, unsigned len
)
2303 struct isl_upoly
*up
;
2305 isl_assert(ctx
, len
>= 1, return NULL
);
2307 up
= isl_upoly_rat_cst(ctx
, f
[0], denom
);
2308 for (i
= 0; i
< len
- 1; ++i
) {
2309 struct isl_upoly
*t
;
2310 struct isl_upoly
*c
;
2312 if (isl_int_is_zero(f
[1 + i
]))
2315 c
= isl_upoly_rat_cst(ctx
, f
[1 + i
], denom
);
2316 t
= isl_upoly_pow(ctx
, i
, 1);
2317 t
= isl_upoly_mul(c
, t
);
2318 up
= isl_upoly_sum(up
, t
);
2324 __isl_give isl_qpolynomial
*isl_qpolynomial_from_affine(__isl_take isl_dim
*dim
,
2325 isl_int
*f
, isl_int denom
)
2327 struct isl_upoly
*up
;
2332 up
= isl_upoly_from_affine(dim
->ctx
, f
, denom
, 1 + isl_dim_total(dim
));
2334 return isl_qpolynomial_alloc(dim
, 0, up
);
2337 __isl_give isl_qpolynomial
*isl_qpolynomial_from_constraint(
2338 __isl_take isl_constraint
*c
, enum isl_dim_type type
, unsigned pos
)
2342 struct isl_upoly
*up
;
2343 isl_qpolynomial
*qp
;
2349 isl_int_init(denom
);
2351 isl_constraint_get_coefficient(c
, type
, pos
, &denom
);
2352 isl_constraint_set_coefficient(c
, type
, pos
, c
->ctx
->zero
);
2353 sgn
= isl_int_sgn(denom
);
2354 isl_int_abs(denom
, denom
);
2355 up
= isl_upoly_from_affine(c
->ctx
, c
->line
[0], denom
,
2356 1 + isl_constraint_dim(c
, isl_dim_all
));
2358 isl_int_neg(denom
, denom
);
2359 isl_constraint_set_coefficient(c
, type
, pos
, denom
);
2361 dim
= isl_dim_copy(c
->bmap
->dim
);
2363 isl_int_clear(denom
);
2364 isl_constraint_free(c
);
2366 qp
= isl_qpolynomial_alloc(dim
, 0, up
);
2368 qp
= isl_qpolynomial_neg(qp
);
2372 __isl_give
struct isl_upoly
*isl_upoly_subs(__isl_take
struct isl_upoly
*up
,
2373 unsigned first
, unsigned n
, __isl_keep
struct isl_upoly
**subs
)
2376 struct isl_upoly_rec
*rec
;
2377 struct isl_upoly
*base
, *res
;
2382 if (isl_upoly_is_cst(up
))
2385 if (up
->var
< first
)
2388 rec
= isl_upoly_as_rec(up
);
2392 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
2394 if (up
->var
>= first
+ n
)
2395 base
= isl_upoly_pow(up
->ctx
, up
->var
, 1);
2397 base
= isl_upoly_copy(subs
[up
->var
- first
]);
2399 res
= isl_upoly_subs(isl_upoly_copy(rec
->p
[rec
->n
- 1]), first
, n
, subs
);
2400 for (i
= rec
->n
- 2; i
>= 0; --i
) {
2401 struct isl_upoly
*t
;
2402 t
= isl_upoly_subs(isl_upoly_copy(rec
->p
[i
]), first
, n
, subs
);
2403 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
2404 res
= isl_upoly_sum(res
, t
);
2407 isl_upoly_free(base
);
2416 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
2417 * in "qp" by subs[i].
2419 __isl_give isl_qpolynomial
*isl_qpolynomial_substitute(
2420 __isl_take isl_qpolynomial
*qp
,
2421 enum isl_dim_type type
, unsigned first
, unsigned n
,
2422 __isl_keep isl_qpolynomial
**subs
)
2425 struct isl_upoly
**ups
;
2430 qp
= isl_qpolynomial_cow(qp
);
2433 for (i
= 0; i
< n
; ++i
)
2437 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_dim_size(qp
->dim
, type
),
2440 for (i
= 0; i
< n
; ++i
)
2441 isl_assert(qp
->dim
->ctx
, isl_dim_equal(qp
->dim
, subs
[i
]->dim
),
2444 isl_assert(qp
->dim
->ctx
, qp
->div
->n_row
== 0, goto error
);
2445 for (i
= 0; i
< n
; ++i
)
2446 isl_assert(qp
->dim
->ctx
, subs
[i
]->div
->n_row
== 0, goto error
);
2448 first
+= pos(qp
->dim
, type
);
2450 ups
= isl_alloc_array(qp
->dim
->ctx
, struct isl_upoly
*, n
);
2453 for (i
= 0; i
< n
; ++i
)
2454 ups
[i
] = subs
[i
]->upoly
;
2456 qp
->upoly
= isl_upoly_subs(qp
->upoly
, first
, n
, ups
);
2465 isl_qpolynomial_free(qp
);
2469 __isl_give isl_basic_set
*add_div_constraints(__isl_take isl_basic_set
*bset
,
2470 __isl_take isl_mat
*div
)
2478 bset
= isl_basic_set_extend_constraints(bset
, 0, 2 * div
->n_row
);
2481 total
= isl_basic_set_total_dim(bset
);
2482 for (i
= 0; i
< div
->n_row
; ++i
)
2483 if (isl_basic_set_add_div_constraints_var(bset
,
2484 total
- div
->n_row
+ i
, div
->row
[i
]) < 0)
2491 isl_basic_set_free(bset
);
2495 /* Extend "bset" with extra set dimensions for each integer division
2496 * in "qp" and then call "fn" with the extended bset and the polynomial
2497 * that results from replacing each of the integer divisions by the
2498 * corresponding extra set dimension.
2500 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial
*qp
,
2501 __isl_keep isl_basic_set
*bset
,
2502 int (*fn
)(__isl_take isl_basic_set
*bset
,
2503 __isl_take isl_qpolynomial
*poly
, void *user
), void *user
)
2507 isl_qpolynomial
*poly
;
2511 if (qp
->div
->n_row
== 0)
2512 return fn(isl_basic_set_copy(bset
), isl_qpolynomial_copy(qp
),
2515 div
= isl_mat_copy(qp
->div
);
2516 dim
= isl_dim_copy(qp
->dim
);
2517 dim
= isl_dim_add(dim
, isl_dim_set
, qp
->div
->n_row
);
2518 poly
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_copy(qp
->upoly
));
2519 bset
= isl_basic_set_copy(bset
);
2520 bset
= isl_basic_set_add(bset
, isl_dim_set
, qp
->div
->n_row
);
2521 bset
= add_div_constraints(bset
, div
);
2523 return fn(bset
, poly
, user
);
2528 /* Return total degree in variables first (inclusive) up to last (exclusive).
2530 int isl_upoly_degree(__isl_keep
struct isl_upoly
*up
, int first
, int last
)
2534 struct isl_upoly_rec
*rec
;
2538 if (isl_upoly_is_zero(up
))
2540 if (isl_upoly_is_cst(up
) || up
->var
< first
)
2543 rec
= isl_upoly_as_rec(up
);
2547 for (i
= 0; i
< rec
->n
; ++i
) {
2550 if (isl_upoly_is_zero(rec
->p
[i
]))
2552 d
= isl_upoly_degree(rec
->p
[i
], first
, last
);
2562 /* Return total degree in set variables.
2564 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial
*poly
)
2572 ovar
= isl_dim_offset(poly
->dim
, isl_dim_set
);
2573 nvar
= isl_dim_size(poly
->dim
, isl_dim_set
);
2574 return isl_upoly_degree(poly
->upoly
, ovar
, ovar
+ nvar
);
2577 __isl_give
struct isl_upoly
*isl_upoly_coeff(__isl_keep
struct isl_upoly
*up
,
2578 unsigned pos
, int deg
)
2581 struct isl_upoly_rec
*rec
;
2586 if (isl_upoly_is_cst(up
) || up
->var
< pos
) {
2588 return isl_upoly_copy(up
);
2590 return isl_upoly_zero(up
->ctx
);
2593 rec
= isl_upoly_as_rec(up
);
2597 if (up
->var
== pos
) {
2599 return isl_upoly_copy(rec
->p
[deg
]);
2601 return isl_upoly_zero(up
->ctx
);
2604 up
= isl_upoly_copy(up
);
2605 up
= isl_upoly_cow(up
);
2606 rec
= isl_upoly_as_rec(up
);
2610 for (i
= 0; i
< rec
->n
; ++i
) {
2611 struct isl_upoly
*t
;
2612 t
= isl_upoly_coeff(rec
->p
[i
], pos
, deg
);
2615 isl_upoly_free(rec
->p
[i
]);
2625 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
2627 __isl_give isl_qpolynomial
*isl_qpolynomial_coeff(
2628 __isl_keep isl_qpolynomial
*qp
,
2629 enum isl_dim_type type
, unsigned t_pos
, int deg
)
2632 struct isl_upoly
*up
;
2638 isl_assert(qp
->div
->ctx
, t_pos
< isl_dim_size(qp
->dim
, type
),
2641 g_pos
= pos(qp
->dim
, type
) + t_pos
;
2642 up
= isl_upoly_coeff(qp
->upoly
, g_pos
, deg
);
2644 c
= isl_qpolynomial_alloc(isl_dim_copy(qp
->dim
), qp
->div
->n_row
, up
);
2647 isl_mat_free(c
->div
);
2648 c
->div
= isl_mat_copy(qp
->div
);
2653 isl_qpolynomial_free(c
);
2657 /* Homogenize the polynomial in the variables first (inclusive) up to
2658 * last (exclusive) by inserting powers of variable first.
2659 * Variable first is assumed not to appear in the input.
2661 __isl_give
struct isl_upoly
*isl_upoly_homogenize(
2662 __isl_take
struct isl_upoly
*up
, int deg
, int target
,
2663 int first
, int last
)
2666 struct isl_upoly_rec
*rec
;
2670 if (isl_upoly_is_zero(up
))
2674 if (isl_upoly_is_cst(up
) || up
->var
< first
) {
2675 struct isl_upoly
*hom
;
2677 hom
= isl_upoly_pow(up
->ctx
, first
, target
- deg
);
2680 rec
= isl_upoly_as_rec(hom
);
2681 rec
->p
[target
- deg
] = isl_upoly_mul(rec
->p
[target
- deg
], up
);
2686 up
= isl_upoly_cow(up
);
2687 rec
= isl_upoly_as_rec(up
);
2691 for (i
= 0; i
< rec
->n
; ++i
) {
2692 if (isl_upoly_is_zero(rec
->p
[i
]))
2694 rec
->p
[i
] = isl_upoly_homogenize(rec
->p
[i
],
2695 up
->var
< last
? deg
+ i
: i
, target
,
2707 /* Homogenize the polynomial in the set variables by introducing
2708 * powers of an extra set variable at position 0.
2710 __isl_give isl_qpolynomial
*isl_qpolynomial_homogenize(
2711 __isl_take isl_qpolynomial
*poly
)
2715 int deg
= isl_qpolynomial_degree(poly
);
2720 poly
= isl_qpolynomial_insert_dims(poly
, isl_dim_set
, 0, 1);
2721 poly
= isl_qpolynomial_cow(poly
);
2725 ovar
= isl_dim_offset(poly
->dim
, isl_dim_set
);
2726 nvar
= isl_dim_size(poly
->dim
, isl_dim_set
);
2727 poly
->upoly
= isl_upoly_homogenize(poly
->upoly
, 0, deg
,
2734 isl_qpolynomial_free(poly
);
2738 __isl_give isl_term
*isl_term_alloc(__isl_take isl_dim
*dim
,
2739 __isl_take isl_mat
*div
)
2747 n
= isl_dim_total(dim
) + div
->n_row
;
2749 term
= isl_calloc(dim
->ctx
, struct isl_term
,
2750 sizeof(struct isl_term
) + (n
- 1) * sizeof(int));
2757 isl_int_init(term
->n
);
2758 isl_int_init(term
->d
);
2767 __isl_give isl_term
*isl_term_copy(__isl_keep isl_term
*term
)
2776 __isl_give isl_term
*isl_term_dup(__isl_keep isl_term
*term
)
2785 total
= isl_dim_total(term
->dim
) + term
->div
->n_row
;
2787 dup
= isl_term_alloc(isl_dim_copy(term
->dim
), isl_mat_copy(term
->div
));
2791 isl_int_set(dup
->n
, term
->n
);
2792 isl_int_set(dup
->d
, term
->d
);
2794 for (i
= 0; i
< total
; ++i
)
2795 dup
->pow
[i
] = term
->pow
[i
];
2800 __isl_give isl_term
*isl_term_cow(__isl_take isl_term
*term
)
2808 return isl_term_dup(term
);
2811 void isl_term_free(__isl_take isl_term
*term
)
2816 if (--term
->ref
> 0)
2819 isl_dim_free(term
->dim
);
2820 isl_mat_free(term
->div
);
2821 isl_int_clear(term
->n
);
2822 isl_int_clear(term
->d
);
2826 unsigned isl_term_dim(__isl_keep isl_term
*term
, enum isl_dim_type type
)
2834 case isl_dim_out
: return isl_dim_size(term
->dim
, type
);
2835 case isl_dim_div
: return term
->div
->n_row
;
2836 case isl_dim_all
: return isl_dim_total(term
->dim
) + term
->div
->n_row
;
2841 isl_ctx
*isl_term_get_ctx(__isl_keep isl_term
*term
)
2843 return term
? term
->dim
->ctx
: NULL
;
2846 void isl_term_get_num(__isl_keep isl_term
*term
, isl_int
*n
)
2850 isl_int_set(*n
, term
->n
);
2853 void isl_term_get_den(__isl_keep isl_term
*term
, isl_int
*d
)
2857 isl_int_set(*d
, term
->d
);
2860 int isl_term_get_exp(__isl_keep isl_term
*term
,
2861 enum isl_dim_type type
, unsigned pos
)
2866 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, type
), return -1);
2868 if (type
>= isl_dim_set
)
2869 pos
+= isl_dim_size(term
->dim
, isl_dim_param
);
2870 if (type
>= isl_dim_div
)
2871 pos
+= isl_dim_size(term
->dim
, isl_dim_set
);
2873 return term
->pow
[pos
];
2876 __isl_give isl_div
*isl_term_get_div(__isl_keep isl_term
*term
, unsigned pos
)
2878 isl_basic_map
*bmap
;
2885 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, isl_dim_div
),
2888 total
= term
->div
->n_col
- term
->div
->n_row
- 2;
2889 /* No nested divs for now */
2890 isl_assert(term
->dim
->ctx
,
2891 isl_seq_first_non_zero(term
->div
->row
[pos
] + 2 + total
,
2892 term
->div
->n_row
) == -1,
2895 bmap
= isl_basic_map_alloc_dim(isl_dim_copy(term
->dim
), 1, 0, 0);
2896 if ((k
= isl_basic_map_alloc_div(bmap
)) < 0)
2899 isl_seq_cpy(bmap
->div
[k
], term
->div
->row
[pos
], 2 + total
);
2901 return isl_basic_map_div(bmap
, k
);
2903 isl_basic_map_free(bmap
);
2907 __isl_give isl_term
*isl_upoly_foreach_term(__isl_keep
struct isl_upoly
*up
,
2908 int (*fn
)(__isl_take isl_term
*term
, void *user
),
2909 __isl_take isl_term
*term
, void *user
)
2912 struct isl_upoly_rec
*rec
;
2917 if (isl_upoly_is_zero(up
))
2920 isl_assert(up
->ctx
, !isl_upoly_is_nan(up
), goto error
);
2921 isl_assert(up
->ctx
, !isl_upoly_is_infty(up
), goto error
);
2922 isl_assert(up
->ctx
, !isl_upoly_is_neginfty(up
), goto error
);
2924 if (isl_upoly_is_cst(up
)) {
2925 struct isl_upoly_cst
*cst
;
2926 cst
= isl_upoly_as_cst(up
);
2929 term
= isl_term_cow(term
);
2932 isl_int_set(term
->n
, cst
->n
);
2933 isl_int_set(term
->d
, cst
->d
);
2934 if (fn(isl_term_copy(term
), user
) < 0)
2939 rec
= isl_upoly_as_rec(up
);
2943 for (i
= 0; i
< rec
->n
; ++i
) {
2944 term
= isl_term_cow(term
);
2947 term
->pow
[up
->var
] = i
;
2948 term
= isl_upoly_foreach_term(rec
->p
[i
], fn
, term
, user
);
2952 term
->pow
[up
->var
] = 0;
2956 isl_term_free(term
);
2960 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial
*qp
,
2961 int (*fn
)(__isl_take isl_term
*term
, void *user
), void *user
)
2968 term
= isl_term_alloc(isl_dim_copy(qp
->dim
), isl_mat_copy(qp
->div
));
2972 term
= isl_upoly_foreach_term(qp
->upoly
, fn
, term
, user
);
2974 isl_term_free(term
);
2976 return term
? 0 : -1;
2979 __isl_give isl_qpolynomial
*isl_qpolynomial_from_term(__isl_take isl_term
*term
)
2981 struct isl_upoly
*up
;
2982 isl_qpolynomial
*qp
;
2988 n
= isl_dim_total(term
->dim
) + term
->div
->n_row
;
2990 up
= isl_upoly_rat_cst(term
->dim
->ctx
, term
->n
, term
->d
);
2991 for (i
= 0; i
< n
; ++i
) {
2994 up
= isl_upoly_mul(up
,
2995 isl_upoly_pow(term
->dim
->ctx
, i
, term
->pow
[i
]));
2998 qp
= isl_qpolynomial_alloc(isl_dim_copy(term
->dim
), term
->div
->n_row
, up
);
3001 isl_mat_free(qp
->div
);
3002 qp
->div
= isl_mat_copy(term
->div
);
3006 isl_term_free(term
);
3009 isl_qpolynomial_free(qp
);
3010 isl_term_free(term
);
3014 __isl_give isl_qpolynomial
*isl_qpolynomial_lift(__isl_take isl_qpolynomial
*qp
,
3015 __isl_take isl_dim
*dim
)
3024 if (isl_dim_equal(qp
->dim
, dim
)) {
3029 qp
= isl_qpolynomial_cow(qp
);
3033 extra
= isl_dim_size(dim
, isl_dim_set
) -
3034 isl_dim_size(qp
->dim
, isl_dim_set
);
3035 total
= isl_dim_total(qp
->dim
);
3036 if (qp
->div
->n_row
) {
3039 exp
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
3042 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
3044 qp
->upoly
= expand(qp
->upoly
, exp
, total
);
3049 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + total
, extra
);
3052 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
3053 isl_seq_clr(qp
->div
->row
[i
] + 2 + total
, extra
);
3055 isl_dim_free(qp
->dim
);
3061 isl_qpolynomial_free(qp
);
3065 /* For each parameter or variable that does not appear in qp,
3066 * first eliminate the variable from all constraints and then set it to zero.
3068 static __isl_give isl_set
*fix_inactive(__isl_take isl_set
*set
,
3069 __isl_keep isl_qpolynomial
*qp
)
3080 d
= isl_dim_total(set
->dim
);
3081 active
= isl_calloc_array(set
->ctx
, int, d
);
3082 if (set_active(qp
, active
) < 0)
3085 for (i
= 0; i
< d
; ++i
)
3094 nparam
= isl_dim_size(set
->dim
, isl_dim_param
);
3095 nvar
= isl_dim_size(set
->dim
, isl_dim_set
);
3096 for (i
= 0; i
< nparam
; ++i
) {
3099 set
= isl_set_eliminate(set
, isl_dim_param
, i
, 1);
3100 set
= isl_set_fix_si(set
, isl_dim_param
, i
, 0);
3102 for (i
= 0; i
< nvar
; ++i
) {
3103 if (active
[nparam
+ i
])
3105 set
= isl_set_eliminate(set
, isl_dim_set
, i
, 1);
3106 set
= isl_set_fix_si(set
, isl_dim_set
, i
, 0);
3118 struct isl_opt_data
{
3119 isl_qpolynomial
*qp
;
3121 isl_qpolynomial
*opt
;
3125 static int opt_fn(__isl_take isl_point
*pnt
, void *user
)
3127 struct isl_opt_data
*data
= (struct isl_opt_data
*)user
;
3128 isl_qpolynomial
*val
;
3130 val
= isl_qpolynomial_eval(isl_qpolynomial_copy(data
->qp
), pnt
);
3134 } else if (data
->max
) {
3135 data
->opt
= isl_qpolynomial_max_cst(data
->opt
, val
);
3137 data
->opt
= isl_qpolynomial_min_cst(data
->opt
, val
);
3143 __isl_give isl_qpolynomial
*isl_qpolynomial_opt_on_domain(
3144 __isl_take isl_qpolynomial
*qp
, __isl_take isl_set
*set
, int max
)
3146 struct isl_opt_data data
= { NULL
, 1, NULL
, max
};
3151 if (isl_upoly_is_cst(qp
->upoly
)) {
3156 set
= fix_inactive(set
, qp
);
3159 if (isl_set_foreach_point(set
, opt_fn
, &data
) < 0)
3163 data
.opt
= isl_qpolynomial_zero(isl_qpolynomial_get_dim(qp
));
3166 isl_qpolynomial_free(qp
);
3170 isl_qpolynomial_free(qp
);
3171 isl_qpolynomial_free(data
.opt
);
3175 __isl_give isl_qpolynomial
*isl_qpolynomial_morph(__isl_take isl_qpolynomial
*qp
,
3176 __isl_take isl_morph
*morph
)
3180 struct isl_upoly
*up
;
3182 struct isl_upoly
**subs
;
3185 qp
= isl_qpolynomial_cow(qp
);
3190 isl_assert(ctx
, isl_dim_equal(qp
->dim
, morph
->dom
->dim
), goto error
);
3192 subs
= isl_calloc_array(ctx
, struct isl_upoly
*, morph
->inv
->n_row
- 1);
3196 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
)
3197 subs
[i
] = isl_upoly_from_affine(ctx
, morph
->inv
->row
[1 + i
],
3198 morph
->inv
->row
[0][0], morph
->inv
->n_col
);
3200 qp
->upoly
= isl_upoly_subs(qp
->upoly
, 0, morph
->inv
->n_row
- 1, subs
);
3202 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
)
3203 isl_upoly_free(subs
[i
]);
3206 mat
= isl_mat_diagonal(isl_mat_identity(ctx
, 1), isl_mat_copy(morph
->inv
));
3207 mat
= isl_mat_diagonal(mat
, isl_mat_identity(ctx
, qp
->div
->n_row
));
3208 qp
->div
= isl_mat_product(qp
->div
, mat
);
3209 isl_dim_free(qp
->dim
);
3210 qp
->dim
= isl_dim_copy(morph
->ran
->dim
);
3212 if (!qp
->upoly
|| !qp
->div
|| !qp
->dim
)
3215 isl_morph_free(morph
);
3219 isl_qpolynomial_free(qp
);
3220 isl_morph_free(morph
);
3224 static int neg_entry(void **entry
, void *user
)
3226 isl_pw_qpolynomial
**pwqp
= (isl_pw_qpolynomial
**)entry
;
3228 *pwqp
= isl_pw_qpolynomial_neg(*pwqp
);
3230 return *pwqp
? 0 : -1;
3233 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_neg(
3234 __isl_take isl_union_pw_qpolynomial
*upwqp
)
3236 upwqp
= isl_union_pw_qpolynomial_cow(upwqp
);
3240 if (isl_hash_table_foreach(upwqp
->dim
->ctx
, &upwqp
->table
,
3241 &neg_entry
, NULL
) < 0)
3246 isl_union_pw_qpolynomial_free(upwqp
);
3250 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_sub(
3251 __isl_take isl_union_pw_qpolynomial
*upwqp1
,
3252 __isl_take isl_union_pw_qpolynomial
*upwqp2
)
3254 return isl_union_pw_qpolynomial_add(upwqp1
,
3255 isl_union_pw_qpolynomial_neg(upwqp2
));
3258 static int mul_entry(void **entry
, void *user
)
3260 struct isl_union_pw_qpolynomial_match_bin_data
*data
= user
;
3262 struct isl_hash_table_entry
*entry2
;
3263 isl_pw_qpolynomial
*pwpq
= *entry
;
3266 hash
= isl_dim_get_hash(pwpq
->dim
);
3267 entry2
= isl_hash_table_find(data
->upwqp2
->dim
->ctx
,
3268 &data
->upwqp2
->table
,
3269 hash
, &has_dim
, pwpq
->dim
, 0);
3273 pwpq
= isl_pw_qpolynomial_copy(pwpq
);
3274 pwpq
= isl_pw_qpolynomial_mul(pwpq
,
3275 isl_pw_qpolynomial_copy(entry2
->data
));
3277 empty
= isl_pw_qpolynomial_is_zero(pwpq
);
3279 isl_pw_qpolynomial_free(pwpq
);
3283 isl_pw_qpolynomial_free(pwpq
);
3287 data
->res
= isl_union_pw_qpolynomial_add_pw_qpolynomial(data
->res
, pwpq
);
3292 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_mul(
3293 __isl_take isl_union_pw_qpolynomial
*upwqp1
,
3294 __isl_take isl_union_pw_qpolynomial
*upwqp2
)
3296 return match_bin_op(upwqp1
, upwqp2
, &mul_entry
);